Novel Insights Into Rheumatoid Arthritis Through Characterization of Concordant Changes in DNA Methylation and Gene Expression in Synovial Biopsies of Patients With Differing Numbers of Swollen Joints

In this study, we sought to characterize synovial tissue obtained from individuals with arthralgia and disease-specific auto-antibodies and patients with established rheumatoid arthritis (RA), by applying an integrative multi-omics approach where we investigated differences at the level of DNA methylation and gene expression in relation to disease pathogenesis. We performed concurrent whole-genome bisulphite sequencing and RNA-Sequencing on synovial tissue obtained from the knee and ankle from 4 auto-antibody positive arthralgia patients and thirteen RA patients. Through multi-omics factor analysis we observed that the latent factor explaining the variance in gene expression and DNA methylation was associated with Swollen Joint Count 66 (SJC66), with patients with SJC66 of 9 or more displaying separation from the rest. Interrogating these observed differences revealed activation of the immune response as well as dysregulation of cell adhesion pathways at the level of both DNA methylation and gene expression. We observed differences for 59 genes in particular at the level of both transcript expression and DNA methylation. Our results highlight the utility of genome-wide multi-omics profiling of synovial samples for improved understanding of changes associated with disease spread in arthralgia and RA patients, and point to novel candidate targets for the treatment of the disease.

Introduction

Rheumatoid arthritis (RA) is a complex, multifactorial, and chronic autoimmune disease that primarily affects the synovial tissue in joints (1). It affects about 1% of the population and manifests with significant unmet medical need (2). Investigating pathogenesis during different stages of disease, or across a spectrum of disease severity, is critical to optimize appropriate therapeutic interventions that affect disease progression (3, 4).

Diagnosis of established RA coincides with the development of painful and swollen joints, although circulating auto-antibodies can be detected up to 10 years before diagnosis (5, 6). Clinically manifested joint swelling reflects synovial tissue inflammation (synovitis), which is characterized by infiltration into the synovium of multiple immune cell types [including T cells, B cells, macrophages, plasma cells, dendritic cells, natural killer (NK) cells, and mast cells (7)], with up to 18 distinct infiltrating cell populations being reported on in a recent single cell transcriptomics analysis (8). As disease progresses, synovial fibroblasts adopt an increasingly aggressive and invasive phenotype, promoting further inflammation and joint damage together with other processes induced by the inflammatory environment, such as the differentiation of bone-resorbing osteoclasts (9, 10). Disease progression in early RA is often associated with the involvement of an increasing number of inflamed joints, but the mechanisms responsible for this spread of disease are poorly understood. Moreover, differences in the rate of disease manifestation and variability of response to therapy indicate that different pathophysiological mechanisms are implicated in disease development and progression compared to disease etiology (11).

An increasing body of evidence indicates that epigenetic modifications play an important role in the regulation of RA pathogenesis (12). Several array-based studies have reported widespread differences in DNA methylation among peripheral blood cells from RA patients and controls, suggesting that epigenetic modifications in circulating cells associate with disease (13, 14). However, wider conclusions may be limited by the unknown correlation of these effects to synovial cells directly at the site of inflammation. Epigenetic modifications have also been implicated in modulating the function of synovial fibroblasts in RA, through comparisons of DNA methylation patterns of cultured cells isolated from RA and osteoarthritis patients (15–20). Such studies have identified DNA methylation patterns that distinguish RA from other forms of arthritis, along with regulatory elements and biomarkers related to the pathological phenotype of RA. However, to identify novel candidate genes for therapeutic interventions that affect disease progression, it is important to study samples from patients across different stages of disease. Two such studies using cultured fibroblast-like synoviocytes found small but statistically significant patterns of hypomethylation in patients with longstanding RA compared to those with early RA, suggesting that the DNA methylome could be associated with the transformation of synovial fibroblasts into invasive cells capable of joint destruction and the resulting disease progression (21, 22).

Here, we gain insights into RA heterogeneity by combining whole-genome DNA methylome and transcriptome analyses from the same synovial biopsies and stratifying patients according to the number of swollen joints, a clinical parameter reflective of disease evolution. We find that swollen joint count based on 66 joints (SJC66), which reflects the amount and spread of inflamed synovial tissue, is associated with major changes in gene transcription and DNA methylation at promoters. By cross-interrogating differentially methylated genomic regions and their associated genes, we reveal novel candidate loci associated with the spread of the disease across joints.

Materials and Methods

Sample Description

Synovial tissue was collected via a mini-arthroscopic procedure as described previously (23) from patients at the Amsterdam University Medical Center, University of Amsterdam. A total of 17 samples were obtained from three cohorts. The first cohort contained individuals that had either arthralgia and/or a positive family history for RA, but without arthritis (as determined by a clinician), and that were positive for IgM Rheumatoid Factor (IgM-RF) and/or Anti-Citrullinated Protein Antibodies (ACPA) (Pre-synoviomics; n = 4) (24). The second cohort consisted of individuals that at inclusion were Disease Modifying Anti-Rheumatic Drug (DMARD)-naïve with early arthritis, as defined by a disease duration of less than 1 year (Synoviomics; n = 9) (25, 26). The third cohort contained samples from patients with established RA on active treatment with a disease duration of more than one year who had at least one swollen joint suitable for synovial tissue sampling (Synoviomics II; n = 4) (27). For all analyses, samples were prepared simultaneously to mitigate batch effects.

All subjects provided written informed consent and the collection and use of the samples received Institutional Review Board review and approval. Characteristics of patients included in this study are listed in Table S1.

RNA-Sequencing and Whole Genome Bisulphite Sequencing

Flash-frozen synovial tissue biopsies were utilized to simultaneously isolate RNA and DNA using an AllPrep DNA/RNA Mini kit (Qiagen), with QIAshredder spin columns (Qiagen) used to disrupt the tissue. RNA samples were quantified and their integrity assessed using Qubit RNA Broad Range Assay Kit (Thermo Fisher Scientific) and an Agilent 2100 Bioanalyzer RNA 6000 Nano Kit (Agilent Technologies), respectively. Depending on sample yield, DNA samples were quantified using Qubit DNA BR or Qubit DNA HS kits (Thermo Fisher Scientific).

RNA-Seq libraries were generated from 150 ng of total RNA. The TruSeq^® Stranded Total RNA LT was used with a Ribo-Zero™ Human/Mouse/Rat kit (Illumina), following the ‘Low Sample’ protocol except for two modifications. Firstly, the time of the ‘Elution 2 – Frag – Prime’ program was reduced from 8 to 6 min to increase the length of the RNA fragments. Secondly, 11 instead of 15 cycles were used to enrich the DNA fragments. Libraries were quantified with a KAPA Library Quantification Kit (KAPA Biosystems) on a QuantStudio 12 K Flex Real-Time PCR System (Thermo Fisher Scientific). Five samples were multiplexed per lane and libraries were clustered and sequenced using HiSeq^® PE Cluster Kit v3 – cBot™ and HiSeq^® SBS Kit v3 kits (Illumina). Paired-end sequencing (2 x 76 bp) was performed using a HiSeq 1500 (Illumina) to a depth of ~40 M reads per sample.

