Post-GWAS functional studies reveal an RA-associated CD40-induced NF-kB signal transduction and transcriptional regulation network targeted by class II HDAC inhibitors.

Currently, it remains difficult to identify which single nucleotide polymorphisms (SNPs) identified by genome-wide association studies (GWAS) are functional and how various functional SNPs (fSNPs) interact and contribute to disease susceptibility. GWAS have identified a CD40 locus that is associated with rheumatoid arthritis (RA). We previously used two techniques developed in our laboratory, single nucleotide polymorphism-next-generation sequencing (SNP-seq) and flanking restriction enhanced DNA pulldown-mass spectrometry (FREP-MS), to determine that the RA risk gene RBPJ regulates CD40 expression via a fSNP at the RA-associated CD40 locus. In the present work, by applying the same approach, we report the identification of six proteins that regulate RBPJ expression via binding to two fSNPs on the RA-associated RBPJ locus. Using these findings, together with the published data, we constructed an RA-associated signal transduction and transcriptional regulation network (STTRN) that functionally connects multiple RA-associated risk genes via transcriptional regulation networks (TRNs) linked by CD40-induced nuclear factor kappa B (NF-kB) signaling. Remarkably, this STTRN provides insight into the potential mechanism of action for the histone deacetylase inhibitor givinostat, an approved therapy for systemic juvenile idiopathic arthritis. Thus, the generation of disease-associated STTRNs based on post-GWAS functional studies is demonstrated as a novel and effective approach to apply GWAS for mechanistic studies and target identification.

Introduction

Rheumatoid arthritis (RA) is an autoimmune disease that affects 0.5–1% of the population. It is a long-term, progressive and disabling disease, resulting in destructive inflammation, swelling and pain in the joints and other tissues and organs. The pathogenesis of RA remains incompletely understood and there is no effective cure (1).

As a complex disease, RA susceptibility is believed to be associated with multiple genes in combination with environmental factors. While genome-wide association studies (GWAS) have identified and validated more than 100 genetic loci that are associated with RA (2), post-GWAS functional studies have proven to be extremely difficult owing to the fact that the single nucleotide polymorphisms (SNPs) used by GWAS to pinpoint these genetic loci are simply a proxy for large stretches of DNA (haplotypes) that contain many other SNPs in linkage disequilibrium (LD), most of which reside in non-coding DNA (3). Similarly, even though a handful of environmental factors have been documented to have relatively consistent associations with RA (4), delineating how each of the environmental factors contributes to RA pathogenesis is equally difficult without understanding the biological mechanisms underlying the contribution of RA risk genes to the disease.

CD40 is a member of the TNF-receptor superfamily, which is mainly expressed on antigen presenting cells such as B cells, dendritic cells and macrophages. Upon ligation of CD40 with its ligand CD40L, a multitude of signaling pathways are activated, including the nuclear factor kappa B (NF-kB) pathway, which mediate a broad range of immune and inflammatory responses (5). Recent large GWAS have validated seven RA associated loci that are involved in CD40-induced NF-kB signaling; they are: CD40, RBPJ, TNFAIP3, TRAF1, TRAF6, NFKBIE and REL (6). Previously, we demonstrated that the disease-associated CD40 locus is a gain-of-function allele with more CD40 expression on the risk allele, which, upon CD40L activation, results in a high level of activation of NF-kB p65 (RelA) by phosphorylation (7). In order to understand the mechanisms underlying the contribution of RA-associated functional SNPs (fSNPs), recently, we developed two novel techniques: single nucleotide polymorphism-next generation sequencing (SNP-seq) and flanking restriction enhanced DNA pulldown-mass spectrometry (FREP-MS) (8). By coupling SNP-seq with FREP-MS, we identified RBPJ as one of the four transcriptional regulators modulating CD40 expression via three fSNPs on the CD40 locus (8). Interestingly, RBPJ itself is also an RA risk gene (6,9).

The CD40–CD40L axis is involved in many diseases, including autoimmune diseases and cardiovascular diseases (10,11). Blocking the interaction between CD40 and CD40L has proved effective in preventing these diseases or in ameliorating their symptoms. This suggests that CD40 is a potential target for drug development (11). However, there is currently no approved drug that inhibits CD40–CD40L signaling since early efforts induced significant side effects (7). Human genetics suggests that inhibiting CD40-mediated intracellular signaling could be effective in humans (7). However, human genetics also suggests that CD40 is not the only risk gene linked to RA. Indeed, it may very well be a combination of multiple risk genes, which seemingly have no connection, that form a presently uncharacterized genetic network that ultimately determines a person’s overall immune response and susceptibility to RA.

Systems biology is a study of biological systems that integrate key elements such as DNA, RNA, proteins and cells as a whole network (12). Based on the data generated by GWAS, we hypothesize that there are disease-associated risk gene transcriptional regulation networks (TRNs) linked by signal transduction pathways in a cell type-specific fashion. We refer to this network as signal transduction and transcriptional regulation network (STTRN). Obviously, the pathogenesis and/or susceptibility of a complex disease is not because of the signal transduction pathway itself, but the perturbation of a set of risk gene expression, presumably by transcriptional regulation via the risk alleles of many causative fSNPs. In such a network, change in one single gene expression in response to an environmental factor could alter the expression of all the risk genes in this network. If all the alleles in this risk network predispose toward disease risk, a destructive signal, possibly initiated by an environmental factor such as smoking, would be amplified. However, at the current stage, our knowledge on gene transcriptional regulation, especially on disease-associated risk gene transcriptional regulation, is extremely limited; therefore, it is difficult to construct a disease-associated STTRN. GWAS, as a tool to study complex disease, provides a window of opportunity to collect the disease-associated functional information required to build such a STTRN.

