Article Text
PDF
Abstract
Objectives Synovial fibroblasts (SFs) are one of the major components of the inflamed synovium in rheumatoid arthritis (RA). We aimed to gain insight into the pathogenic mechanisms of SFs through elucidating the genetic contribution to molecular regulatory networks under inflammatory condition.
Methods SFs from RA and osteoarthritis (OA) patients (n=30 each) were stimulated with eight different cytokines (interferon (IFN)-α, IFN-γ, tumour necrosis factor-α, interleukin (IL)-1β, IL-6/sIL-6R, IL-17, transforming growth factor-β1, IL-18) or a combination of all 8 (8-mix). Peripheral blood mononuclear cells were fractioned into five immune cell subsets (CD4+ T cells, CD8+ T cells, B cells, natural killer (NK) cells, monocytes). Integrative analyses including mRNA expression, histone modifications (H3K27ac, H3K4me1, H3K4me3), three-dimensional (3D) genome architecture and genetic variations of single nucleotide polymorphisms (SNPs) were performed.
Results Unstimulated RASFs differed markedly from OASFs in the transcriptome and epigenome. Meanwhile, most of the responses to stimulations were shared between the diseases. Activated SFs expressed pathogenic genes, including CD40 whose induction by IFN-γ was significantly affected by an RA risk SNP (rs6074022). On chromatin remodelling in activated SFs, RA risk loci were enriched in clusters of enhancers (super-enhancers; SEs) induced by synergistic proinflammatory cytokines. An RA risk SNP (rs28411362), located in an SE under synergistically acting cytokines, formed 3D contact with the promoter of metal-regulatory transcription factor-1 (MTF1) gene, whose binding motif showed significant enrichment in stimulation specific-SEs. Consistently, inhibition of MTF1 suppressed cytokine and chemokine production from SFs and ameliorated mice model of arthritis.
Conclusions Our findings established the dynamic landscape of activated SFs and yielded potential therapeutic targets associated with genetic risk of RA.
- synovitis
- fibroblasts
- arthritis
- rheumatoid
- cytokines
- autoimmune diseases
Statistics from Altmetric.com
Key messages
What is already known about this subject?
In rheumatoid arthritis (RA), a variety of dysregulated molecules from immune cells and mesenchymal cells drive disease progression. Synovial fibroblasts (SFs), the most abundant resident mesenchymal cells in the inflamed synovium, produce a variety of pathogenic molecules including interleukin 6.
Genome-wide association studies have identified more than 100 RA susceptibility loci. To gain insight into the pathogenic mechanisms of SFs, understanding the genetic contribution to molecular regulatory networks under inflammatory condition is crucial.
What does this study add?
Integrated analyses of activated SFs revealed an overview of relationships between pathogenic gene expressions, epigenomic modulations and RA susceptibility loci.
Chromatin remodelling induced by synergistic proinflammatory cytokines were associated with RA heritability, and some transcription factors (metal-regulatory transcription factor-1, runt-related transcription factor 1) could be crucial for the structural rearrangement and the formation of inflammatory arthritis.
How might this impact on clinical practice or future developments?
Our findings established the dynamic landscape of activated SFs and yielded potential therapeutic targets associated with genetic risk of RA.
Introduction
Rheumatoid arthritis (RA) causes persistent synovitis leading to disabling joint destruction. Current treatment strategies that target cytokines (eg, tumour necrosis factor-α (TNF-α), interleukin (IL)-6), cell surface proteins (eg, CD20, CD80/86) or signalling molecules (eg, Janus kinase) have brought a paradigm shift in RA treatment. However, achieving sustained remission is still challenging even with such agents.1 Although the concepts of targeting multiple molecules have been proposed, combination or bispecific antibodies (anti-TNF-α and anti-IL-1β or anti-IL-17) failed to improve therapeutic efficacy.2 These findings imply that some unknown factors play critical roles in the progression of synovitis.
In the pathogenesis of RA, the activities of a variety of dysregulated molecules in immune cells and mesenchymal cells are orchestrated by genetic and environmental factors.3 To date, more than 100 RA susceptibility loci have been identified in genome-wide association studies (GWAS).4 Recent genetic studies of autoimmune diseases have reported that the majority (>90%) of these risk variants are located in non-coding regions and regulate gene expression in a cell type-specific manner,5 partly in an environment-specific fashion.6 An integrated understanding of the risk variants’ contribution to gene regulatory networks is crucial to gain insight into the pathogenic mechanisms of RA.
Synovial fibroblasts (SFs), the most abundant resident mesenchymal cells in the synovium, are major local effectors in the initiation and perpetuation of destructive joint inflammation through their production of a variety of pathogenic molecules including IL-6.3 Previous multiomics data of unstimulated SFs have proposed activated pathways in RASFs.7 However, a comprehensive picture of SFs’ contribution to RA pathogenesis has largely remained elusive, perhaps due to their complex features that mutate in response to the proinflammatory milieu.8 To date, a number of single cytokines that induce the inflammatory behaviour of SFs have been reported (eg, interferon (IFN)-γ and IL-17 from T cells and TNF-α, IL-1β, IFN-α, IL-18 and transforming growth factor-β1 (TGF-β1) from monocytes).9 To make matters more complicated, in in vivo, SFs are expected to be exposed to a more complex environment. Data show that some cytokine combinations (eg, TNF-α and IL-17) synergistically enhance the expression of inflammatory molecules.10 Those findings emphasise the need to simultaneously analyse the mechanisms underlying the accelerated inflammatory behaviour of SFs in the presence of single cytokines and cytokine combinations.
Here, we used integrative methods to analyse genomic, transcriptomic and epigenomic features of RASFs in the presence of various proinflammatory cytokines in RA joints.11 Analyses of activated SFs revealed an overview of relationships between pathogenic gene expressions, epigenomic modulations and RA susceptibility loci. Chromatin remodelling induced by synergistic proinflammatory cytokines were associated with RA heritability, and some transcription factors (metal-regulatory transcription factor-1 (MTF1), runt-related transcription factor 1 (RUNX1)) could be crucial for the structural rearrangement and the formation of inflammatory arthritis.
Results
RASFs display distinctive transcriptomic and epigenomic signatures
We stimulated cultured SFs from RA and osteoarthritis (OA) patients (n=30 each) with eight different cytokines (IFN-α, IFN-γ, TNF-α, IL-1β, IL-6/sIL-6R, IL-17, TGF-β1, IL-18) or a combination of all 8 (8-mix). We also fractionated peripheral blood mononuclear cells (PBMCs) from the same patients into five major immune cell subsets (CD4+ T cells, CD8+ T cells, B cells, natural killer (NK) cells, monocytes) (figure 1).
Experimental design for integrative analysis of activated SFs from RA and OA patients. Our study design included SFs stimulated by eight different factors plus a combination of all the factors. Specifically, cells were treated for 24 hours with one of the following: IFN-α 100 U/mL, IFN-γ 200 U/mL, TNF-α 10 ng/mL, IL-1β 10 ng/mL, IL-6/sIL-6R 200 ng/mL, IL-17 10 ng/mL, TGF-β1 10 ng/mL or IL-18 100 ng/mL or 8-mix, a mixture of the above eight cytokines. In addition, we used five freshly isolated PBMC populations (CD4+ T cells, CD8+ T cells, B cells, NK cells, monocytes) from the same patient cohort. RNA sequencing of individual samples from RA and OA patients (n=30 per each) was carried out, and ChIP sequencing and Hi-C analysis were conducted with pooled samples. SNP genotyping array was performed in all patients. IFN, interferon; IL, interleukin; NK cells, natural killer cells; NS, non-stimulated; OA, osteoarthritis; PBMCs, peripheral blood mononuclear cells; RA, rheumatoid arthritis; SFs, synovial fibroblasts; SNP, single nucleotide polymorphism; TGF-β1, transforming growth factor-β1; TNF-α, tumour necrosis factor-α.
Under non-stimulatory condition, RASFs and OASFs showed profound differences at the level of transcriptome (figure 2A) and epigenome (online supplemental figure 1A). Genes highly expressed in RASFs were enriched in cell adhesion associated pathways (online supplemental figure 1B), which is consistent with previous reports,12 13 and the differences were largely preserved after stimulation (figure 2B, (online supplemental figure 1C)). Although most of the genes induced by each cytokine stimulation were shared between the diseases, some genes showed different behaviours (online supplemental figure 1D). For instance, SOCS5, coding a cytokine signalling regulator, was induced by IL-1β in OASFs, though it was highly expressed even under non-stimulatory condition in RASFs (figure 2C). Epigenomic data supported the activity of an enhancer located in SOCS5 intron region in non-stimulated RASFs (figure 2D), which might be due to long-time cytokine exposure. Moreover, CXCL11, coding a chemokine associated with T cell chemotaxis, was highly induced by IFN-γ in RASFs compared with OASFs (figure 2E), and the corresponding difference in epigenome modification were seen (figure 2F). Genes differentially expressed between RASFs and OASFs, as well as those induced by each cytokine stimulation, are listed in online supplemental table 1). In addition, the relationship between the inflammatory environment and the subpopulation of SFs, which has been attracting attention in recent years, is discussed in online supplementary note 1.14
Supplemental material
Distinctive transcriptomic and epigenomic signatures in RASFs. (A) A volcano plot comparing RASFs and OASFs under non-stimulated condition. Red points mark the genes with significantly increased or decreased expression in RASFs (FDR<0.01). (B) Principal component analysis (PCA) of gene expression levels for the top 1000 variable genes. Samples projected onto PC1/PC2 (left) or PC3/PC4 (right). Numbers in parentheses indicate contribution ratio (percentage of variation) of the first 4 PCs. Arrows link the centroid of indicated groups and adjusted to start from the origin. (C) Transcript abundances of SOCS5 from RNA sequencing data in SFs under non-stimulation and stimulation by IL-1β. boxes, IQR; whiskers, distribution; dots, outliers. (D) Organisation of transcriptional regulatory regions around SOCS5 gene in SFs under non-stimulation and stimulation by IL-1β. Boxed area indicate a putative enhancer. data were visualised using the Integrative Genomics Viewer (IGV). (E) Transcript abundances of CXCL11 from RNA sequencing data in SFs under non-stimulation and stimulation by IFN-γ. boxes, IQR; whiskers, distribution; dots, outliers. (F) Organisation of transcriptional regulatory regions around CXCL11 gene in SFs under non-stimulation and stimulation by IFN-γ. Boxed area indicates a putative enhancer. Data were visualised using the IGV. IFN, interferon; IL, interleukin; NS, non-stimulated; n.s., not significant; OA, osteoarthritis; PCA, principal component analysis; RA, rheumatoid arthritis; SFs, synovial fibroblasts; TGF-β1, transforming growth factor-β1; TNF-α, tumour necrosis factor-α.
Gene regulatory effects of RA risk loci are modulated by an environmental perturbation
We next performed cis-eQTL (expression quantitative trait locus) analysis to evaluate the effect of genetic variants on gene expressions. Together with tissue-by-tissue analysis, we also used Meta-Tissue software (online supplementary materials) for a meta-analysis across SFs under 10 stimulations and 5 PBMC subsets. To improve analytical ability, we jointly analysed the RA and OA samples, and afterwards we compared the effect sizes between the diseases for significant eQTLs. We considered variants with FDR <0.1 in single tissue analysis or m value >0.9 in meta-analysis as significant eQTLs.15 16
As a result, 3245–4118 genes in SFs and 2557–2828 genes in PBMCs had significant eQTLs (online supplemental figure 2A). In total, 2368 genes showed eQTL effects only in SFs (online supplemental figure 2B). Although similar eQTL effect sizes were observed in RASFs and OASFs, when the analyses were limited to differential peak regions between the diseases, some loci showed genome-wide significant difference (online supplemental figure 2C). In total, 12 loci showed significant interaction with the disease term (FDR <0.1, (online supplemental table 2). One example is rs35293523-LPAR1 in non-stimulated SFs (figure 3A, RA: p=2.1 × 10-6, OA: p=0.114, interaction FDR=0.0078). This locus is located in H3K27ac peak which is significantly greater in RASFs (figure 3B, p=2.93×10-8), indicating its enhancer activity in RA. As this single nucleotide polymorphism (SNP) is in tight LD with bone mineral density associated GWAS variant (rs12338959, r2=0.99 in East Asian (EAS), r2=0.79 in European (EUR) population),17 epigenomic modulation of this locus could be associated with osteoporosis, a common comorbidity of RA.18
Supplemental material
Disease-specific function of RA genetic risk loci. (A) Expression of LPAR1 in OASFs (left) and RASFs (right) under non-stimulated condition. Individuals are stratified according to the rs35293523 genotype. Nominal p values in eQTL mapping are shown. (B) Organisation of transcriptional regulatory regions around LPAR1 gene and positional relationship of RA risk loci (blue triangle, rs35293523) in OASFs and RASFs under non-stimulated condition. Data were visualised using the Integrative Genomics Viewer. eQTL, expression quantitative trait locus; OA, osteoarthritis; RA, rheumatoid arthritis; SFs, synovial fibroblasts.
On the other hand, some of the RA GWAS loci had eQTL effects in SFs in a stimulation dependent manner (online supplemental table 3), (online supplemental figure 3). One example is rs6074022, which is in tight LD (r2=0.95 in EUR, r2=0.9 in EAS population) with an established RA risk SNP rs4810485.19 rs6074022 had robust eQTL effects on CD40 in SFs, especially under IFN-γ or 8-mix stimulations (figure 4A,B). Importantly, the presence of an active regulatory region at rs6074022 was inferred only under these conditions (figure 4C). Although CD40 is also expressed by B cells (figure 4D), the eQTL effect was not observed at this locus (figure 4B).
Supplemental material
Stimulation-specific function of RA genetic risk loci. (A) A dot plot of rs6074022-CD40 Cis-eQTL meta-analysis posterior probability m values versus tissue-by-tissue analysis -log10 p value. The grey solid line (m-value=0.9) corresponds to the significance threshold in this study. (B) Expression of CD40 in RASFs (filled circles) and OASFs (open circles) stimulated by IFN-γ (left), 8-mix (middle) and B cells (right) from each individual plotted according to the rs6074022 genotype. Nominal p values in eQTL mapping are shown. (C) Transcriptional regulatory regions around the CD40 gene and positional relationship of rs6074022 (blue triangle) in stimulated SFs and PBMCs. IRF1 biding sites were obtained from the public epigenome browser ChIP-Atras. Data were visualised using the Integrative Genomics Viewer. (D) Transcript abundance of CD40 from RNA sequencing data in stimulated SFs and PBMCs. (E) A volcano plot of differential gene expression analysis comparing the presence or absence of CD40 ligand (CD40L) for IFN-γ-stimulated SFs. Orange and blue points mark the genes with significantly increased or decreased expression, respectively, for the addition of CD40L (FDR<0.01). Boxes, IQR; whiskers, distribution; dots, outliers in (B) and (D). CPM, counts per million; eQTL, expression quantitative trait locus; IFN, interferon; IL, interleukin; NK cells. natural killer cells; NS, non-stimulated; n.s., not significant; OA, osteoarthritis; PBMCs, peripheral blood mononuclear cells; RA, rheumatoid arthritis; SFs, synovial fibroblasts; TGF-β1, transforming growth factor-β1; TNF-α, tumour necrosis factor-α.
The biological role of the CD40-CD40L pathway in SFs has been discussed.20 21 We performed transcriptomic analysis of RASFs stimulated with a 2-trimer form of the CD40 ligand and IFN-γ. As a result, some cytokines (eg, IL6) and chemokines (eg, CCL5, CXCL10) were significantly upregulated by the ligation of CD40L (figure 4E). Taken together, we conjecture that CD40 expression in SFs is influenced by genetic and environmental predisposition, and the CD40-CD40L pathway might have a pathogenic role in RASFs.
Genetic risk of RA accumulate in the transcriptomic and epigenomic perturbations induced by synergistic proinflammatory cytokines
As acquired epigenomic differences between RA and OA are not necessarily associated with inherited genome, we next sought to reveal the condition most closely associated with RA susceptibility.
First, to elucidate the transcriptomic changes that is associated with RA genetic risk, we performed gene-set enrichment analysis using MAGMA software (online supplementary materials). This analysis indicated that perturbed gene sets subject to IFN-α, IFN-γ and 8-mix stimulation significantly overlapped with RA risk loci in both EUR and EAS populations (online supplemental figure 4AB). Those data contrasted with the non-significant association of transcriptome differences between the diseases with RA genetic risk. These findings indicate that there is an accumulation of RA genetic risk in the pathways that are perturbed under specific stimulatory conditions in cultured SFs.
Next, to elucidate the epigenomic changes which is associated with genetic risk, we assessed the enrichment of GWAS top-associated loci in regulatory regions including super-enhancers (SEs). SEs are large clusters of enhancers collectively bound by an array of transcription factors (TFs) to define cell identity, and they are hotspots for disease susceptibility.22–24 Although the significant overlap of SEs in Th cells and B cells with RA risk loci has been reported,22 SFs have not been examined. Here we compared the enrichment of RA risk loci to SEs and typical-enhancers (TEs), and risk loci showed significant enrichment with SEs in CD4+ T cells and B cells, as well as with SEs in 8-mix stimulated SFs (figure 5A). Some risk loci showed overlap with SEs which appear uniquely under 8-mix treatment (figure 5B,C). When we performed a similar analysis using risk loci for type 1 diabetes mellitus (a representative non-articular autoimmune disease), only SEs in CD4+ T cells and B cells showed significant enrichment (online supplemental figure 4C). The number or width of 8-mix SEs were comparable to those in other stimulations (online supplemental figure 5). Consequently, SFs might behave as key players in RA pathogenesis especially under combination of cytokines.
Enrichment of RA genetic risk in SFs SEs under eight cytokine stimulation. (A) Enrichment of RA risk loci in transcriptional regulatory regions of stimulated SFs and PBMCs. Active enhancers were classified into SEs and TEs following standard rose algorithms. The red solid lines and the black solid lines are the cutoffs for Bonferroni significance and p=0.05, respectively. (B) A Circus plot showing the overlap of SEs in SFs under different stimulatory conditions. only the regions unique to each condition or common to all of the conditions are depicted. (C) A Circus plot showing the overlap of RA risk loci and SEs in SFs under different stimulatory conditions. IFN, interferon; IL, interleukin; NK cells, natural killer cells; NS, non-stimulated; OA, osteoarthritis; PBMCs, peripheral blood mononuclear cells; RA, rheumatoid arthritis; SE, super-enhancer; SFs, synovial fibroblasts; TE, typical enhancer; TGF-β1, transforming growth factor-β1; TNF-α, tumour necrosis factor-α.
SEs induced by cytokine mixtures regulate genes crucial for RA pathogenesis
While the 8-mix is an artificial stimulatory condition, its unique impact on epigenomes (figure 5B), synergistic enhancement of arthritic gene expressions (online supplemental figure 6), and significant association with RA genetic risk (online supplemental figure 4A,B, figure 5A) induced us to interpret that it reflects some aspect of SF behaviour which is relevant to RA pathogenesis. The contribution of each cytokine to 8-mix condition is discussed in online supplementary note 2.
In order to characterise the genes regulated by 8-mix SEs, we combined the three-dimensional (3D) genome architectures (chromatin loops detected by Hi-C analysis), the position of SEs, promoter regions (defined with H3K4me3 ChIP sequencing analysis). We annotated ‘SE-contacted genes’ such that one side of Hi-C loop anchors overlapped an SE, the other side coincided with the transcriptional start site (TSS) and coexisted with the H3K4me3 peak (figure 6A). SEs were highly overlapped with Hi-C loop anchors than were TEs or H3K4me1 peaks (online supplemental figure 7A). When the TSS and H3K27ac peak was connected by a Hi-C loop, the variation of mRNA expression showed a significant correlation with the H3K27ac peak variation (online supplemental figure 7B). These results underscore the validity of connecting active enhancer marks and TSS by Hi-C loops as previous reports.25 26 When we compared the expression of SE-contacted genes and TE-contacted genes, the former showed significantly higher expression than the latter (online supplemental figure 7C).24
RA pathogenic factors regulated by SEs in SFs treated with eight cytokines. (A) A schematic image of ‘SE-contacted genes’. (B) A Venn diagram representing the overlap of TE-contacted (top) or SE-contacted (bottom) genes in SFs under different stimulatory conditions. Red, blue and black text highlight genes whose contacted SEs overlap with RA risk loci, cytokines and chemokines and transcription factors, respectively. NS, non-stimulated; RA, rheumatoid arthritis; SE, super-enhancer; SFs, synovial fibroblasts; TE, typical enhancer; TNF-α, tumour necrosis factor-α; TSS, transcriptional start site.
Next, we compared genes contacted by either SEs or TEs of SFs under three different conditions: non-stimulated, TNF-α or the 8-mix. The proportions of overlap between these conditions were smaller in SE-contacted genes (9.7%) than TE-contacted genes (14.5%) (figure 6B), indicating that the stimulation-specific expression profile and SEs formation are associated. SE-contacted genes included a number of TFs (eg, MTF1, RUNX1), cytokines (eg, IL6) and chemokines (eg, CCL5, CCL8) (figure 6B, (online supplemental table 4)).
Supplemental material
IL6 is a representative example of an 8-mix SE-contacted gene. Although this gene is regulated by an SE (almost 30 kb long) that exists upstream of the TSS in non-stimulated or TNF-α stimulated SFs, this SE elongates to 70 kb long under the 8-mix and an additional Hi-C loop emerges (figure 7A). When we inhibited SE formation with JQ1, a BRD4 inhibitor, the increased IL-6 expression under the 8-mix was disturbed (figure 7B). The necessity of SEs for elevated IL-6 production under synergistic inflammation was inferred.
Representative SE-contacted genes in SFs treated with eight cytokines. (A, C, E) Organisation of transcriptional regulatory regions around IL6 (A), RUNX1 (C) and MTF1 (E) genes and positional relationship of RA risk loci (blue triangle, rs8133848 for RUNX1 and rs28411352 for MTF1) and chromatin conformation in stimulated SFs. Data were visualised using the Integrative Genomics Viewer. (B) Expression of IL-6 in SFs (n=10) treated with JQ1. The mRNA and protein expression were quantified by qRT-PCR (left) and ELISA (right), respectively. Horizontal crossbars, mean; error bars, SD. P values were determined using one-way ANOVA followed by Tukey’s multiple comparison test (**p<0.01, ***p<0.001). (D, F) Transcript abundances of RUNX1 isoform (RUNX1b) (D) and MTF1 (F) from RNA sequencing data for stimulated SFs and PBMCs. Boxes, IQR; whiskers, distribution; dots, outliers. ANOVA, analysis of variance; CPM, counts per million; IFN, interferon; IL, interleukin; MTF1, metal-regulatory transcription factor-1; NK cells, natural killer cells; NS, non-stimulated; OA, osteoarthritis; PBMCs, peripheral blood mononuclear cells; RA, rheumatoid arthritis; RUNX1, runt-related transcription factor 1; SE, super-enhancer; SFs, synovial fibroblasts; TGF-β1, transforming growth factor-β1; TNF-α, tumour necrosis factor-α; TPM, transcripts per million.
Another example is RUNX1, a master-regulator involved in hematopoiesis.27–29 In our study, RA risk locus rs8133843 overlapped with an 8-mix SE that exists upstream of RUNX1 (figure 7C). A Hi-C loop was formed with the promoter immediately above the second exon of RUNX1 only in the 8-mix, and the RUNX1 expression was higher in the 8-mix compared with others (figure 7D).
MTF1, a zinc finger TFs, is another example. RA risk locus rs28411352 overlapped with an 8-mix unique SE that exists upstream of the MTF1 (figure 7E). The Hi-C loop was detected with the promoter only under the 8-mix. MTF1 expression was upregulated in the 8-mix (figure 7F).
TFs associated with stimulation-induced SE formation control arthritis progression
Finally, we searched for candidate modulators that were crucial for SE formation, especially in the 8-mix. In the previous study, key TFs for SE formation were reported to be controlled by SEs themselves, forming a self-regulatory network.24 From this perspective, we used motif analysis to focus on SE-contacted TFs that were also enriched in 8-mix SEs (figure 8A). Among SE-contacted genes, TFs such as SNAI1, TCF4 and MTF1 showed significant motif enrichment in 8-mix SEs. MTF1 was the only example that also showed the overlap of 8-mix SE and RA risk variant (figure 7E). Although the RUNX1 motif was not enriched in 8-mix SEs compared with non-stimulated SEs, its motif was significantly enriched in SEs compared with the background sequence, both in 8-mix and without stimulation (p=1.0×10-16 and 1.0×10-20, respectively). In vitro validation analysis showed that the expression of 8-mix SE-contacted genes was significantly suppressed by MTF1 and RUNX1 knockdown (p=2.0×10-3 and 4.3×10−4, respectively) (figure 8B, (online supplemental figure 8). The effect of MTF1 knockdown was more pronounced in 8-mix SE-contacted genes than TE-contacted genes, as anticipated from its motif enrichment in SEs. In in vitro assay, the increased 8-mix SE-contacted gene expressions (IL-6, CCL5) from stimulated RASFs was reduced by treatment with APTO-253, a MTF1 inhibitor (figure 8C,D).30 Furthermore, in a collagen induced arthritis (CIA) model, APTO-253 demonstrated significant preventive (online supplemental figure 9) and therapeutic activity (figure 8E,F) on arthritis formation. Collectively, these results indicated that certain TFs play critical roles in the formation of epigenomic structures induced by synergistic proinflammatory cytokines and support MTF1 inhibition as a promising therapy candidate for RA (online supplemental figure 10).
Transcription factors associated with stimulation-induced SE formation and arthritis progression. (A) Table depicts transcription factor binding motifs enriched at SEs in stimulated SFs. Following are summarised: attribution to SE-contacted genes, relative enrichment p values to TEs in each stimulatory condition or to SEs of non-stimulated condition. (B) Expression of SE- or TE-contacted genes in SFs stimulated by eight cytokines in cells depleted of specified transcription factors (TCF4, SNAI1, MTF1 and RUNX1) relative to control SFs. Boxes, IQR; whiskers, distribution. P values, paired t-test (*p<0.05, ***p<0.001). (C, D) Expression of IL-6 and CCL5 in SFs (n=13) treated with APTO-253. The mRNA and protein level were quantified by qRT-PCR (C) and ELISA (D), respectively. Horizontal crossbars, mean; error bars, SD. P values, one-way ANOVA followed by Tukey’s multiple comparison test (*p<0.05, **p<0.01, ***p<0.001). (E, F) Therapeutic effect of APTO-253 on CIA model. Following the onset of arthritis, mice were intravenously injected with either control or 15 mg/kg APTO-253 for twice per day for two consecutive days per week. Clinical (E) and pathological scores (F). Dots, mean; error bars, SD. P values, Mann-Whitney U test (*p<0.05, ***p<0.001). AVOVA, analysis of variance; CIA, collagen-induced arthritis; IL, interleukin; MTF1, metal-regulatory transcription factor-1; NS, non-stimulated; RUNX1, runt-related transcription factor 1; SE, super-enhancer; SFs, synovial fibroblasts; TE, typical enhance; TNF-α, tumour necrosis factor-α.
Discussion
In this study, we were able to describe the dynamic landscape of activated SFs and their contribution to RA pathogenesis. RASFs demonstrated distinct epigenomic and transcriptomic features from OASFs, which might be due to acquired epigenomic modifications from long-standing inflammation. Importantly, heritable RA risk was enriched in stimulation dependent epigenomic structures in SFs, a topic that has not been addressed deeply in previous reports.
Recent large-scale eQTL studies enhanced our understanding of complex diseases.31 32 Considering the importance of studying disease-relevant tissues for functional understanding of GWAS variants,15 33 we conducted cis-eQTL analysis of SFs, which are major local effector cells in arthritic joints. One example of eQTL-eGene pairs in SFs was the association of an RA risk SNP (rs4810485) with CD40 expression. Although no significant eQTL effects of this locus in B cells was observed in our data, a previous report showed the protein-level QTL association of CD40 in CD19+ B cells and this locus.34 On the contrary, a larger eQTL study by the GTEx consortium that showed genome-wide significant eQTL effects of rs6074022 on CD40 in lung (mainly fibroblasts) (p=3.1×10-26) and cultured fibroblasts (p=8.4×10-16), but not in EBV-transformed lymphocytes (p=0.04).35 Naturally, although the possible importance of CD40 signals in B cells for RA pathogenesis cannot be neglected, our results shed light on the role of CD40-CD40L signals in a genetic network of synovitis.
The GWAS-SEs enrichment analysis indicated the importance of SFs under synergistic stimulations in the development of RA (analogous to CD4+ T cells and B cells). During cytokine synergy, Hi-C analysis suggested that there were dynamic conformational changes in 3D structures involving SEs and the promoter of pathological molecules. We found marked expression of 8-mix stimulated SE-contacted genes (ie, IL6, RUNX1, MTF1).
SEs can collapse when their cofactors (eg, BET family) are perturbed.36 On the other hand, selective modulation of disease-associated SEs in a cell type-specific manner may have better safety profiles than pan-SEs inhibitors (ie, JQ1). In the present study, we analysed TFs that have the potential to be selective SE modulators in activated SFs. Our results suggested that MTF1 participates in SE formation, putatively making a feedback loop to maintain the epigenomic machinery. MTF1 regulates gene expression in response to zinc and various stresses.37 In the setting of disease, MTF1 could contribute to tumour metastasis and chemoresistance.38 In the present study, APTO-253, a MTF1 inhibitor,30 demonstrated inhibitory activity on the expression of pathogenic molecules (IL-6, CCL5) from RASFs and antiarthritic activity in a CIA model. Previously, the zinc-ZIP8-MTF1 axis was identified as a catabolic regulator of cartilage destruction.39 Furthermore, intra-articular injection of adenovirus expressing-MTF1 in an OA mouse model promoted the expression of various inflammatory molecules in SFs. This evidence supports the essential role of MTF1 as a genetic hub for joint destruction and inflammation.
There are some limitations of this study. First, the number of patients included in the cis-eQTL analysis was limited owing to sample accessibility, resulting in putatively large false negatives. Second, previous reports have shown that culture procedure could affect the phenotype of SFs.39–42 It should be noted that the results presented were obtained from the early passage SFs, and using directly isolated cells may have led to less biased analysis. Third, the stimulation endpoint was only 24 hours, and more detailed effects of each cytokine could have been analysed if expression information at different stimulation endpoints had been available. Finally, isolated PBMCs were not artificially stimulated with cytokines ex vivo. Thus, it is not clear whether the eQTL difference between SFs and PBMCs is attributable solely to cell type difference.
Overall, our multiomics approach using activated local cells in joints shed light on the concept of SF-targeted therapy from the perspective of epigenome remodelling related to genetic risk and would be beneficial in searching for novel drug targets.
Acknowledgments
We would like to thank G Inoue, M Abe, K Myouzen and K Kobayashi for their technical assistance at the Laboratory for Autoimmune Diseases, and Y Momozawa and M Kubo for providing technical advice at the Laboratory for Genotyping Development, RIKEN. We received generous support from all the physicians who participated in sample collection at the Department of Orthopaedic Surgery, the University of Tokyo.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
Handling editor Josef S Smolen
HT and MO contributed equally.
Contributors HT, SS, KI, YKo, KY and KF designed the research project. MO conducted bioinformatics analysis on the advice of KI, AS performed RNA sequencing. HT performed ChIP sequencing and other in vitro experiments. TS and KS performed Hi-C. HT, YT, HI, JH, YKa and ST contributed human samples. SS performed figure editing. HT and MO wrote the manuscript with critical inputs from YKo, KY and KF.
Funding This research was supported by funding from Takeda Pharmaceutical (YKo, KY and KF), the Ministry of Health, Labour and Welfare, Ministry of Education, Culture, Sports, Science and Technology KAKENHI Grant-in-Aid for Scientific Research (B) (18H02846) and Grant-in-Aid for Scientific Research (C) (17K09972) from the Japan Society for the Promotion of Science.
Competing interests None declared.
Patient consent for publication Not required.
Provenance and peer review Not commissioned; externally peer reviewed.
Data availability statement Data are available in a public, open access repository. The datasets generated during this study are available at the National Bioscience Database Center (NBDC) with the study accession code hum0207.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.