Whole Genome Bisulphite Sequencing (WGBS) libraries were generated using an EpiGnome Methyl-Seq Library Preparation Kit (Epicentre, now Illumina) from 100 ng of sample DNA. Bisulphite conversion was performed using an EZ DNA Methylation-Lightning Kit (Zymo Research). Each bisulphite conversion reaction contained 500 pg of unmethylated lambda DNA (Promega), which was used as a control to verify that bisulphite conversion efficiencies were at least 98%. Libraries were quantified using a KAPA Library Quantification Kit (KAPA Biosystems) on a 7900HT Real-Time PCR System (Thermo Fisher Scientific). Six samples were multiplexed per lane and libraries were clustered and sequenced using HiSeq^® PE Cluster Kit v4 – cBot™ and HiSeq^® SBS Kit v4 kits (Illumina). Paired-end sequencing (2 x 125 bp) was performed using a HiSeq 1500 (Illumina) to a depth of ~600M reads per sample. To provide sufficient coverage each batch was sequenced over 2 high output runs.

Exploratory Analyses of DNA Methylation and Gene Expression

To concurrently explore the DNA methylation and gene expression data, Multi-Omics Factor Analysis (MOFA) was applied (v1.0.0) (28). In short, MOFA performs unsupervised dimensionality reduction simultaneously across multiple data modalities from the same sample through a small number of inferred latent factors, enabling the detection of co-ordinated changes between the different data modalities (28). Here, we used the 5,000 most variable CpGs and genes as input and with the number of latent factors set to 9, the tolerance to 0.1 and the factor threshold to 0.02.

Gene Expression Data Analysis

Quality assessment of the raw reads was performed using FastQC (v0.11.2) (29). Adapters were removed from reads using Scythe (v0.991) (30) and sequences were quality-trimmed using Sickle (v1.33) (31) using a quality threshold of 20. Alignment to the GRCh38 build of the human genome was performed using Kallisto (v0.44.0) (32).

All subsequent analyses were performed in R (v3.5.0) (33). Gene-level counts were generated from the transcript abundances using tximport (v1.12.0) (34). Allosome-associated genes were removed to mitigate obvious sex effects. Differential Gene Expression (DGE) analyses were conducted using DESeq2 (v1.22.2) (35) where SJC66_high was compared with SJC66_low whilst correcting for sex, age and DMARD usage using the following design formula:

$$GeneExpression\sim sex+age+DMARDusage+SJC{66}_{dichotomized}$$

DMARD usage was a binary variable defined by the usage of any medication: conventional DMARDs (cDMARD, including methotrexate) or biological DMARDs (bDMARD). As Kallisto provided abundance levels for individual transcripts, Gene Differential Expression (GDE) analyses were also conducted to identify genes where particular transcripts were differentially expressed (36). In short, differential expression analyses was performed using DESeq2 and the resulting p-values were combined using the Lancaster aggregation method found in aggregation (v1.0.1) (36), where observations were weighted by the base expression. DGEs and GDEs were defined as genes with a false discovery rate (FDR)-adjusted p-value less than 0.05.

DNA Methylation Data Analysis

Quality assessment of the raw reads was performed using FastQC (v0.11.2) (29). Adapter and quality trimming was performed using Skewer (v0.1.123) (37) and a quality filter of 20. To assess bisulphite conversion rates, Bismark (v0.14.1) (38) was used to align the reads to the genome of the phage lambda, and again for alignment to the GRCh38 build of the human genome. Post-alignment filtering of unmapped reads, reads aligning at multiple locations and reads with a mapping score lower than 10 was carried out using SAMtools (39).

All subsequent analyses were performed in R (v3.5.0). CpG loci located on the allosomes were removed to mitigate the sex effect. The differential methylation analyses were performed using dmrseq (v1.2.5) (40), where we contrasted SJC66_high with SJC66_low whilst correcting for sex, age and DMARD usage using the following design formula:

$$Methylation\sim sex+age+DMARDusage+SJC{66}_{dichotomized}$$

Differentially Methylated Regions (DMRs) were annotated using ChipPeakAnno (v3.16.1) (41) to genes if the DMR was located within 5,000 bp upstream or 1,000 bp downstream of the gene as obtained from Gencode (v29) (42).

Integrated DNA Methylation and Gene Expression Analysis

Integrated analyses were based on the DMRs and GDEs found through the separate DNA methylation and gene expression analyses. The overlap between DMRs and GDEs were called Genes displaying both Differential Expression and Methylation (GDEMs). For each GDEM the median percentage methylation was calculated for all constituent CpGs per sample and correlated with the log₂ transformed expression counts to obtain the Pearson correlation coefficient. Confidence intervals (95%) were calculated through 1000 bootstraps for each GDEM. In short, 17 samples were randomly drawn from the original samples with replacement, whereupon the Pearson correlation coefficient was calculated. This process was repeated 1000 times to generate the empirical distribution function, which was then used to estimate the confidence intervals (43). The aforementioned bootstrapping approach was performed using the boot (v1.3) package (44). For inferential purposes, p-values were calculated by means of a permutation approach. In short, per GDEM, 1000 sets of consecutive CpGs equal to the length of the observed DMR were sampled and correlated with the gene expression signal as described above, after which the proportion of correlation coefficients higher than the observed correlation coefficient was calculated yielding the p-value.

Functional Enrichment, Cell Type Enrichment and Protein-Protein Interaction Network Analyses

Gene set enrichment analyses were performed using fgsea (v1.8.0) (45) using the Metabase pathways terms (46) as reference. Metabase pathway terms with an FDR-adjusted p-value less than 0.05 were considered significant.

Cell proportions were imputed using xCell (v1.1.0) (47), where transcripts per million were used to estimate the proportion of each of the 64 immune and stromal cell types. Subsequent linear regressions were performed to calculate the p-values to assess statistical significance. Again, we compared SJC66_high with SJC66_low also correcting for age, sex and DMARD usage.

Protein-protein interaction (PPI) network analyses were performed using the STRING (v11) database (48), in order to identify whether a set of genes was over-represented for interactions. In short, the PPI analysis returned networks of genes where the encoded proteins interacted, co-expressed or co-evolved with one another, based on text mining, curated databases, and experimental data.

Data Availability

The datasets generated and analyzed for this study can be found in the ArrayExpress repository under accession number E-MTAB-6638 and E-MTAB-6684 for WGBS and RNA-Seq, respectively. All code is hosted on GitHub at https://github.com/enricoferrero/BTCURE.

Results

Swollen Joint Count Is Associated With the Latent Factor Explaining Variance in Gene Expression and DNA Methylation