In this study, using SNP-seq and FREP-MS, together with allele-imbalanced DNA pulldown-western blot (AIDP-Wb), a novel DNA pulldown assay recently developed in our laboratory to specifically detect the allele-imbalanced binding of a known protein to a fSNP (13), we screened a library containing 156 SNPs collected from the seven aforementioned RA loci in LD with r2 > 0.8 on the CD40-induced NF-kB signal pathway. We validated two fSNPs on the RBPJ locus together with six regulatory proteins that modulate the RBPJ expression. Based on these findings, as well as previous publications, we constructed and verified an RA-associated STTRN involving the CD40-induced NF-kB signaling, using both RNAi gene knockdown and histone deacetylase (HDAC) inhibition in human B cells and fibroblast-like synoviocytes (FLS).

Results

Identification and validation of two fSNPs on the RA-associated RBPJ locus

Previously, we introduced SNP-seq to identify four RA-associated fSNPs on the RA-associated CD40 locus (8). In the present work, we applied SNP-seq to screen a DNA library containing 156 SNPs in LD with r2 > 0.8 collected from the seven RA risk loci involved in the CD40-induced NF-kB pathway, including CD40, RBPJ, TRAF1, TRAF6, NFKBIE, REL and TNFAIP3 (6), as outlined in Supplementary Material, Figure S1A. The quality of the screen is evidenced by high reproducibility between two technical repeats with a correlation coefficient demonstrating an r2 = 0.8893 (Supplementary Material, Fig. S1B). As a result, we identified 69 candidate fSNPs with a P-value < 0.05 and a Slope > 0.05 and <−0.05 (Supplementary Material, Table S1). Within these 69 candidate fSNPs, we recovered the four fSNPs (rs4810485, rs1883832, rs6032664 and rs6056926) on the RA-associated CD40 locus which we previously identified and confirmed (8). Among the remaining 65 candidate fSNPs, we identified six SNPs (rs10517086, rs17630466, rs4473645, rs35944082, rs874040 and rs6448434) among 13 SNPs on an RA-associated RBPJ locus in LD with rs874040 (r2 > 0.8). Among these six candidate fSNPs, five SNPs except for rs6448434 were demonstrated as fSNPs by electrophoretic mobility shift assay (EMSA) as each identified fSNP showed allele-imbalanced gel shifting (Fig. 1A). To further validate these five fSNPs, as a proof of evidence, we performed a luciferase reporter assay on two of these five fSNPs; they were rs17630466 and rs35944082. As noted, an allele-imbalanced luciferase activity between the two alleles of each of these two fSNP was evident (Fig. 1B). In both cases, the risk allele G has lower luciferase activity than the non-risk allele A, suggesting that both SNPs are likely bind to suppressive regulators.

Figure 1

Identification and characterization of rs17630446 and rs35944082 in the RA-associated RBPJ locus. (A) EMSA showing allele-imbalanced gel shifting on six candidate fSNPs on the RBPJ locus in LD with rs874040. Red nucleotides indicate the risk alleles on each SNP. Red arrows indicate the allele-imbalanced shifting bands. The data represent three biological replicates (n = 3). (B) Luciferase reporter assays demonstrating fSNPs rs17630466 and rs35944082. The data represent a combination of six biological replicates (n = 6). RLU: relative luminescent unit. (C) AIDP showing allele-imbalanced binding of LEF1 and MSH6 specifically bound to rs17630466 and IRF2, APOBEC3C, MEF2B and MEF2D to rs35944082. PARP-1 was used as internal loading control. The data represent three biological replicates (n = 3). (D) ChIP assays demonstrate the specific binding of LEF1 to rs17630466 and MEF2D to rs35944082. The data represent a combination of three biological replicates (n = 3). Ct: cycle threshold; IgG: immunoglobulin G (E) Motif analyses reveal the sequence similarity between the LEF1 binding sequence and the sequence from rs17630466 and between the IRF2 and MEF binding sequence and the sequence from rs35944082. Red nucleotides indicate SNPs. sh, shRNA; NS, non-significant.

We also performed in silico analysis on the 13 SNPs in LD with rs874040 (r2 > 0.8) using HaploReg 4.1, a web-based tool for epigenetic and functional annotation of genetic variants (14). Using Haploreg, we scored each of the 13 SNPs on a scale of 0–5 to reflect the number of positive annotations for promoter and enhancer histone methylation generated by ChIP-seq, DNase hypersensitivity by DNase-seq and predicted protein binding and alteration in binding motifs (Supplementary Material, Table S2). Among the six candidate SNPs we identified, both rs35944082 and rs874040 scored 5 (together with one non-fSNP), rs4473645 scored 3, rs17630466 and rs6448434 scored 2 and rs10517086 scored 1, suggesting a moderate enrichment of fSNPs in this analysis.

Identification of regulatory proteins that specifically bind to fSNPs rs17630466 and rs35944082 modulating RBPJ expression

All 156 SNPs in the library we screened are in the non-coding region of human DNA; therefore, they do not change protein coding sequences. Instead, they presumably exert their functions by binding to unknown regulatory proteins to control risk gene expression (15–20). To understand how the two fSNPs, rs17630466 and rs35944082, identified on the RBPJ locus modulate RBPJ expression, we performed FREP-MS on these two fSNPs using the G allele (risk allele) from both rs17630466 and rs35944082. We chose the G alleles since they show more protein binding in the allele-imbalanced gel shifting for both fSNPs (Fig. 1A). To increase specificity, we performed FREP-MS with these two fSNPs in parallel so that we can use these two fSNPs as specificity controls for each other (8). To further increase the fidelity, we performed this FREP-MS analysis with each SNP in duplicate. By doing so, based on the specific peptide spectrum counts shown in Supplementary Material, Table S3, we identified two proteins, LEF1 and MSH6, that specifically bind to rs17630466 and four proteins, IRF2, APOBEC3C, MEF2B and MEF2D, that bind to rs35944082. While all these proteins were reported as transcriptional regulators, MSH6 is a DNA mismatch repair protein that has never been documented as a transcription factor (21). Moreover, among these proteins, only MEF2 could be identified as a protein specifically binding to fSNP rs35944082 using in silico analysis (HaploReg 4.1; Supplementary Material, Table S2).

To demonstrate that LEF1 and MSH6 bind to rs17630466 and IRF2, APOBEC3C, MEF2B and MEF2D bind to rs35944082, respectively, we applied AIDP-Wb (13). Using this method, we observed an allele-imbalanced binding of LEF1 and MSH6 to rs17630466 and IRF2, APOBEC3C, MEF2B and MEF2D to rs35944082, respectively (Fig. 1C). In both cases, the risk allele G binds more of these factors than the non-risk allele A. The noted allele-imbalanced binding of these regulatory proteins to rs17630466 and rs35944082 showed in this assay further validates that rs17630466 and rs35944082 are indeed bona fide fSNPs. To further validate the specific binding of these proteins to the two fSNPs, as a proof of evidence, we performed chromatin immunoprecipitation (ChIP) assay on rs17630466 using an anti-LEF1 antibody and on rs35944082 using an anti-MEF2D antibody. A significant decrease in the binding of LEF1 to a DNA fragment containing rs17630466 and MEF2D to a DNA fragment containing rs35944082 was observed when we compared BL2 cells, a human B cell line, with either LEF1 or MEF2D shRNA knockdown to a scrambled shRNA control (Fig. 1D). In addition, the binding of LEF1 to rs17630466, as well as IRF2 and MEF2B/D to rs35944082, was further supported by the sequence similarity between these protein binding motifs (22,23) (ISMARA) and the sequences from the two fSNPs (Fig. 1E). Together, these data support that there is the specific binding of LEF1 and MSH6 to rs17630466 and the specific binding of IRF2, APOBEC3C, MEF2B and MEF2D to rs35944082.

To demonstrate that LEF1 and MSH6 as well as IRF2, APOBEC3C, MEF2B and MEF2D are the transcriptional regulators modulating RBPJ expression, we knocked down each of these proteins in BL2 cells using lentiviral shRNA. As a result, we observed that knockdown of each of these six proteins increases the expression of RBPJ at both the protein level, detected by western blots (Fig. 2A), and the mRNA level, detected by qPCR (Fig. 2B). These data, consistent with the data from the allele-imbalanced luciferase reporter assay presented in Figure 1B, suggest that all these six proteins function as RBPJ suppressors, with the risk allele G binding more of these proteins and having lower luciferase activity than the non-risk allele A. Previously, we showed that RBPJ is a transcription suppressor of CD40 via binding to the fSNP rs4810485 on the RA-associated CD40 locus (8). Consistent with this data, we also observed a decreased CD40 expression in all these six knockdown BL2 cell lines (Fig. 2), which further demonstrates the suppressive function of these six proteins in regulating RBPJ expression.

Figure 2

Western blots and qPCR assays demonstrate MSH6 and LEF1 as well as MEF2B, MEF2D, APOBEC3C and IRF2 as transcription suppressors of RBPJ in human BL2 cells. (A) Western blots and qPCR assays showing that knockdown of MSH6 and LEF1 with shRNA results in upregulation of RBPJ and downregulation of CD40 in BL2 cells. (B) Western blots and qPCR assays showing that knockdown of MEF2B, MEF2D, APOBEC3C and IRF2 with shRNA leads to upregulation of RBPJ and downregulation of CD40 in BL2 cells. ***P-value < 0.001. The western blot data represent three biological replicates (n = 3). The qPCR data represent a combination of three biological replicates (n = 3). sh, shRNA Ct: cycle threshold.

Class II HDAC inhibitors trigger a TRN from HDAC7 to MEF2D, RBPJ and CD40 in human B cells

MEF2D was recently reported to be associated with a systemic lupus erythematosus (SLE) sub-phenotype in Swedish cohorts (24). Like RA, SLE is also an autoimmune disease, resulting from the immune system mistakenly attacking its own body; therefore, these conditions share several clinical and pathological features as well as underlying biological mechanisms (25). Also, previously, it was reported that the expression of MEF2D is regulated by HDACs since treating human myocytes with the class II HDAC inhibitor MC1568 decreases MEF2D expression (26). MC1568, a derivative of (aryloxopropenyl)pyrrolyl hydroxyamide, is a novel, potent and specific inhibitor of class II HDACs. To determine if MC1568 regulates MEF2D expression in human B cells, we treated BL2 cells with different concentrations of MC1568. Significant downregulation of MEF2D was observed at both the mRNA and protein levels (Fig. 3A and C). To confirm this data, we used another HDAC inhibitor, ITF2357, also known as givinostat, to treat the BL2 cells. ITF2357 is a potent inhibitor of both class I and class II HDACs as well as an approved drug in the European Union for the treatment of systemic juvenile idiopathic arthritis (SJIA) (27). Our results showed that the treatment with ITF2357 also decreased the expression of MEF2D in BL2 cells (Fig. 3B and C). These data demonstrate that MEF2D, as a SLE risk gene, is regulated by class II HDACs. Again, consistent with our findings in the MEF2D shRNA knockdown BL2 cells (Fig. 2), we also observed a significant upregulation of RBPJ and downregulation of CD40 expression in either MC1568- or ITF2357-treated BL2 cells (Fig. 3A–C).