We profiled the DNA methylome and transcriptome of 17 synovial tissue samples (Table S1) using whole genome bisulphite sequencing (WGBS) and RNA-sequencing (RNA-Seq). We initially attempted to link gene expression changes to disease duration as well as to cross-patient variations in the Disease Activity Score of 28 joints (DAS28) variables, but these analyses resulted in weak and non-biologically relevant signals, and were ultimately deemed inconclusive for this set of samples (data not shown). Principal Component Analysis (PCA) indicated that variation in both DNA methylation and gene expression were independently correlated with the swollen joint count in 66 joints (SJC66; Figures 1A, B). Comparison of the first principal component of both the DNA methylation and gene expression data suggested agreement with samples 33, 19, 3, 12A and 25 broadly separating from the other samples for both modalities (Figure 1C). While sample 33 appeared to be an outlier based on the DNA methylation data the removal thereof did not alter the correlation substantially (Figure S1). To further explore DNA methylation and gene expression at a genome-wide level in an integrative fashion, we performed variance decomposition using multi-omics factor analysis (MOFA) (28). MOFA infers a set of latent factors that capture sources of variability across different measured -omic modalities of the same samples. We found that most variance (approx. 70%) was better explained by gene expression as compared to DNA methylation (Figure 1D). Further decomposition of the variance identified 8 latent factors, with LF1 explaining 40% of the variance in gene expression, whereas variance in methylation was more evenly distributed amongst all eight latent factors. Focusing on LF1, we observed a marked separation between samples with SJC66 of 9 and above compared to those with SJC66 of 8 or less (Figures 1E and S2A, B). While most samples were obtained from the knee, two were obtained from the ankle. Interrogation of the first latent factor did not indicate any correlation with the source of the sample (Figure S2C). Importantly, sex was also strongly associated with LF1 (p-value = 8.8E-03; Figure S2D) where samples with SJC66 of 9 and above were mostly males and most samples with SJC66 equal to or fewer than 8 were mostly females. Therefore, subsequent comparative analyses accounted for the imbalance in sex by removal of allosome-associated genes as well as including sex as covariate in linear models. females. Therefore, subsequent comparative analyses accounted for the imbalance in sex by removal of allosome-associated genes as well as including sex as covariate in linear models.

FIGURE 1

Figure 1 Exploratory analysis. Analysis of SJC66 regressed onto the first principal component of (A) DNA methylation and (B) gene expression annotated with the Pearson r-squared (R²) and the p-value. (C) Comparison of the first principal components for DNA methylation on the x-axis and gene expression together coloured for SJC66. Trendlines represent the mean and 95% confidence intervals. (D) Results obtained from multi-omics factor analysis, where the figure at the top depicts a bar chart representing the percentage variance explained by gene expression and DNA methylation, which has been decomposed into the separate latent factors (LF) at the bottom. (E) Decomposition of the variance observed for LF1 into the loadings assigned by MOFA to the top 20 weighted features for gene expression (top; colours represent log₂(count)) and DNA methylation (bottom; colours represent percentage methylation).

Differences in DNA Methylation and Gene Expression in Patients With High and Lower Numbers of Swollen Joints Are Associated With Immune Response and Cell Adhesion Pathways

Having observed genome-wide differences through exploratory analyses in both DNA methylation and gene expression data that associated with SJC66, we next investigated which regions and genes were differentially methylated and expressed. At this point we investigated DNA methylation and gene expression separately. While samples with SJC66 = 0 were included for the aforementioned exploratory analyses, they were excluded from subsequent comparative analyses as they were medically not a homogeneous group consisting of very early arthritis as well as late arthritis with no swelling following medication. The swollen samples were stratified according to the separation observed in the exploratory analyses (Figure 1E) and we compared SJC66_high (SJC66 ≥ 9) with SJC66_low (SJC66 < 8) correcting for age, sex and DMARD usage. Comparative methylation analysis identified 3,536 DMRs (Figure 2A and Table S2), where 2,140 were hypomethylated and 1,396 were hypermethylated. The majority of DMRs were located within 1 Kb of a TSS or in distal regions (Figure 2B). Notably, the most statistically significant DMRs spanned regions larger than 10 Kb (Table S2). For example, the top two DMRs were located within the promoter regions of the MIRLET7B and MIR3619 genes (MIRLET7BHG; 22:46,072,724-46,090,732; Figure 2C) and microRNA 10B (MIR10B; 2:176,147,225-176,162,267; Figure 2D) and exceeded 15 Kb in length. While the MIRLET7BHG-associated DMR was hypermethylated, the MIR10B-associated DMR was hypomethylated when comparing SJC66_high with SJC66_low. Functional analysis of the DMRs identified several pathways with evidence of differential methylation. Among the top 10, we observed that the NGF/TrkA MAPK pathway, immune response pathways including antigen presentation by MHC class II, macrophage and dendritic cell phenotype shifts, as well as CCR3 signalling in eosinophils, were hypermethylated (Figure 2E and Table S3).

FIGURE 2

Figure 2 SJC66_high versus SJC66_low differential methylation analysis. (A) Volcano plots depicting the -log₁₀(p-value) on the y-axis relative to the difference in % methylation when comparing the SJC66_high with SJC66_low samples on the x-axis. Colours represent the non-significant (orange) and significant regions (blue). Statistical significance is further depicted by a dashed horizontal line. (B) At the top: distribution of the DMRs relative to the transcription start site (TSS) and at the bottom: upset plot indicating the distribution of the DMRs across the different genetic features. (C, D) DMRs found at genomic co-ordinates (C) 22:46,072,724-46,090,732 and (D) 2:176,147,225-176,162,267 are located in the promoters of MIRLET7BHG and MIR10B, respectively. DNA methylation is visualized as the percentage methylation (y-axis) with a smoothed trendline per sample. (E) The 10 most enriched MetaCore gene sets associated to the promoter-bound DMRs. Tick marks represent gene ranks relative to the direction of the methylation effect and is summarized by the normalized expression score (NES).

Comparing SJC66_high with SJC66_low at the gene expression level identified 142 DGEs, of which 106 up-regulated and 36 down-regulated (Figure 3A and Table S4). The top 2 DGEs, chemokine ligand 13 (CCL13) (Figure 3B) and C-Type Lectin Domain Containing 10A (CLEC10A) (Figure 3C), were both found to be more highly expressed in the SJC66_high samples. Functional analyses revealed a striking similarity with the differential methylation results, where genes associated with antigen presentation by MHC class II as well as the NGF/TrkA MAPK pathways were overexpressed (Figure 3D and Table S5).

FIGURE 3

Figure 3 SJC66_high versus SJC66_low differential expression analysis. (A) Volcano plots depicting the -log₁₀(p-value) on the y-axis relative to the difference in log₂ fold change when comparing the SJC66_high with SJC66_low samples on the x-axis. Colours represent the non-significant genes (orange), the significant genes (blue) and the significant genes with a log₂ fold change of larger than 1 (green; “sig. interesting”). Statistical significance is further depicted by a dashed horizontal line. Mean and standard error bars superposed onto a jitterplot for the two most differentially expressed genes (B) CCL13 and (C) CLEC10A. (D) The 10 most enriched MetaCore gene sets associated to the DGEs. Tick marks represent gene ranks relative to the direction of the expression effect and are summarized by the normalized expression score (NES).

Transcripts Associated With Differentially Methylated and Expressed Genes Display Concordant Expression and Identify Nodal Points for Key Interactions

Having observed pathway-concordant differences in DNA methylation and gene expression, we were interested in identifying specific genes that were both differentially methylated and expressed. Although DNA methylation has traditionally been associated with gene-level expression, emerging evidence shows that it also regulates alternative splicing (49–51). Therefore, we complemented our Differential Gene Expression (DGE) analysis, which masks transcript-level dynamics, with Gene Differential Expression (GDE) analysis, which identifies genes that display transcript-level differences, by combining the p-values of individual transcripts associated with a single gene (36). We identified 290 genes that displayed perturbations in the expression of their transcripts (GDEs; Table S6). Combining the GDEs with the DMRs yielded 97 unique DMRs associated with 59 unique genes, which we termed Genes Differentially Expressed and Methylated (GDEMs) (Figures 4A, B). We observed that most transcripts, those associated with 52 GDEMs, were concordantly differentially expressed. Transcripts of the remaining 7 GDEMs were found to display opposite directions of expression, namely Cytohesin 1 (CYTH1), Eukaryotic Translation Initiation Factor 4 Gamma 1 (EIF4G1), a putative monooxygenase (KIAA1191), Kinesin Light Chain 1 (KLC1), Cell Division Cycle 16 (CDC16), Fas-Activated Serine/Threonine Kinase (FASTK) and G protein coupled receptor 132 (GPR132; Figure 4B). Investigation of the methylation status of the gene overlaid onto the transcript location indicated that the observed DMRs for KIAA1191, KLC1, CDC16, FASTK and GPR132 were located in the promoter shared by most transcripts (Figure 4C). The mechanisms that could cause opposite direction of expression in these genes are currently unclear although studies in Arabidopsis have identified a methylation reader complex that can enhance rather than suppress gene transcription in the presence of methylation (52). Our work therefore supports the rationale for further validation and mechanistic studies.

FIGURE 4

Figure 4 Differential methylation and alternative splicing. (A) Comparison of the -log₁₀(p-values) obtained from gene differential expression analysis and differential methylation analysis on the x- and y-axis, respectively. Grey, green, blue and purple represent the genes: unchanged for methylation and gene expression (“NS”), differentially expressed (“GDE”), differentially methylated (“DMR”), differentially methylated and expressed (“GDEM”), respectively. (B) Visualization of the Wald statistic (a statistic calculated by DESeq2 representing the effect size relative to the variance) calculated for the individual transcripts (represented in dots) associated with each of the GDEMs. Grey and blue dots represent non-significant and significantly differentially expressed transcripts, respectively. Large red circles left of the gene symbols indicate genes with transcripts that display opposite effect sizes. (C) Summary visualization of the DNA methylation and gene expression signals for GDEMs: CYTH1, EIF4G11, KIAA1191, KLC1, CDC18, FASTK and GPR132. The “Transcripts” track represents individual transcripts coloured by the Wald statistic. Stars in red indicate differentially expressed transcripts. The “M_difference” track represents the difference in percentage methylation when comparing SJC66_high with SJC66_low. The “DMR” track represents the locations of the DMRs as found by dmrseq.

Functionally, we observed that the 59 GDEMs were primarily over-represented for immune response-associated pathways, specifically T-lymphocyte associated pathways (Figure 5A). We interrogated which of the GDEMs encoded known interaction partners by querying the STRING database for documented interactions (48). Almost half (46%) of the GDEMs encoded for interacting proteins with the most connected GDEMs being ITGB2 (11 interactions) and LCP2 (9 interactions) (Figure 5C). Both ITGB2 and LCP2 were differentially methylated at multiple locations, with the largest visible differences occurring within the TSS and downstream thereof (Figures 5D, E). Transcript-wise, ITGB2-transcripts ENST00000498666 and ENST00000397852 as well as LCP2 transcript ENST00000046794 were most differentially expressed. Expression Quantitative Trait Methylation (eQTM) analyses confirmed strong inverse correlations between the differences in methylation in the promoter regions of ITGB2 (21:44918461-44921815) and LCP2 (5:170295513-170298924) with transcripts ENST00000498666 (r = -0.9; p-value = 1E-04) and ENST00000046794 (r = -0.9; p-value = 1E-04), respectively (Table 1). While the association between the ITGB2 promoter DMR and ENST00000397852 expression was non-significant (p-value = 0.2353), the correlation coefficient remained high (r = -0.87) (Table S7). Nonetheless, given the centrality of ITGB2 and LCP2 among the GDEMs would make them interesting candidates in future targeted studies.

FIGURE 5

Figure 5 Functional analyses GDEMs. (A) The 10 most overrepresented MetaCore gene sets found for the GDEMs depicting the ratio of pathway-associated GDEMs relative to the total number of pathway-associated genes. The size and colour intensity are proportional to the number of pathway-associated GDEMs and the -log₁₀(adjusted p-values), respectively. (B) GDEMs were analysed for known interactions by querying the STRING database. Represented here is the protein-protein interaction network, where the thickness of the edge is proportional to the STRINGdb evidence score and the colour proportional to the Wald statistic, respectively. (C) Bargraph representing the number of interactions per protein. Summary visualization of the DNA methylation and gene expression signals for GDEMs: (D) LCP2, (E) ITGB2. The “Transcripts” track represents individual transcripts coloured by the Wald statistic. The “M_difference” track represents the difference in percentage methylation when comparing SJC66_high with SJC66low. The “DMR” track represents the locations of the DMRs as found by dmrseq.

TABLE 1

Table 1 Differentially methylated and expressed genes.

Estimated Cellular Composition Suggests Lower Proportion of Neuronal Cells in SJC66_high

Systematic differences in DNA methylation and gene expression could reflect changes in the cellular composition. To this end, we estimated the cellular composition using the transcriptome data as input for xCell, which is capable of estimating enrichment scores of 64 immune and stromal cell types (47). By comparing the estimated proportions from SJC66_high with SJC66_low we identified significant differences for 6 cell types: neurons, dendritic cells (DCs: all, conventional and immature), megakaryocytes and platelets (Figure 6 and Table S8). Expectedly, higher enrichment scores for DCs (all subtypes) and platelets were estimated for the SJC66_high samples. By contrast, lower proportions of neuron and megakaryocyte signatures were observed for the SJC66_high samples. Notably, the difference in neuronal enrichment was found to be the most statistically significant with an almost fourfold difference when comparing SJC66_high with SJC66_low.

FIGURE 6

Figure 6 Cellular composition estimation. Cellular proportions were estimated for 64 cell types using the xCell algorithm. Of these, 6 cell types were found to be significantly different when comparing SJC66_high with SJC66_low. Visualized as scatter- and cross-bar plots are the significantly differentially represented cell types, namely the neurons, conventional dendritic cells (cDC), dendritic cells (DC), immature dendritic cells (iDC), megakaryocytes, and platelets.