Figure 3

Effect of class II HDAC inhibitors on MEF2D as well as RBPJ and CD40 expression in BL2 cells. (A, B) Western blots showing downregulation of MEF2D, upregulation of RBPJ and decreased expression of CD40 in BL2 cells treated with HDAC inhibitors MC1568 and ITF2357. The data represent three biological replicates (n = 3). (C) qPCR showing downregulation of MEF2D, upregulation of RBPJ and decreased expression of CD40 in BL2 cells treated with HDAC inhibitors MC1568 and ITF2357. The data represent a combination of three biological replicates (n = 3). (D, E). qPCR showing that knockdown of HDAC7 and HDAC4 results in downregulation of MEF2D and upregulation of RBPJ, which, in turn, decrease CD40 expression in human primary FLS. The data represent a combination of three biological replicates (n = 3). ***P-value < 0.001. si, siRNA.

Both HDAC4 and HDAC7 belong to class II HDACs and can be inactivated by MC1568 (26,28). To identify which one of these two HDACs is involved in the regulation of MEF2D, we performed siRNA knockdown on HDAC7 and HDAC4 in human FLS. Downregulation of each of these genes resulted in decreased MEF2D expression, increased RBPJ expression and subsequent upregulation of CD40 expression (Fig. 3D and E), suggesting that both HDAC7 and HDAC4 can regulate the MEF2D expression, possibly in the same transcriptional complex. Coincidently, HDAC7 is also an RA risk gene (29).

Together, our findings reveal an RA-associated TRN in human cells linking HDAC4/HDAC7 to MEF2D, RBPJ and CD40. Moreover, this network is targeted by class II HDAC inhibitors.

An RA-associated STTRN involving CD40-induced NF-kB signaling in human B cells

Previously, we demonstrated that, with increased cell surface CD40 expression (as is the case for carriers of the RA risk allele), there is an increased activation of the classical NF-kB pathway as measured by the phosphorylation of NF-kB p65 in human B cells (7). To determine if the downregulation of CD40 induced by MEF2D inactivation affects NF-kB p65 activation, we first treated both the scrambled shRNA control and the MEF2D shRNA knockdown BL2 cells with MEGACD40L, a trimeric CD40 ligand, at 64 ng/ml for 15 min as we described previously (7). Activation of NF-kB p65 was evaluated with an anti-phosphorylated p65 antibody. Consistent with our previous results, we observed a reduction of phosphorylated NF-kB p65 in the MEF2D knockdown BL2 cells (Fig. 4A). We also checked the activation of p65 by MEGACD40L in the BL2 cells treated with either MC1568 or ITF2357. In both cases, we also observed reduced phospo-NF-kB p65 in treated cells (Fig. 4B and C).

Figure 4

Inactivation of NF-kB p65 in the MEF2D knockdown- and class II HDAC inhibitor-treated BL2 cells. (A) Western blots showing reduced phosphorylation of NF-kB p65 following MEF2D knockdown in BL2 cells activated by MEGACD40L. The data represent three biological replicates (n = 3). (B, C) Western blots showing reduced phosphorylation of NF-kB p65 in BL2 cells treated with MC1568 or ITF2357. Cells were activated with MEGACD40L following shRNA knockdown and MC1568 or ITF2357 treatment. The data represent three biological replicates (n = 3). p-p65: phosphorylated p65 (D) qPCR showing expression of the eight RA risk genes as well as IL6 and IL1β in the MEF2D knockdown BL2 cells activated by MEGACD40L. The data represent a combination of three biological replicates (n = 3). (E) qPCR showing expression of the eight RA risk genes as well as IL6 and IL1β in the BL2 cells treated with MC1568 or ITF2357. BL2 cells were activated by MEGACD40L after treatment. The data represent a combination of three biological replicates (n = 3). **P-value < 0.01 and ***P-value < 0.001. sh, shRNA.

NF-kB has long been known to play an important role in autoimmune diseases such as RA (30,31). CD40-induced NF-kB signaling in disease pathogenesis is highlighted by GWAS analysis in RA (2,32). Besides its role in the regulation of various pro-inflammatory genes, NF-kB also plays a critical role in regulating the expression of many other genes involved in cell survival, activation and differentiation of innate immune cells and inflammatory T cells (33). Based on TRRUST V2, a manually curated database of human and mouse transcriptional regulatory networks (34), we identified 17 NF-kB p65 targeted genes that are associated with high risk RA (GWAS catalog, 2109). Among these, eight were reported to be regulated by CD40 in human B cells: CD83 (35), CDK4/CDK6 (36), c-FLIP (37), ICAM1 (38), TNFAIP3 (39), IRF4 (40) and TRAF1 (41). We therefore performed qPCR on these eight RA risk genes, together with IL6 and IL1β, the two key proinflammatory cytokines, in the MEF2D knockdown BL2 cells treated with MEGACD40L. As a result, we observed a significant downregulation in the expression of all these eight RA risk genes plus IL6 and IL1β in the MEF2D knockdown BL2 cells compared with the scrambled controls (Fig. 4D). We also checked the expression of these eight RA risk genes plus IL6 and IL1β in BL2 cells activated by MEGACD40L after either ITF2357 or MC1568 treatment. While we observed a similar downregulation of all these eight genes plus IL6 and IL1β in the ITF2357-treated cells, a different expression pattern of these genes was identified in the MC1568-treated cells (Fig. 4E). The expression of CDK4 and CDK6 was not significantly changed, and the expression of ICAM1, IL6 and IL1β, instead of being downregulated in the ITF2357-treated cells, was significantly upregulated in the MC1568-treated cells. Of note, only ITF2357 is approved for clinical use.

Together, these data demonstrate that there is an RA-associated CD40-induced STTRN. In this STTRN, an RA risk gene TRN from HDAC7 to MEF2D, RBPJ and CD40 is linked via the CD40-induced NF-kB signal transduction to another RA risk gene TRN from NF-kB p65 to the eight RA risk genes, including IL6 and IL1β.

An RA-associated CD40-induced NF-kB STTRN in primary human FLS

FLS are the main stromal cells of the joint synovium; therefore, they play a critical role in the etiology (e.g. cartilage destruction) of RA (42). To determine whether the RA-associated CD40-induced NF-kB STTRN exists in FLS, primary human FLS isolated from a patient with osteroarthritis (OA) were used owing to the unavailability of normal FLS. To make use of this FLS, we generated all our results by comparing the treated FLS with the same untreated FLS control. Primary human FLS were treated with 1 μm ITF2357 for 24 h and then activated by 64 ng/ml MEGACD40L for 15 min. Proteins were isolated from both treated samples and untreated controls for western blot analysis to assess the expression of MEF2D, RBPJ, CD40 and phosphorylated NF-kB p65. Consistent with our results found in the ITF2357-treated human BL2 cells, we observed a similar downregulated MEF2D, increased RBPJ, decreased CD40 and subsequent decreased phosphorylation of NF-kB p65 in the ITF2357-treated FLS (Fig. 5A). To further explore the link between ITF2357 treatment and NF-kB signaling, we repeated the qPCR for the eight RA risk genes as well as IL6 and IL1β regulated by NF-kB p65 in FLS. A similar downregulation of these genes was observed in the ITF2357-treated FLS as was observed in the ITF2357-treated BL2 cells (Fig. 5B). Together, these data demonstrate the existence of a conserved RA-associated CD40-induced NF-kB STTRN in both human B cells and FLS.

Figure 5

Demonstration of an RA-associated CD40-induced NF-kB STTRN in FLS. (A) Western blots showing that ITF2357 (1 μm) treatment in FLS cells results in inactivation of NF-kB p65 through downregulation of CD40, which is caused by downregulation of MEF2D and upregulation of RBPJ. The data represent three biological replicates (n = 3). p-p65: phosphorylated p65 (B) qPCR showing downregulation of the risk genes CDK4, CDK6, c-FLIP, ICAM1, CD83, TNFAIP3, IRF4 and TRAF1 as well as IL6 and IL1β in the ITF2357-treated FLS. FLS were activated by MEGACD40L after ITF2357 treatment. The data represent a combination of three biological replicates (n = 3). (C, D) Scratch-wound assay showing decreased cell migration in the ITF2357-treated and the HDAC4/HDAC7 knockdown FLS. The data represent three biological replicates (n = 3). (E) Cell cycle analysis showing a G1 arrest in the ITF2357-treated and the HDAC4/HDAC7 knockdown FLS. The data represent three biological replicates (n = 3). (F) Caspase3/7 activity analysis indicating increased apoptosis in the ITF2357-treated and the HDAC4/HDAC7 knockdown FLS. The data represent three biological replicates (n = 3). ***P-value < 0.001. si: siRNA.

ITF2357, as a prescribed drug, decreases cell migration, blocks cell cycle progression and increases apoptosis via the RA-associated CD40-induced NF-kB STTRN

In RA, synovial lining increases through enhanced cell migration and proliferation, and/or decreased apoptosis of fibroblast-like or macrophage-like synoviocytes, which is responsible for the synovial hyperplasia that contributes to the destruction of cartilage and bone (43). Therefore, reversal of these cellular activities has been proposed as potential options for the therapeutic intervention for RA. Among the eight risk genes that are regulated by ITF2357 in the RA-associated CD40-induced NF-kB STTRN, ICAM1 is a cell adhesion molecule involved in cell migration (44). CDK4 and CDK6 play key roles in cell proliferation, and inhibition of CDK4/6 expression in different RA mice models ameliorates arthritis probably by suppressing synovial hyperplasia and joint destruction (45). c-FLIP is a master anti-apoptotic regulator (46). To assess whether the ITF2357-treated or the HDAC4/HDAC7 knockdown FLS reduces FLS migration, we performed a scratch-wound assay to measure basic cell migration (47). In both cases, impaired cell migration showed by delayed wound gap filling was observed (Fig. 5C and D), presumably caused by the decreased expression of ICAM1 upon either ITF2357 treatment or the HDAC4/HDAC7 siRNA knockdown. To determine if the ITF2357-treated or the HDAC4/HDAC7 siRNA knockdown FLS alter cell proliferation, we performed cell cycle analysis and observed a significant G0/G1 phase arrest and a corresponding decrease of S and G2/M (Fig. 5E), data that is consistent with the downregulation of CDK4 and CDK6 in these cells after treatment with ITF2357 or the knockdown of HDAC4/HDAC7. We also assessed apoptosis induced by an anti-Fas IgM monoclonal antibody (mAb) in the ITF2357-treated or the HDAC4/HDAC7 knockdown FLS. Consistent with a previous report (48), in both ITF2357-treated or the siRNA knockdown FLS, we detected a significant increase of apoptosis that is showed by an increased caspase 3/7 activity (Fig. 5F). Together, while these data reveal the possible therapeutic effects of ITF2357 on RA, they also demonstrate the existence of the RA-associated CD40-induced NF-kB STTRN.

Discussion