Discussion

In this study, we highlight insights into RA progression by combining the outputs of parallel whole-genome DNA methylome and transcriptome analyses on extracted preparations of synovial biopsies, from auto-antibody positive individuals, early arthritis patients and patients with established RA, stratified by the number of swollen joints. We observed that synovia from patients with a higher number of swollen joints (SJC66 ≥ 9) were different at the level of DNA methylation and gene expression from synovia from patients with a lower number of swollen joints. Specifically comparing SJC66_high with SJC66_low revealed 3536 DMRs and 142 DGEs, with both datasets primarily enriched for pathways associated with immune responses. The most significant difference in methylation was found spanning the promoter regions of MIRLET7B and MIR10B. Interestingly, mouse miR-let-7b has been shown to provoke arthritic joint inflammation by remodeling naïve myeloid cells into M1 macrophages via TLR-7 ligation (53) and can augment disease severity (54). MIR10B has been shown to regulate Th17 cells in patients with ankylosing spondylitis (55) but no studies have specifically associated it with RA. At the level of transcription, CCL13 and CLEC10A were found to be the most differentially expressed. CCL13 (MCP-4) is an extensively studied chemokine that is thought to be involved with RA pathogenesis and disease progression (56–58). By contrast, not much is known about the role of CLEC10A in RA besides it being highly expressed on immature dendritic cells (DCs), monocyte-derived DCs and alternatively activated macrophages (59), as well as having been observed in the inflamed synovium of patients with active RA (60). Chemokines CCL13, CCL8, CXCL11, CXCL10, and CXCL9 regulate the recruitment of leukocytes into tissue and have therefore been implicated in the pathogenesis of RA (61). Differential methylation was observed in the vicinity of the promoter for CCL13, CXCL11, and CXCL9. Such results support a role for epigenetic/transcriptional processes in the spread of pathology to additional joints. While definitive mechanisms of joint spreading remain elusive, possible roles for immune cell migration due to chemokine expression (62, 63) could be further evaluated based on our data.

Altogether, we observed that 3% of DMRs associate with 20% of the differentially expressed GDEs. It is not surprising that not all DMRs could be linked to genes as a large number are found in distal intergenic (18.6%) or intronic (28.3%) regions, making any functional inference challenging. Of the genes that presented both differentially expression and methylation, protein-protein interaction networks indicated that half encoded for interacting proteins, suggesting that the observed GDEMs function together. The most interconnected GDEMs appeared to be ITGB2 and LCP2, with multiple regions of differential methylation observed surrounding both genes. While we observed transcript differential expression, most transcripts belonging to ITGB2 and LCP2 behaved similarly and all displayed reasonably strong inverse correlations with the DMRs located in the promoter area. ITGB2 encodes an integrin, which would typically be involved in cell-surface mediated signalling. We observed that the gene encoding ITGB2’s interaction partner integrin alpha L (ITGAL) was also differentially methylated and expressed. ITGB2 and ITGAL together form lymphocyte function-associated antigen (LFA-1), which interacts with intracellular adhesion molecule 1 (ICAM-1 or CD54) resulting in an enhanced immune cell influx into the synovial tissue (64, 65). Inhibiting LFA-1 has been reported to reduce inflammation and joint destruction in murine models of arthritis (66).

Functionally, the 59 GDEMs were primarily over-represented for immune response-associated pathways, specifically T-lymphocyte associated ones. It would be fascinating to understand why the DNA methylome and transcriptome between patients with a SJC66 of 9 and above relative to patients with a SJC66 of 8 present with such a sudden split instead of a gradual difference. It is clear that inflammation is likely an important factor contributing to the observed differences as patients with a high SJC66 also generally present higher levels of inflammation as expressed through using erythrocyte sedimentation rate (ESR) or concentration C-reactive protein (CRP). Our transcriptomic data indeed suggested an increased proportion of immune cells among SJC66_high samples, as would be expected while pathology develops and cells migrate into the affected joints. Previous work has shown that clinically manifest arthritis in established RA is associated with increased infiltration of leukocytes. In synovial tissue samples from clinically involved joints, scores for infiltration by DCs are consistently higher than in clinically uninvolved joints obtained simultaneously from the same RA patients (67). Importantly, when comparing different clinically inflamed joints from the same RA patient simultaneously, leukocyte infiltration in one inflamed joint was shown to be representative of that in other inflamed joints, supporting the notion that leukocytes migrate from one joint to another (68). Indeed, there is continuous influx of leukocytes into the joints in established RA (69). We postulate that if synovial leukocytes exhibit properties that would facilitate cell migration, arthritis might spread from one inflamed joint to another. The results presented here support a disease mechanism in which, after development of clinically established RA, inflammatory and cell adhesion-associated processes play a key role in the progression of RA to greater joint involvement (70–72). Interestingly, we also observed differences of neuronal signatures suggesting a lower relative enrichment of neurons among SJC66_high samples. In addition, enrichment analyses on the DMRs suggested hypermethylation of genes encoding nociception receptors, which are typically associated with peripheral sensory neurons (73, 74). A similar decrease in neuronal signature has previously been associated with RA severity, where the authors noted a potential role in the maladaptive response towards damage (75). This is consistent with a more general loss of anti-inflammatory control by the nervous system in RA (76). Important to note is that our observations are based on estimates made by xCell, which can only calculate enrichment scores based on signatures rather than absolute values of cells. We are therefore unable to discern whether a population increased in size or a different population had decreased. Ideally, a similar estimate would have been generated based on the DNA methylation data, but the currently available reference datasets do not include the cell types measured available in xCell.

There are two limitations of this study, namely the fact that sex confounds the separation between SJC66_high and SJC66_low and the limited sample size. While we have sought to mitigate the confounding effect of sex by removing genes and CpGs located on the allosomes as well as by including sex as a covariate in our analyses, we acknowledge that we cannot fully eliminate the possibility that a sex effect is present. Accordingly, validation studies would be necessary where the DMEGs are verified using a larger, independent cohort while controlling for an interaction effect between sex and SJC66. The observed differences in transcript expression could be validated using a quantitative PCR approach with primers designed specifically against particular transcripts. Similarly, for validating the DMRs, targeted bisulphite-sequencing using primers for the regions of interest would be a cost-effective approach.

In conclusion, our study constitutes an exploratory analysis of whole genome DNA methylation and gene expression data performed on primary synovial tissue material from auto-antibody positive arthralgia patients without arthritis as well as patients with early and established RA patients. Where previous studies investigated cells from patients with RA versus disease controls and were potentially limited by their use of cultured cells, we focus on an integrative analysis of epigenetic marks and alternative splicing associated with swelling spread, providing novel insights into the mechanisms of disease progression towards more severe phenotypes. Nonetheless, further validation is necessary if the identified target genes are to be used for monitoring or treatment of the swelling and associated inflammatory processes in the joints of RA patients.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ebi.ac.uk/arrayexpress/, E-MTAB-6638; https://www.ebi.ac.uk/arrayexpress/, E-MTAB-6684.