Previously, we demonstrated that RBPJ suppresses CD40 expression via binding to an RA-associated fSNP rs4810485 on the CD40 locus in human B cells and FLS (8). Interestingly, RBPJ itself is also an RA risk gene. Based on GWAS catalog (2019), two LDs on the RBPJ locus are associated with RA. In LD with rs874040, we identified fSNPs rs17630466 and rs35944082 regulating RBPJ expression via binding to LEF1 and MSH6 as well as to IRF2, APOBEC3C, MEF2B and MEF2D, respectively. All these proteins were demonstrated to be the transcriptional suppressors, downregulating RBPJ expression in human B cells and primary human FLS; therefore, more binding of these proteins to the risk allele G of either rs17630466 or rs35944082 will result in a decreased expression of RBPJ. As RBPJ is a transcriptional suppressor of CD40 (8), decreased expression of RBPJ will lead to upregulation of CD40, which, in turn, can activate NF-kB p65 upon CD40L ligation. Activation of NF-kB p65 plays an important role in the pathogenesis of autoimmune diseases, such as RA, MS and SLE (30,31). Together, these data explain the contribution of the RA-associated RBPJ locus to the susceptibility of RA.

Based on our findings in this work as well as in previous publications (6–8), we constructed an RA-associated CD40-induced NF-kB STTRN as shown in Figure 6. In this STTRN, transcriptional regulation of four RA risk genes HADC7, MEF2D, RBPJ and CD40 is functionally linked to the transcription of other eight RA risk genes CD83, CDK4/CDK6, c-FLIP, ICAM1, TNFAIP3, IRF4 and TRAF1 as well as two proinflammatory genes IL6 and IL1β via CD40-induced NF-kB signaling. We verified this STTRN by inactivating HDAC4/HDAC7 either using HDAC inhibitor ITF2357 or using RNAi knockdown of HDAC4/HDAC7. We show that inactivation of HDAC4/HDAC7 can result in decreased expression of MEF2D, which, in turn, downregulates CD40 by suppressing RBPJ expression. Downregulation of CD40 can further lead to the inactivation of NF-kB p65 by dephosphorylation, which suppresses the expression of at least 10 genes as the direct targets of NF-kB p65, including the eight RA risk genes as well as IL6 and IL1β. We also show that inactivation of HDAC4/HDAC7can result in decreased cell migration and proliferation and increased apoptosis, presumably via the STTRN by inactivating ICAM1, CDK4/6 and c-FLIP. Upregulation of proinflammatory cytokines IL6 and IL1β was reported in patients with RA (49,50), and increased cell migration and proliferation and decreased apoptosis in FLS were also reported to be responsible for the synovial hyperplasia that contributes to the destruction of cartilage and bone in RA (43). Therefore, inactivation of IL6 and IL1β, as well as ICAM1, CDK4/6 and c-FLIP expression, by ITF2357 reveals a novel mechanism that shows the pharmacological effect of ITF2357 as an approved drug for the treatment of SJIA. Other possible pharmacological effects of ITF2357 were reported to show that ITF2357 can also decrease disease activity by increasing the number of Treg cells and decreasing the number of TH17 cells in the arthritis mouse model (51). While both MC1568 and ITF2357 inactivate NF-kB p65 via the RA-associated CD40-induced NF-kB STTRN, MC1568 alters the expression of NF-kB p65 target genes in a different pattern, such as upregulation of ICAM1, IL6 and IL1β expression. These differences explain why ITF2357, not MC1568, can be used as a prescribed drug for SJIA.

Figure 6

Scheme showing incomplete RA-associated CD40-induced NF-kB STTRN in human B cells and FLS. Twelve RA risk genes were shown in this RA-associated CD40-induced STTRN, including HDAC7, MEF2D, RBPJ, CD40 as well as CD83, CDK4, CDK6, c-FLIP, ICAM1, TNFAIP3, IRF4 and TRAF1. Only the transcriptional regulation of CD40 and RBPJ was shown in an RA-associated fashion. TRAF1, NFKBIE, TRAF6, TNFAIP3 and REL are the five other RA risk genes involved in CD40-induced NF-kB signal. RA risk gene PAX5 could be another upstream regulator of this STTRN by modulating LEF1 expression. ITF2357, an approved drug for SJIA, suppresses transcription of the eight RA risk genes as well as IL6 and IL1β via HDAC4/HDAC7, MEF2D, RBPJ and CD40 transcriptional regulation. The possible functional relevancies of these eight RA risk genes as well as IL6 and IL1β are listed. TRAF1, CD40 and TNFAIP3 are in the feedback loop of transcriptional regulation via NF-kB p65.

Many questions remain to be answered in this RA-associated CD40-induced NF-kB STTRN, including: (1) what are the regulatory proteins that modulate RBPJ expression through the other four fSNPs in LD with rs874040 on the RBPJ locus? Are these proteins encoded by RA risk genes or regulated by other RA risk genes? (2) How are the other RA risk genes in this STTRN regulated, especially for HDAC7 and MEF2D? The association of HDAC7 with RA and MEF2D with SLE has never been validated; therefore, more functional studies on these risk genes are required to determine their roles in the RA. Once all the RA-associated STTRNs are fully elucidated, the pathogenesis and susceptibility of RA could be predicted in accordance with an individual’s genotype. Besides, as evidenced by our studies with ITF2357, the RA-associated STTRNs will also help develop rational strategies for the selection of therapeutic targets.