Ethics Statement

The studies involving human participants were reviewed and approved by AMC medisch ethisch toetsingscommissie under MEC 02/152, MEC 05/107, and MEC 07/253 for the Synoviomics, Pre-synoviomics, and Synoviomics II cohorts, respectively. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

EF, HL, DT, CL, MS, DG, PT, and RP conceived the study. KM and GR carried out the laboratory experimental work. EF and CL conceived the analytical design. EF and AL performed the data analysis. HL, DT, CL, MM, PH, and WJ helped supervise the secondary data analyses. EF, AL, KM, HL, and DT led the writing of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

AL was funded by the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. ITN-2014-EID-641665. This study was supported by the Dutch Arthritis Foundation (06-1-303 and 11-1-407), IMI BeTheCure (115142), and FP7 Euro-TEAM consortium (305549).

Conflict of Interest

AL, KM, GR, HL, DT, CL, DG, and RP were employed by GSK when this study was conducted. EF and PT were employed by Novartis and Candel Therapeutics when this study was conducted, respectively. WJ was financially supported by GSK and Mead Johnson.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank the BTCure research consortium for providing access to the samples, Victor Neduva (Target and Pathway Validation, Target Sciences, GSK) for transferring the raw data files, Erika Cule (Research Statistics, GSK) for input into statistical methodology, and Yee Ying Chang (Clinical Genetics, Genome Diagnostics Laboratory, Amsterdam UMC, University of Amsterdam) for implementing the eQTM function in R. This manuscript has been released as a Pre-Print at https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3576744 (77).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.651475/full#supplementary-material

References

1. Tak PP, Bresnihan B. The pathogenesis and prevention of joint damage in rheumatoid arthritis: Advances from synovial biopsy and tissue analysis. Arthritis Rheum (2000) 43:2619–33. doi: 10.1002/1529-0131(200012)43:12<2619::AID-ANR1>3.0.CO;2-V

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Taylor PC, Moore A, Vasilescu R, Alvir J, Tarallo M. A structured literature review of the burden of illness and unmet needs in patients with rheumatoid arthritis: a current perspective. Rheumatol Int (2016) 36:685–95. doi: 10.1007/s00296-015-3415-x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature (2014) 506:376–81. doi: 10.1038/nature12873

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Wijbrandts CA, Tak PP. Prediction of Response to Targeted Treatment in Rheumatoid Arthritis. Mayo Clin Proc (2017) 92:1129–43. doi: 10.1016/j.mayocp.2017.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Nielen MMJ, van Schaardenburg D, Reesink HW, van de Stadt RJ, van der Horst-Bruinsma IE, de Koning MHMT, et al. Specific autoantibodies precede the symptoms of rheumatoid arthritis: a study of serial measurements in blood donors. Arthritis Rheum (2004) 50:380–6. doi: 10.1002/art.20018

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Kolfenbach JR, Deane KD, Derber LA, O’Donnell CI, Gilliland WR, Edison JD, et al. Autoimmunity to peptidyl arginine deiminase type 4 precedes clinical onset of rheumatoid arthritis. Arthritis Rheum (2010) 62:2633–9. doi: 10.1002/art.27570

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Tak PP, Smeets TJ, Daha MR, Kluin PM, Meijers KA, Brand R, et al. Analysis of the synovial cell infiltrate in early rheumatoid synovial tissue in relation to local disease activity. Arthritis Rheum (1997) 40:217–25. doi: 10.1002/art.1780400206

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Zhang F, Wei K, Slowikowski K, Fonseka CY, Rao DA, Kelly S, et al. Defining inflammatory cell states in rheumatoid arthritis joint synovial tissues by integrating single-cell transcriptomics and mass cytometry. Nat Immunol (2019) 20:928–42. doi: 10.1038/s41590-019-0378-1

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Tak PP, Zvaifler NJ, Green DR, Firestein GS. Rheumatoid arthritis and p53: how oxidative stress might alter the course of inflammatory diseases. Immunol Today (2000) 21:78–82. doi: 10.1016/S0167-5699(99)01552-2

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ospelt C. Synovial fibroblasts in 2017. RMD Open (2017) 3:e000471. doi: 10.1136/rmdopen-2017-000471

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Pitzalis C, Kelly S, Humby F. New learnings on the pathophysiology of RA from synovial biopsies. Curr Opin Rheumatol (2013) 25:334–44. doi: 10.1097/BOR.0b013e32835fd8eb

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Nemtsova MV, Zaletaev DV, Bure IV, Mikhaylenko DS, Kuznetsova EB, Alekseeva EA, et al. Epigenetic Changes in the Pathogenesis of Rheumatoid Arthritis. Front Genet (2019) 10:570. doi: 10.3389/fgene.2019.00570

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ospelt C, Gay S, Klein K. Epigenetics in the pathogenesis of RA. Semin Immunopathol (2017) 39:409–19. doi: 10.1007/s00281-017-0621-5

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Nair N, Wilson AG, Barton A. DNA methylation as a marker of response in rheumatoid arthritis. Pharmacogenomics (2017) 18:1323–32. doi: 10.2217/pgs-2016-0195

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Whitaker JW, Shoemaker R, Boyle DL, Hillman J, Anderson D, Wang W, et al. An imprinted rheumatoid arthritis methylome signature reflects pathogenic phenotype. Genome Med (2013) 5:40. doi: 10.1186/gm444

PubMed Abstract | CrossRef Full Text | Google Scholar

16. de la Rica L, Urquiza JM, Gómez-Cabrero D, Islam ABMMK, López-Bigas N, Tegnér J, et al. Identification of novel markers in rheumatoid arthritis through integrated analysis of DNA methylation and microRNA expression. J Autoimmun (2013) 41:6–16. doi: 10.1016/j.jaut.2012.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Whitaker JW, Boyle DL, Bartok B, Ball ST, Gay S, Wang W, et al. Integrative omics analysis of rheumatoid arthritis identifies non-obvious therapeutic targets. PLoS One (2015) 10:1–14. doi: 10.1371/journal.pone.0124254

CrossRef Full Text | Google Scholar

18. Ekwall A-KH, Whitaker JW, Hammaker D, Bugbee WD, Wang W, Firestein GS. The Rheumatoid Arthritis Risk Gene LBH Regulates Growth in Fibroblast-like Synoviocytes. Arthritis Rheumatol (Hoboken NJ) (2015) 67:1193–202. doi: 10.1002/art.39060

CrossRef Full Text | Google Scholar

19. Hammaker D, Whitaker JW, Maeshima K, Boyle DL, Ekwall A-KH, Wang W, et al. LBH Gene Transcription Regulation by the Interplay of an Enhancer Risk Allele and DNA Methylation in Rheumatoid Arthritis. Arthritis Rheumatol (Hoboken NJ) (2016) 68:2637–45. doi: 10.1002/art.39746

CrossRef Full Text | Google Scholar

20. Ham S, Bae J-B, Lee S, Kim B-J, Han B-G, Kwok S-K, et al. Epigenetic analysis in rheumatoid arthritis synoviocytes. Exp Mol Med (2019) 51:22. doi: 10.1038/s12276-019-0215-5

CrossRef Full Text | Google Scholar

21. Ai R, Whitaker JW, Boyle DL, Tak PP, Gerlag DM, Wang W, et al. DNA Methylome Signature in Synoviocytes From Patients With Early Rheumatoid Arthritis Compared to Synoviocytes From Patients With Longstanding Rheumatoid Arthritis. Arthritis Rheumatol (2015) 67:1978–80. doi: 10.1002/art.39123

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Karouzakis E, Raza K, Kolling C, Buckley CD, Gay S, Filer A, et al. Analysis of early changes in DNA methylation in synovial fibroblasts of RA patients before diagnosis. Sci Rep (2018) 8:7370. doi: 10.1038/s41598-018-24240-2

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Gerlag DM, Tak PP. How to perform and analyse synovial biopsies. Best Pract Res Clin Rheumatol (2013) 27:195–207. doi: 10.1016/j.berh.2013.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

24. de Hair MJH, van de Sande MGH, Ramwadhdoebe TH, Hansson M, Landewé R, Van Der Leij C, et al. Features of the synovium of individuals at risk of developing rheumatoid arthritis : Implications for understanding preclinical rheumatoid arthritis. Arthritis Rheumatol (2014) 66:513–22. doi: 10.1002/art.38273

PubMed Abstract | CrossRef Full Text | Google Scholar

25. de Hair MJH, Harty LC, Gerlag DM, Pitzalis C, Veale DJ, Tak PP. Synovial tissue analysis for the discovery of diagnostic and prognostic biomarkers in patients with early arthritis. J Rheumatol (2011) 38:2068–72. doi: 10.3899/jrheum.110426

PubMed Abstract | CrossRef Full Text | Google Scholar

26. van de Sande MGH, de Hair MJH, Schuller Y, van de Sande GPM, Wijbrandts CA, Dinant HJ, et al. The features of the synovium in early rheumatoid arthritis according to the 2010 ACR/EULAR classification criteria. PLoS One (2012) 7:e36668. doi: 10.1371/journal.pone.0036668

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Aletaha D, Neogi T, Silman AJ, Funovits J, Felson DT, Bingham CO, et al. 2010 Rheumatoid arthritis classification criteria: an American College of Rheumatology/European League Against Rheumatism collaborative initiative. Ann Rheum Dis (2010) 69S:1580–8. doi: 10.1136/ard.2010.138461

CrossRef Full Text | Google Scholar

28. Argelaguet R, Velten B, Arnol D, Dietrich S, Zenz T, Marioni JC, et al. Multi-Omics Factor Analysis—a framework for unsupervised integration of multi-omics data sets. Mol Syst Biol (2018) 14:e8124. doi: 10.15252/msb.20178124

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Andrews S. FastQC: a quality control tool for high throughput sequence data. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

Google Scholar

30. Buffalo V. Scythe: a Bayesian adapter trimmer.

Google Scholar

31. Joshi N. Sickle: a windowed adaptive trimming tool for FASTQ files using quality.

Google Scholar

32. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol (2016) 34:525–7. doi: 10.1038/nbt.3519

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Development Core Team R. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing (2008). Available at: http://www.r-project.org/.

Google Scholar

34. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research (2016) 4:1521. doi: 10.12688/f1000research.7563.2

CrossRef Full Text | Google Scholar

35. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Yi L, Pimentel H, Bray NL, Pachter L. Gene-level differential analysis at transcript-level resolution. Genome Biol (2018) 19:52. doi: 10.1186/s13059-018-1419-z

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Jiang H, Lei R, Ding S-W, Zhu S. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinf (2014) 15:182. doi: 10.1186/1471-2105-15-182

CrossRef Full Text | Google Scholar

38. Krueger F, Andrews SR. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics (2011) 27:1571–2. doi: 10.1093/bioinformatics/btr167

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. Genome Project Data Processing S. The Sequence Alignment/Map format and SAMtools. Bioinformatics (2009) 25:2078–9. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Korthauer K, Chakraborty S, Benjamini Y, Irizarry RA, User G. Detection and accurate false discovery rate control of differentially methylated regions from whole genome bisulfite sequencing. Biostatistics (2018) 00:1–17. doi: 10.1093/biostatistics/kxy007

CrossRef Full Text | Google Scholar

41. Zhu LJ, Gazin C, Lawson ND, Pagès H, Lin SM, Lapointe DS, et al. ChIPpeakAnno: a Bioconductor package to annotate ChIP-seq and ChIP-chip data. BMC Bioinf (2010) 11:237. doi: 10.1186/1471-2105-11-237

CrossRef Full Text | Google Scholar

42. Frankish A, Diekhans M, Ferreira A-M, Johnson R, Jungreis I, Loveland J, et al. GENCODE reference annotation for the human and mouse genomes. Nucleic Acids Res (2019) 47:D766–73. doi: 10.1093/nar/gky955

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Efron B, Tibshirani RJ. An introduction to the bootstrap. New York, N.Y.; London: Chapman & Hall (1993).

Google Scholar

44. Canty A, Ripley B. boot. (2021).

Google Scholar

45. Sergushichev AA. An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation. bioRxiv (2016), 060012. doi: 10.1101/060012

CrossRef Full Text | Google Scholar

46. Clarivate. MetaCore.

Google Scholar

47. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res (2019) 47:D607–13. doi: 10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Shayevitch R, Askayo D, Keydar I, Ast G. The importance of DNA methylation of exons on alternative splicing. RNA (2018) 24:1351–62. doi: 10.1261/rna.064865.117

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Lev Maor G, Yearim A, Ast G. The alternative role of DNA methylation in splicing regulation. Trends Genet (2015) 31:274–80. doi: 10.1016/J.TIG.2015.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Linker SM, Urban L, Clark SJ, Chhatriwala M, Amatya S, McCarthy DJ, et al. Combined single-cell profiling of expression and DNA methylation reveals splicing regulation and heterogeneity. Genome Biol (2019) 20:30. doi: 10.1186/s13059-019-1644-0

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Harris CJ, Scheibe M, Wongpalee SP, Liu W, Cornett EM, Vaughan RM, Jacobsen SE. A DNA methylation reader complex that enhances gene transcription. Science (2018) 362(4619):1182–6. doi: 10.1126/science.aar7854

CrossRef Full Text | Google Scholar

53. Kim SJ, Chen Z, Essani AB, Elshabrawy HA, Volin MV, Volkov S, et al. Identification of a Novel Toll-like Receptor 7 Endogenous Ligand in Rheumatoid Arthritis Synovial Fluid That Can Provoke Arthritic Joint Inflammation. Arthritis Rheumatol (2016) 68:1099–110. doi: 10.1002/art.39544

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Umar S, Palasiewicz K, Van Raemdonck K, Volin MV, Romay B, Amin MA, et al. IRAK4 inhibition: a promising strategy for treating RA joint inflammation and bone erosion. Cell Mol Immunol (2020). doi: 10.1038/s41423-020-0433-8