In addition, as we showed in this report, we have developed a sequential technique to perform post-GWAS functional studies by coupling SNP-seq with FREP-MS and AIDP-Wb. SNP-seq is an unbiased high-throughput screening technique to identify fSNPs based on a type IIS restriction enzyme Bpm I protection assay (8). Once we confirm the fSNPs by using EMSA and luciferase reporter assay, FREP-MS can be applied to identify proteins that specifically bind to the fSNP by enzymatically separating the fSNP sequence together with its binding proteins from the remaining DNA-nuclear extract (NE)-bead complex (8,52). After the specific binding of a regulatory protein to a fSNP is validated using AIDP-Wb (13) and ChIP assay, RNAi knockdown is performed to assess the function of the fSNP-bound proteins in regulating risk gene expression. While this approach is relatively efficient in identifying and characterizing disease-associated fSNPs, its limitation is obvious. All these techniques are in vitro assays using cell lines or cellular extracts, and they often do not fully represent what is occurring in vivo. Especially for RNAi, it is purely a loss-of-function technique that cannot identify the allele-imbalanced regulatory mechanism as a non-coding fSNP does. In consideration of these premises, we think that using CRISPR/cas9 gene editing to introduce the different alleles of a fSNP in mice would be the best approach to study non-coding fSNPs. However, this approach is currently a time-consuming and technically difficult process; therefore, we believe that a functional validation of a fSNP in vitro is a guarantee for success for in vivo knockin studies.

In summary, we have shown an example of how to utilize post-GWAS functional studies to collect functional information for translating GWAS data into biological mechanisms. Based on these functional data, an RA-associated CD40-induced NF-kB STTRN was identified and validated. As CD40-induced NF-kB signaling is involved in many diseases, applying the same approach to identify each disease-associated, cell-type-specific CD40-induced NF-kB STTRN may help identify unique therapeutic targets for a range of diseases.

Materials and Methods

Cells and culture

The human B cell line BL2 was purchased from DSMZ and cultured in RPMI 1640 medium. The primary FLS were isolated from synovial tissue samples and cultured in DMEM medium. Both media were supplemented with 10% fetal bovine serum.

Fibroblast-like synoviocytes

Synovial tissue samples were collected from OA patients undergoing knee arthroplasty. The use of human materials was approved by University of Pittsburgh with the IRB number: STUDY18100138, and written informed consent was obtained from all individuals before the operative procedure. FLS were isolated from the synovial tissue samples according to a protocol previously described (53).

Primers and antibodies

All primers used in this work are listed in Supplementary Material, Table S4 and were purchased from IDT, except for CD40. The CD40 primers were purchased from Genecopoeia with Cat#: HQP022955. All antibodies used are listed in Supplementary Material, Table S5 with the vendor’s information.

SNP-seq

SNP-seq was performed as previously described (8). The SNP-seq library was made by the construct sequence shown in Supplementary Material, Table S4. The library was generated by Ultramer® DNA oligonucleotides (ITD) free of mutations. For fSNP screening, 100 ng of library DNA was amplified by PCR using primers bioseq and G3 for 15 cycles with AccuPrime Taq (Thermo Fisher Scientific) at 95°C for 90 s; 58°C for 90 s and 72°C for 40 s. After gel purification, 10 ng of biotinylated DNA was attached to 4 μl streptavidin-Dynabeads (Invitrogen) according to the manufacturer’s protocol. The DNA-beads were then incubated with 100 μg NE for 1 h at room temperature in LightShift Chemiluminescent EMSA Kit reaction buffer (Thermo Fisher Scientific). After washing and separation, the DNA-NE-beads were digested with 2 μl Bpm I (NEB) for 30 min at 37°C. After another wash and separation, the DNA was amplified again with bioseq and G3 and re-attached to the Dynabeads for the next SNP-seq cycle. Seven cycles were performed in total with four buffer-treated controls and four NE-treated samples. Next-generation sequencing library was prepared using DNAs from cycles 1, 4 and 7 (8). Sequencing was performed using Miseq for 50 bp single end sequencing.

For quality control, we eliminated the SNPs, from which complete sequencing data were not available across cycles 1, 4 and 7, and we also eliminated the SNPs that carry mutations in the Bpm I binding sites. To identify fSNPs, we first calculated the ratio of the sequence count between the two alleles for each SNP and for each of the four replicates at cycle 7. We identified SNPs that demonstrate allele-imbalanced protection at cycle 7, with a significant difference in the ratio between the four controls and the four samples using Student’s t-test and with a P-value < 0.05. Second, we identified SNPs for which the allele-imbalanced protection increased with the increasing cycle number, indicating progressive enrichment. For each SNP, we calculated the average ratio of the sequence count between alleles from the four buffer-treated controls and the four NE-treated samples and normalized the ratio of the sample with the control. We next used the normalized ratio from cycles 1, 4 and 7 to calculate for Slope and identified SNPs with a Slope > 0.05 and <−0.05 as an empirical cut-off point. We then identified candidate fSNPs with a P-value < 0.05 and a Slope > 0.05 and <−0.05.

Electrophoretic mobility shift assay

EMSA was performed using the LightShift Chemiluminescent EMSA Kit (Thermo Fisher) according to the manufacturer’s instructions. For the probe, a 31 bp SNP fragment with the SNP centered in the middle was made by annealing two oligos. The double-stranded oligos were then biotinylated using the Biotin 3′ End DNA Labeling Kit (Thermo Fisher). The data represent three biological repeats.

FREP-MS

FREP-MS was previously described (8). In brief, ~10 μg of the purified FREP construct DNA was conjugated to 150 μl streptavidin coupled Dynabeads (Life Technologies) according to the manufacturer’s instructions. The DNA-beads were then washed and mixed with 1 mg NE in a buffer from the LightShift™ Chemiluminescent EMSA Kit at RT for 1 h. After separation and washing, the DNA-NE-beads were digested with 5 μl EcoR I (100 units/μl NEB) at 37°C for 30 min to remove the 3′ DNA plus the proteins that bound to this non-SNP region. After separation and washing, the DNA-NE-beads were subsequently digested with 5 μl BamH I (100 units/μl NEB) at 37°C for 45 min to release the fSNP sequence plus the fSNP-bound proteins. The supernatant was collected for protein complex identification by mass spectrometry.