CrossRef Full Text | Google Scholar

55. Chen L, Al-Mossawi MH, Ridley A, Sekine T, Hammitzsch A, De Wit J, et al. MiR-10b-5p is a novel Th17 regulator present in Th17 cells from ankylosing spondylitis. Ann Rheum Dis (2017) 76:620–4. doi: 10.1136/annrheumdis-2016-210175

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Hintzen C, Quaiser S, Pap T, Heinrich PC, Hermanns HM. Induction of CCL13 expression in synovial fibroblasts highlights a significant role of oncostatin M in rheumatoid arthritis. Arthritis Rheum (2009) 60:1932–43. doi: 10.1002/art.24602

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Yamaguchi A, Nozawa K, Fujishiro M, Kawasaki M, Suzuki F, Takamori K, et al. CC motif chemokine ligand 13 is associated with rheumatoid arthritis pathogenesis. Mod Rheumatol (2013) 23:856–63. doi: 10.3109/s10165-012-0752-4

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Rump L, Mattey DL, Kehoe O, Middleton J. An initial investigation into endothelial CC chemokine expression in the human rheumatoid synovium. Cytokine (2017) 97:133–40. doi: 10.1016/j.cyto.2017.05.023

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Hoober JK. ASGR1 and Its Enigmatic Relative, CLEC10A. Int J Mol Sci (2020) 21:4818. doi: 10.3390/ijms21144818

CrossRef Full Text | Google Scholar

60. McHugh J. Synovial macrophage populations linked to RA remission. Nat Rev Rheumatol (2020) 16:471. doi: 10.1038/s41584-020-0481-6

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Miyabe Y, Miyabe C, Iwai Y, Luster AD. Targeting the Chemokine System in Rheumatoid Arthritis and Vasculitis. JMA J (2020) 3:182–92. doi: 10.31662/jmaj.2020-0019

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Elemam NM, Hannawi S, Maghazachi AA. Role of Chemokines and Chemokine Receptors in Rheumatoid Arthritis. ImmunoTargets Ther (2020) 9:43–56. doi: 10.2147/itt.s243636

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Yu X, Song Z, Rao L, Tu Q, Zhou J, Yin Y, et al. Synergistic induction of CCL5, CXCL9 and CXCL10 by IFN-γ and NLRs ligands on human fibroblast-like synoviocytes—A potential immunopathological mechanism for joint inflammation in rheumatoid arthritis. Int Immunopharmacol (2020) 82:1–6. doi: 10.1016/j.intimp.2020.106356

CrossRef Full Text | Google Scholar

64. Lowin T, Straub RH. Integrins and their ligands in rheumatoid arthritis. Arthritis Res Ther (2011) 13:244. doi: 10.1186/ar3464

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Tak PP, Thurkow EW, Daha MR, Kluin PM, Smeets TJ, Meinders AE, et al. Expression of adhesion molecules in early rheumatoid synovial tissue. Clin Immunol Immunopathol (1995) 77:236–42. doi: 10.1006/clin.1995.1149

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Suchard SJ, Stetsko DK, Davis PM, Skala S, Potin D, Launay M, et al. An LFA-1 (alphaLbeta2) small-molecule antagonist reduces inflammation and joint destruction in murine models of arthritis. J Immunol (2010) 184:3917–26. doi: 10.4049/jimmunol.0901095

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Kraan MC, Versendaal H, Jonker M, Bresnihan B, Post WJ, Hart B, et al. A symptomatic synovitis precedes clinically manifest arthritis. Arthritis Rheum (1998) 41:1481–8. doi: 10.1002/1529-0131(199808)41:8<1481::AID-ART19>3.0.CO;2-O

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Kraan MC, Reece RJ, Smeets TJM, Veale DJ, Emery P, Tak PP. Comparison of synovial tissues from the knee joints and the small joints of rheumatoid arthritis patients: Implications for pathogenesis and evaluation of treatment. Arthritis Rheum (2002) 46:2034–8. doi: 10.1002/art.10556

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Thurlings RM, Wijbrandts CA, Bennink RJ, Dohmen SE, Voermans C, Wouters D, et al. Monocyte scintigraphy in rheumatoid arthritis: The dynamics of monocyte migration in immune-mediated inflammatory disease. PLoS One (2009) 4:1–6. doi: 10.1371/journal.pone.0007865

CrossRef Full Text | Google Scholar

70. Tsubaki T, Arita N, Kawakami T, Shiratsuchi T, Yamamoto H, Takubo N, et al. Characterization of histopathology and gene-expression profiles of synovitis in early rheumatoid arthritis using targeted biopsy specimens. Arthritis Res Ther (2005) 7:R825–36. doi: 10.1186/ar1751

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Lequerré T, Bansard C, Vittecoq O, Derambure C, Hiron M, Daveau M, et al. Early and long-standing rheumatoid arthritis: distinct molecular signatures identified by gene-expression profiling in synovia. Arthritis Res Ther (2006) 11:R105. doi: 10.1186/ar2744

CrossRef Full Text | Google Scholar

72. Choy E. Understanding the dynamics: pathways involved in the pathogenesis of rheumatoid arthritis. Rheumatol (Oxford) (2012) 51(Suppl 5):v3–11. doi: 10.1093/rheumatology/kes113

CrossRef Full Text | Google Scholar

73. Dubin AE, Patapoutian A. Nociceptors: The sensors of the pain pathway. J Clin Invest (2010) 120:3760–72. doi: 10.1172/JCI42843

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Miller RE, Tran PB, Obeidat AM, Raghu P, Ishihara S, Miller RJ, et al. The Role of Peripheral Nociceptive Neurons in the Pathophysiology of Osteoarthritis Pain. Curr Osteoporos Rep (2015) 13:318–26. doi: 10.1007/s11914-015-0280-1

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Orange DE, Agius P, DiCarlo EF, Robine N, Geiger H, Szymonifka J, et al. Identification of Three Rheumatoid Arthritis Disease Subtypes by Machine Learning Integration of Synovial Histologic Features and RNA Sequencing Data. Arthritis Rheumatol (2018) 70:690–701. doi: 10.1002/art.40428

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Koopman FA, van Maanen MA, Vervoordeldonk MJ, Tak PP. Balancing the autonomic nervous system to reduce inflammation in rheumatoid arthritis. J Intern Med (2017) 282:64–75. doi: 10.1111/joim.12626

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Ferrero E, Li Yim A, Maratou K, Lewis HD, Royal G, Tough DF, et al. Novel Insights into Rheumatoid Arthritis Through Characterisation of Concordant Changes in DNA Methylation and Gene Expression in Synovial Biopsies of Patients with Differing Numbers of Swollen Joints. SSRN Electron J (2020). doi: 10.2139/ssrn.3576744

CrossRef Full Text | Google Scholar