Luciferase reporter assay

Luciferase reporter assays were performed using pGL3 luciferase reporter basic vector (Cat#: E1751, Promega). Luciferase activity was measured by the Dual-Glo® Luciferase Reporter Assay System (Cat#: E2920, Promega). All the experiments were performed according to the manufacturer’s protocol. Insert target sequences are listed in Supplementary Material, Table S4. All data represent a combination of six independent biological replicates.

Western blot

Whole cell proteins were isolated with RIPA buffer (Sigma). Cytosolic proteins and nuclear proteins were isolated with NE-PER Nuclear and Cytoplasmic Extraction Reagents (Thermo Scientific) according to the manufacturer’s instructions. Western blots were performed using standard protocol. For control, α-tubulin was used unless as indicated. The data represent three biological replicates.

AIDP-western blot

NE was isolated as described above. A 31 bp biotinylated double-stranded DNA containing either one of the two alleles for an SNP was generated from biotinylated primers purchased from IDT by annealing. Two micrograms DNA was attached to 20 μl Dynabeads™ M-280 Streptavidin according to the manufacturer’s instructions. DNA-beads were mixed with ~100 μg NE at RT for 1 h with rotation in 50 μl binding buffer used in the LightShift Chemiluminescent kit. After magnetic separation and wash, DNA-bound proteins were eluted by adding 1× protein sample buffer and were incubated at 95°C for 5 min for western blots analysis. PARP1 was used for internal loading control. The data represent three biological replicates.

ChIP assay

ChIP was performed as described previously (54). Briefly, scrambled shRNA BL2 cells and shRNA knockdown BL2 cells were cross-linked with 1% formaldehyde for 10 min. Sonication was carried out at 30% amplitude with 20 s on and 50 s off for 5 min. Ten micrograms antibody coupled to Dynabeads™ Protein A/G (Thermo Fisher Scientific, Cat#:10001D and 10003D) was used for immunoprecipitation. DNA was purified with Qiagen PCR purification kit after reversal of the cross-link. Rabbit IgG was used as an isotype control antibody (Cell signaling Technology, Cat#: 2729). The purified DNA was used for real-time PCR analysis. The primers that were used are listed in Supplementary Material, Table S4. The data represent three biological replicates.

Real-time PCR

Total RNA was isolated with RNeasy Mini kit (Qiagen). cDNA was synthesized with SuperScript® III Reverse Transcriptase (Invitrogen) after 1 μg RNA sample was treated with DNase I (Invitrogen). All of the procedures were performed following the manufacturer’s protocols. GAPDH was used as a control for all the PCR reactions. All the data represent a combination of three biological replicates unless as indicated.

RNAi knockdown

For shRNA stable knockdown in human BL2 cells, the shRNA lentiviruses were generated by cloning the targeting sequence into pLKO.1 vector (Addgene, Cat#: 10878) according to Addgene’s pLKO.1 protocol. For siRNA transient knockdown in human FLS, siRNAs were purchased from Horizon Discovery, and the knockdown was performed according to the manufacturer’s protocol.

MC1568 and ITF2357 treatment

BL2 cells or FLS were seeded 1 day before the treatment. On day 1, both DMSO (control) and different concentrations of MC1568 (5 and 10 μm) or ITF2357 (0.5 and 1 μm) were added. For immunoblot detection, after 24 h treatment, BL2 cells or FLS were activated either with or without MEGACD40L (64 ng/ml) for 15 min.

Scratch-wound assay

FLS were seeded into a 24 well tissue culture plate at a density of ~60–70% confluence 1 day before the treatment. Then, cells were either treated with 1 μm ITF2357 or were transfected with siRNAs targeting both HDAC4 and HDAC7. After another 24 h incubation, cells were activated by MEGACD40L at 64 ng/ml for 15 min. The monolayer cells were then gently and slowly scratched with a 200 μl pipette tip across the center of the well. After scratching, the detached cells were gently removed by washing and the well was replenished with fresh medium. After another 18 h incubation, photos were taken. The data represent three biological replicates.

Cell proliferation analysis

1X10 (5) FLS were seeded in a six well plate 1 day before the treatment. Cells were then either treated with 0.5 and 1 μm ITF2357 or were transfected with siRNAs targeting both HDAC4 and HDAC7 for 24 h. Cells were then further activated by MEGACD40L at 64 ng/ml for 15 min and were stained with 20 μl of 500 μg/ml propidium iodide (Biolegeng, Cat#: 421301) according to the manufacturer’s protocol. Data were collected by Attune™ NxT Flow Cytometer Software (Thermo Fisher Scientific) and FlowJo version 6. The data represent three biological replicates.

Caspase-Glo® 3/7 assay

1X10 (4) FLS were seeded in a 96 well plate 1 day before the treatment. Then, cells were either treated with 0.5 and 1 μm ITF2357 or were transfected with siRNAs targeting both HDAC4 and HDAC7. After 24 h incubation, cells were further activated by MEGACD40L at 64 ng/ml for 15 min and were induced for apoptosis by exposing the cells to anti-Fas IgM mAb (MyBioSource Cat#: MBS485033MBL) at 1 μg/ml for another 24 h as previously described (48). Apoptosis was evaluated by detecting the caspase 3/7 activities using Caspase-Glo® 3/7 assay system (Cat#: G8091, Promega) according to the manufacturer’s protocol. Fluorescence was read at 485 nm excitation, 520 nm emission, followed by luminescence. The data represent a combination of three biological replicates.

Statistical analysis

P-values were calculated using Student’s t-test with two tails. Error bars represent the median with standard error (S.E.)

Click here for additional data file.