Overall Downregulation of mRNAs and Enrichment of H3K4me3 Change Near Genome-Wide Association Study Signals in Systemic Lupus Erythematosus: Cell-Specific Effects

This study was designed to define gene expression and H3K4me3 histone modifications in T cells, B cells, and monocytes in systemic lupus erythematosus (SLE). Array studies of total peripheral blood mononuclear cells have demonstrated gene expression signatures related to neutrophils, interferon, and other inflammatory pathways. It is not clear how consistent these effects are across different cell types. In this study, RNA-seq and chromatin immunoprecipitation-seq were utilized to identify gene expression patterns and H3K4me3 histone modifications related to promoter activation in SLE. Across the three cell types, there was 55% concordance for gene expression changes related to SLE. Key conserved pathways were ribosome biogenesis among upregulated genes and heat shock response among downregulated genes. ETS family transcription factors (TFs) and STAT1 were revealed as common regulators by position weight matrices. When epigenetic changes were leveraged with gene expression, the pivotal TFs ATF3 and FOS were defined with ATF3 also cross-referencing with gene expression-identified TFs. Genome-wide association study (GWAS) single nucleotide polymorphisms associated with SLE were cross-referenced with both mRNA and H3K4me3 changes in SLE. Baseline mRNA expression and H3K4me3 peak height was higher at sites that cross-referenced with GWAS signals, however, all three cell types exhibited an overall decrease in expression of GWAS-associated RNAs differentially expressed in SLE. H3K4me3 changes in SLE were also enriched in GWAS-associated sites. In summary, the SLE disease process is associated with both shared and cell-specific changes in gene expression and epigenetics. Surprisingly, GWAS-associated RNAs were overall markedly decreased across all three cell types. TF analysis identified ATF3, FOS, STAT1, and ETS family members as critical, all pathways with a recognized relationship to the SLE disease process. GWAS signals clearly mark both cell-type specific changes in SLE as well as concordant changes across all three cell types. Interpretation of single nucleotide polymorphism effects in SLE will require tissue-specific mechanistic studies and therapeutics will require mechanistic studies in multiple cell types.

This study was designed to define gene expression and H3K4me3 histone modifications in T cells, B cells, and monocytes in systemic lupus erythematosus (SLE). Array studies of total peripheral blood mononuclear cells have demonstrated gene expression signatures related to neutrophils, interferon, and other inflammatory pathways. It is not clear how consistent these effects are across different cell types. In this study, RNA-seq and chromatin immunoprecipitation-seq were utilized to identify gene expression patterns and H3K4me3 histone modifications related to promoter activation in SLE. Across the three cell types, there was 55% concordance for gene expression changes related to SLE. Key conserved pathways were ribosome biogenesis among upregulated genes and heat shock response among downregulated genes. ETS family transcription factors (TFs) and STAT1 were revealed as common regulators by position weight matrices. When epigenetic changes were leveraged with gene expression, the pivotal TFs ATF3 and FOS were defined with ATF3 also cross-referencing with gene expression-identified TFs. Genome-wide association study (GWAS) single nucleotide polymorphisms associated with SLE were cross-referenced with both mRNA and H3K4me3 changes in SLE. Baseline mRNA expression and H3K4me3 peak height was higher at sites that cross-referenced with GWAS signals, however, all three cell types exhibited an overall decrease in expression of GWAS-associated RNAs differentially expressed in SLE. H3K4me3 changes in SLE were also enriched in GWAS-associated sites. In summary, the SLE disease process is associated with both shared and cell-specific changes in gene expression and epigenetics. Surprisingly, GWAS-associated RNAs were overall markedly decreased across all three cell types. TF analysis identified ATF3, FOS, STAT1, and ETS family members as critical, all pathways with a recognized relationship to the SLE disease process. GWAS signals clearly mark both cell-type specific changes in SLE as well as concordant changes across all three cell types. Interpretation of single nucleotide polymorphism effects in SLE will require tissue-specific mechanistic studies and therapeutics will require mechanistic studies in multiple cell types. inTrODUcTiOn Systemic lupus erythematosus (SLE) is the quintessential systemic autoimmune disease. Nearly, every organ can be involved in SLE and a great diversity of immunologic features have been described (1). In spite of enormous effort to identify a common pathway that links all clinical and laboratory features, only the interferon signature has been identified consistently across multiple laboratories (2)(3)(4). This major advance in the understanding of SLE arose from the study of peripheral blood mononuclear cells (PBMC) by gene expression arrays and is associated with a broad range of autoimmune diseases (5)(6)(7). Most cells express type I interferon receptors, although transcriptional responses differ between cell types. Thus, responses to interferons could be coordinate among different cells or discordant between cell types.
Other cytokines over-expressed in SLE have receptors that are not ubiquitously expressed. Furthermore, we have identified high levels of circulating endotoxin in SLE which could drive immune activation in a cell-specific manner due to the limited expression of TLR4 (8). We therefore hypothesized that different cell types in peripheral blood are impacted by the SLE disease process in both coordinate and distinct manners and key commonalities could be dissected by using a joint transcriptomic and epigenetic approach in purified cell types. Recognition of distinct cell effects is critical for interpretation of genome-wide association studies (GWASs) where efforts to examine mechanism may require a tissue-specific approach. Several prior studies have used array approaches on purified hematopoietic cells. In these studies, concordance across cell types was modest (9)(10)(11). Interferon response genes exhibited increased expression as a category in all sorted cell types but gene-level expression was often tissue specific. Using RNA-seq, deconvolution was applied, finding that MHC gene expression was most altered in B cells and ribosomal gene expression was most altered in monocytes (12). Thus, there has been little concordance across cell types in different studies. These expression studies suggest that many SLE effects are tissue specific; however, these studies have generally used arrays with analyses restricted to the genes included on the array. RNA-seq offers a more holistic approach and has the advantage of a greater dynamic range.
A compelling reason to focus on epigenetic changes is that these changes are typically more stable in patients over time compared with gene expression. Although epigenetic changes in SLE have received attention in recent years, most studies have focused on DNA methylation (13)(14)(15)(16). Changes in histone modifications are restricted to tighter peaks than that seen with DNA methylation and can be more easily used to focus on DNA motifs associated with transcription factor (TF) binding. We have previously found IRF1 as a pivotal TF by defining H4ac peaks in monocytes from SLE patients (17)(18)(19). We went on to define changes in H3K4me3 at promoters and enhancers in SLE monocytes, linking a type I interferon effect to chromatin and defining a set of TFs identified through position weight matrices (PWMs) (20,21). There are additional studies that have utilized ATAC-seq of SLE B cells, which found that increased accessibility was seen at loci related to B cell activation (22). To best define cell-type specific changes in the transcriptome and epigenome in SLE, we chose to examine changes in H3K4me3, a histone modification associated with active promoters. We paired this with RNA-seq to define changes in gene expression. Finally, we cross-referenced our SLE effects on the transcriptome and epigenome with published GWAS variants. Our findings support a model where many SLE effects are concordant with cell-type specific differences in H3K4me3 and enrichment of GWAS signals. These data have important implications for mechanistic studies of single nucleotide polymorphisms and the development of therapeutics.

MaTerials anD MeThODs samples and cells
The female healthy donor blood samples were obtained from the Center for AIDS Research, under a protocol approved by the University of Pennsylvania Perelman School of Medicine Institutional Review Board (IRB) in accordance with their guidelines. All subjects gave written informed consent in accordance with the Declaration of Helsinki. Samples were handled identically as the shipped SLE samples. The female SLE patient blood samples were from the Johns Hopkins Lupus Cohort (23,24) under a separate IRB protocol in accordance with the recommendations and guidelines of the Johns Hopkins IRB Committee. All patients gave written informed consent in accordance with the Declaration of Helsinki. Six SLE patients and six controls were used ( Table 1). PBMC and monocytes were purified using adherence as previously described (8). Dynal beads (Dynabeads, Invitrogen, Oslo, Norway) were used to separate CD3 T cells and CD19 B cells according to the manufacturer's protocol. Flow cytometry was used to ensure purity of the three cell populations.
Peripheral blood mononuclear cells were stimulated with 3 µg/ml of phytohemagglutinin (PHA) for 24 h followed by cell separation as above. The MYC inhibitor, 10058-F4, was used at a concentration of 50 µM. The inhibitor led to minimal cell death at this concentration and was added 20 min prior to stimulation with PHA.

rna isolation and rna-seq
Total RNA was isolated from isolated T cells, B cells, and monocytes by the Qiagen RNeasy Kit (Valencia, CA, USA) and DNA was removed by column DNase digestion. The RNA quality was defined by RIN and OD260/280. The RNA-seq libraries were made with the Ovation ® Ultralow Library Systems and sequenced on an Illumina HiSeq at BGI@CHOP. A customized workflow was utilized. We aligned sequence reads to the human reference genome (GRCh38) and transcriptome using the STAR program. Read counts of each cell type were converted to log2 scale and normalized by the Lowess method. The limma method was applied to the normalized data to obtained differentially expressed genes.

chromatin immunoprecipitation (chiP) and chiP-seq
Chromatin immunoprecipitation assay experiments were performed as previously described (20,25). The antibody for H3K4me3 was from Active Motif (Carlsbad, CA, USA). The Illumina TruSeq ChIP library preparation kit (San Diego, CA, USA) was used and sequencing was performed on an Illumina HiSeq. The reference genome (hg38) was indexed by the novoindex function of the NovoAlign package. Reads were filtered by SAM fields such as mapq, cigar, and flag and extended to 200 bp long at their 3′ end. Reads were mapped to the 1 kb promoter regions of known genes to get an integer matrix from each cell type. The read count matrices were normalized and analyzed the same way as the RNA-seq data to obtain results of differential H3K4me3. gene set analysis Differential expression or differential H3K4me3 of genes was ranked by their significance using the limma test. Parametric analysis of gene set enrichment (PAGE) method was used to compare each gene set and all the other genes. PAGE reported the significance of the difference with a p-value, false discovery rate, and enrichment score. Positive and negative enrichment scores indicated the overall increase and decrease, respectively, of gene expression or H3K4me3 of the gene set.

TF analysis
2,414 PWMs of TFs from multiple databases (ENCODE, Jasper, TRANSFAC, and UniPROBE) were used as the reference dataset. To identify TF targets, we mapped promoter sequences of known genes to these PWMs and selected genes with the best matches to a PWM as potential target genes regulated by the corresponding TF. On average, each target set includes approximately 1,300 genes. We ran gene set analysis to identify target sets with overall change of expression or H3K4me3 in SLE.

genome-Wide association crossreferencing
Chromatin immunoprecipitation-seq data processed as above were used to analyze enrichment around SNP sites. The average coverage around a 1 kb (−500 to +500) region of SNPs was calculated. The log-transformed coverage was normalized by Loess method between all samples within each cell type. RNA-seq data were processed as previously described. SNPs were mapped to nearest gene TSS for RNA-seq analysis. Gene-level read counts were normalized and group comparison between controls and patients were also done as previously described. 95 SNP with GWAS p-values less than 1E-8 were selected from a study of 1,311 SLE patients (26). Additionally, a set of 46 SNPs (27) curated from two published studies was used (28,29) and this is referred to as LD46 and a set of 25 SNPs curated from a meta-analysis of 7,219 SLE cases was used and referred to as META25 (28). All three sets of SNPs were used to generate a joint set of 143 SNPs used for analyzing RNA-seq data. The ChIP-seq data were analyzed using all SNPs and binning according to the p-value of the association with SLE.
Online visualization of data is provided at http://projectsle2. awsomics.org.

Transcriptomic analysis
There is relatively little known about the effect of SLE in different cell types and whether the interferon effect is common to all hematopoietic cell types. The goal of this transcriptomic analysis was to identify genes with differential expression between control and SLE groups that were concordant in CD3 T cells, CD19 B cells, and monocytes. Control and SLE samples clearly clustered distinctly in a principal component analysis, demonstrating a clear effect of disease (Figure 1). The direction of change for SLE compared with control was the same for the three cell types. Coding genes and long non-coding RNAs exhibited similar patterns while repetitive elements were different.

Differential Expression in Each Cell Type
In CD19 B cells, 1,589 and 1,489 genes had higher and lower expression, respectively in SLE (at least twofold change and p < 0.05). "Ribosome biogenesis" (GO:0042254) was the most upregulated gene set and the MSigDB "TNFa signaling via NFKB" pathway was the most downregulated gene set. In CD3 T cells, 1,313 and 1,173 genes had, respectively, higher and lower expression in SLE. "Ribosome biogenesis" and "TNFa signaling via NFKB" pathway were also the most up-and downregulated gene sets, respectively, as was seen in B cells. In monocytes, 848 and 716 genes had higher and lower expression in SLE, respectively. "Ribosome biogenesis" was again the most upregulated gene set and "nuclear chromatin" (GO:0000790) was the most downregulated gene set.
From this initial comparison, two findings stand out. (1) Upregulated and downregulated gene sets were quantitatively balanced in all three cell types. (2) Surprisingly, the gene set "ribosome biogenesis" was identified as upregulated in all three cell types while heat shock protein genes were uniformly found to be downregulated in each cell type. Individual genes differed; however, the signatures were robust across cell types.

Concordance of Differential Expression Among the Three Cell Types
Genes differentially expressed in SLE were likely to be shared between cell types. On average, genes differentially expressed Each SLE-control cell type pair exhibits a uniform shift in the PC2 axis from control to SLE while the PC1 axis largely discriminates cell type for coding genes and lncRNAs. Repetitive elements still exhibit clear distinctions between cell type and disease state but the direction of change is different from the coding genes. H3K4me3 has a strong cell type dependence but less discrimination between disease and control samples.
in one cell type were 6.1 times more likely to be differentially expressed toward the same direction and 4.9 times less likely to be differentially expressed toward the opposite direction in the other two cell types. Respectively, 95 and 97 genes had higher and lower expression concordantly in all cell three types. GTF2E1 (general TF) was the most consistently upregulated gene, increased by 8.5 (B cell), 4.1 (monocyte), and 6.1 (T cell) fold in SLE. HSPA1A (the gene encoding the heat shock protein, Hsp70) was 197 (B cell), 33.5 (monocyte), and 17.6 (T cell) fold decreased, making it the most consistently downregulated gene.
To better examine the concordance and cell-type specific expression, we performed clustering. Clustering was performed across all 18 RNA-seq samples to identify 8 gene clusters having the distinctive co-expression patterns across cell types (Figure 2). The two largest clusters (Cluster 1, downregulated and Cluster 8, upregulated) were concordant across all three types and accounted for 55% of the total differential gene set. One of the gene sets most over-represented in Cluster 1 (downregulated concordantly) was "MAPK signaling pathway. " Genes related to major biological processes regulated by p38 MAP kinases including cytokine response, differentiation, apoptosis, and autophagy, were also downregulated in all three cell types (Figure 3; Table S1 in Supplementary Material). Negative regulators of MAP kinases, DUSP4/11/12, all had significantly higher expression in all three cell types. Top gene sets over-represented in Cluster 8 (upregulated concordantly) were related to RNA processing, transport, and translation.

Gene Set Analysis
We noted that the interferon signal was barely detectable as defined by gene set enrichment. Because it was surprising that the interferon signature was so subtle, we sought gene expression for interferon-responsive genes and tabulated levels in each cell type. When this was performed, we saw distinctive effects in each cell type (Figure 4). In B cells, only CCL2 was over-expressed in SLE compared with controls with a significant p-value. In T cells, IFIT1, IFIT3, IFITM1, OASL, and RSAD2 were overexpressed in SLE compared with controls. In monocytes, CCL2 and OASL were over-expressed in SLE compared with controls. These over-expressed genes were balanced by a group that was under-expressed in SLE compared with controls. In contrast to the interferon response, heat shock proteins were downregulated comparably in all three cell types (Figure 4).
The most consistently upregulated gene set was ribosome biogenesis. Ribosomes increase with proliferation, cell activation, and enlargement. MYC was also a common pathway among upregulated genes in all three cell types (22-fold increased in T cells, 39-fold increased in B cells, and 113-fold increased in monocytes) and MYC is known to regulate ribosome mass and protein synthesis (30). We therefore examined known cell typeindependent MYC targets in our gene set (31). MYC targets were significantly increased in T cells and B cells (p < 0.0001) although other oncogenes such as MYB and TP53 were not consistently upregulated in our samples ( Figure 5A). MYC has been previously described as upregulated in SLE and in murine models of SLE (32)(33)(34).
To test the relationship between MYC and ribosome biogenesis, PBMC from healthy donors were stimulated with PHA to induce proliferation and the elaboration of cytokines that model immune activation in SLE ( Figure 5B). Cells were then separated and qRT-PCR performed for ribosome biogenesis targets that were identified in the RNA-seq data. PHA stimulation increased the ribosome biogenesis targets, as expected. The MYC inhibitor significantly depressed expression of the ribosome biogenesis targets after PHA stimulation in all three cell types (p < 0.01). These data support a role for MYC in the expression of ribosome biogenesis genes and provides an important conceptual link between two major gene sets with altered expression in SLE.

h3K4me3 analysis
This study was designed to leverage both expression and epigenetics to develop a more complete picture of the pathways dysregulated in SLE and their cell-type concordance. We examined histone marks as a more durable echo of the cells' exposures and milieu. We employed ChIP-seq for H3K4me3, a mark of active promoters.
We first classified promoters to analyze changes in SLE based on potential for expression. Based on the average H3K4me3 in  controls and the level of expression defined by RNA-seq, promoters were classified into three classes ( Figure 6A): (1) The "inactive" promoters had low-H3K4me3 read depth (the left peak), The "active" promoters had high-H3K4me3 read depth (the right peak), and (3) The "poised" promoters had intermediate H3K4me3 read depth. In all three cell types, poised promoters were most likely to have changed H3K4me3 in SLE and none of the promoters changed dramatically from "inactive" to "active" or vice versa ( Figure 6B). Comparably to what was seen in the transcriptomic data, more genes had decreased H3K4me3 than   increased H3K4me3. These data support a model where cells exhibit altered behavior in SLE but have not been "re-wired" to execute completely novel functions.

Differential H3K4me3 in Each Cell Type
Gene set analysis was applied to the change of promoter H3K4me3 in SLE. It identified "cell fate commitment" (GO:0045165) and  "cell chemotaxis" (GO:0060326) as the gene sets with the most increased and decreased H3K4me3, respectively, in B cells. In CD3 T cells, MSigDB "TNFA_SIGNALING_VIA_NFKB" pathway was the gene set with the most increased H3K4me3 and "lymphocyte migration" (GO:0072676) was the gene set with the most decreased H3K4me3. We noted that the "TNFA_SIGNALING_ VIA_NFKB" pathway was downregulated at the level of RNA and upregulated at the level of promoter H3K4me3 in T cells. This could be interpreted as a chromatin environment favorable for the expression of NFκB targets but lacking the acute stimulation for expression or as a state consistent with endotoxin tolerance.

Concordance of H3K4me3 Among the Three Cell Types
H3K4me3 in controls agreed on the classification of over 75% of promoters across the three cell types. Genes with the highest H3K4me3 in a single cell type were associated with cell lineage. For example, genes with the highest promoter H3K4me3 in B cells and T cells were enriched with genes of B cell activation (GO:0042113) and T cell activation (GO:0042110), respectively. Genes with highest H3K4me3 in monocyte were enriched with those related to inflammatory response (GO:0006954) and leukocyte migration (GO:0050900). This analysis demonstrates the fidelity of our ChIP-seq approach by confirming the expected cell type specificity. There were not many promoters having significant H3K4me3 change in all three cell types, with B cells and monocytes being weakly correlated (correlation coefficient = 0.02). Promoters with increased H3K4me3 in all cell types included HLAB, RFX2 (regulatory factor X2), and RARA (retinoic acid receptor). Promoters with decreased H3K4me3 in all cell types included KANSL1 (KAT8 regulatory complex), CX3CR1 (chemokine receptor), and RORC (RAR-related orphan receptor).
We examined H3K4me3 using a clustering approach, as was done for the RNA-seq data (Figure 7). Clusters 1 and 8 had decreased and increased H3K4me3, respectively, in all three cell types. The gene set over-represented most significantly in Cluster 1 (downregulated concordantly) was "regulation of leukocyte mediated immunity" (GO:0002703), which includes genes such as CCR2, IL27RA, STAT5A, STAT6, and VAMP2. The MSigDB "TNFA_SIGNALING_VIA_NFKB" pathway was the gene set most significantly over-represented in Cluster 8 (upregulated concordantly). Clusters 2 and 7 were the two largest clusters, with, respectively, decreased and increased H3K4me3 in both monocyte and T cells. Cluster 2 was over-represented with genes of signaling pathways, such as mTOR, Wnt, and PPAR, and Cluster 7 was most over-represented with DNA repair genes.

Concordant Change in Transcription and H3K4me3
We asked how much the change of promoter H3K4me3 contributed to the differential expression of genes in SLE. Both RNA and H3K4me3 changes related to SLE were likely to be shared across all three cell types (Figures 8A,B). Genes with significant differential expression were more likely to have significant H3K4me3 change toward the same direction in all three cell types (Figure 8C), and less likely to have H3K4me3 change toward the opposite direction. This is consistent with the role of H3K4me3 in gene regulation. The gene clusters defined by transcription patterns were also significantly overlapped to their counterparts defined by H3K4me3 patterns (Figure 8D). However, less than 10% of all genes with significant differential expression in SLE also had significant differential H3K4me3 toward the same direction, suggesting that the differential expression of most genes was not directly associated with H3K4me3 change,

Concordant Change in Gene Sets Defined by Transcription and H3K4me3
Our goal was to leverage this robust data set with both chromatin and RNA level results to understand commonalities and differences between cell effects in SLE. We analyzed themes related to cell biological processes since individual genes and promoters might be differentially regulated in the three cell types. Based on gene set analysis, genes related to RNA transport (KEGG:hsa03013) and leukocyte activation (GO:0045321) had increased and decreased transcription and H3K4me3, respectively, in all three cell types. 57% of the RNA transport genes with increased transcription in B cells also had increased promoter H3K4me3. Table S2 in Supplementary Material lists additional gene sets with concordant and discordant transcription and H3K4me3 changes. We considered whether the effects on histone methylation could reflect dysregulated expression of histone methyltransferases or demethylase enzymes. To this end, we examined a select set of histone methyltransferases and demethylases ( Figure 9A). The expression of PRDM9 was the most dysregulated in SLE; however, alterations were not consistent across cell types. When we examined the overall levels of transcripts for the three cell types using the entire set of enzymes, there was no clear imbalance in transcription of methyltransferases or demethylases (Figure 9B).
Our analysis also identified region-specific effects. There were genomic regions within which genes had concordant transcription and promoter H3K4me3 changes. For example, a ~300 kb region on chromosome 17 had decreased RNA and H3K4me3 in B cells from SLE patients. The CD300 gene family resides in this region. In SLE monocytes, there was a decrease in both transcription and H3K4me3 within a ~400 kb region on chromosome 7, a region that includes members of GIMAP gene family. GIMAP genes are involved in lymphocyte development. Gene variants of GIMAP family members have been associated with atopic and autoimmune disease (35).

TF Motif Analysis
In addition to analyzing gene sets and region effects, we analyzed upstream transcriptional activators, reasoning that the individual genes induced or repressed in the SLE disease process may differ from cell type to cell type but that there might be common signaling pathways and TFs amenable to drug targeting. Gene set analysis of TF targets revealed strong agreement between the cell types using the RNA-seq dataset. Most noticeably, many members of the ETS (E26 transformation specific) TF family had higher target expression in all cell types (Figure 10A). FigUre 8 | Overlapping of genes with significant differential expression or H3K4me3 in systemic lupus erythematosus (SLE). Genes with differential expression or differential H3K4me3 were likely to be shared across all three cell types. (a) Overlapping of genes with differential expression in multiple cell types. (b) Overlapping of genes with differential promoter H3K4me3 in multiple cell types. Concordance of expression and H3K4me3 was defined by cell type and cluster. (c) Overlapping of genes with differential expression and promoter H3K4me3 in each cell types. (D) Overlapping of gene clusters grouped by expression and H3K4me3 data. Odds ratios and labeled p-values were calculated by the Exact test. We compared TF targets using the ChIP-seq dataset as well after adjusting data for the GC content of their promoter sequences. Target sets of TFs such as SP1, ATF3, and HMGA1 had increased H3K4me3 in all cell types (Figure 10B). We wished to determine if there was a common pathway linking H3K4me3 and gene expression. ATF3 was identified as statistically associated across all cell types in both RNA-seq and ChIP-seq data sets (Figure 11). CDK11A, a kinase gene involved in apoptosis, is a target of ATF3 and had at least 60% increase of both expression and H3K4me3 in all three cell types, supporting a functional role for the ATF3 identified by the in silico approach.

genome-Wide association-identified single nucleotide Polymorphisms associated With sle analyzed
In addition to gene expression studies identifying the interferon signature, GWAS has been instrumental in identifying pivotal pathways in SLE such as immune complex processing, type I interferons, and lymphocyte signaling (1,36). Some identified variants have been analyzed with respect to functional differences, however, most GWAS variants are found in regulatory regions not coding regions, thus hampering functional analysis. Additionally, cell-specific effects may complicate efforts to define functional consequences of variants. We therefore crossreferenced our RNA-seq and ChIP-seq data sets with published GWAS data using three sets of variants. We first hypothesized that the GWAS variants would be enriched at locations with higher RNA expression. This should be reflected in higher RNA and H3K4me3 signals near the variants. We compared expression from control cells near the combined set of 143 GWAS variants to all other genes ( Figure 12A). All three cell types have higher expression from the TSS nearest the GWAS variants, supporting a functional role for the variants. Similarly, H3K4me3 was higher surrounding the GWAS variants with the highest p-values for SLE association (displayed on the X-axis) and the two curated sets of variants, LD46 and Meta25, again supporting functionality ( Figure 12B).
We further reasoned that SLE-associated genomic variants might be further enriched with changes in expression related to SLE or changed H3K4me3 related to SLE. We therefore looked at change in expression and change in H3K4me3 near the GWAS variants. Near the GWAS variants, there was an impressive decrease in RNA abundance in all three cell types ( Figure 13A). The change in H3K4me3 near GWAS variants was similarly analyzed and monocytes exhibited decreased H3K4me3 while T cells and B cells exhibited an increase in H3K4me3 near GWAS variants (Figure 13B). The SNPs with the highest p-value (displayed on the X-axis) have the strongest association. These data are consistent in supporting the concept that GWAS variants impact gene regulation.

DiscUssiOn
To understand whether conserved pathways affect multiple cell types in SLE, we performed a detailed analysis of RNA and H3K4me3 changes in SLE. We selected patients with low-disease activity to focus on core facets of the SLE transcriptome and to avoid medication effects. 55% of the genes with altered expression were concordant across the three cell types. The most pivotal modules of gene expression identified in all three cell types in our study were completely fundamental to cell biology. RNA-seq analysis identified MAP kinase and NFκB pathways as strongly downregulated in SLE. Diminished MAP kinases have been previously described in SLE T cells (37,38). We also identified ribosome biogenesis as a concordantly upregulated gene set and linked it to the MYC pathway. Increased MYC expression has been seen in SLE and is known to upregulate ribosome biogenesis (30,(32)(33)(34)39). MAP kinases promote apoptosis via MYC (40)(41)(42). Thus, the lower levels of MAP kinase transcripts in SLE cells would be predicted to have the effect of promoting cell survival, growth, and proliferation, all settings where ribosome biogenesis is increased ( Figure S1 in Supplementary Material). A curious finding was enrichment of NFκB-regulated genes among downregulated RNAs but increased H3K4me3. NFκB-dependent transcriptional activation relies on H3K4 methyltransferases such as MLL and SET families recruited by MKL1 (43). Thus, the pattern we identified is most consistent with a state of endotoxin tolerance (44).
Our ultimate goal was to define pivotal modulators of gene expression. We therefore examined upstream TF motifs and pathways using an in silico approach. The TFs implicated in concordant transcription changes in all three cell types were ETS family members and STAT1, among others. These two TFs have been highlighted because of their recognized role in SLE (45)(46)(47). ATF3 was notable because it reached significance when analyzed across all cell types and in both RNA-seq and ChIP-seq data sets and ATF3 targets exhibited increased expression. ATF3 is inducible by type I interferons in macrophages, providing a potential mechanism for increased function in SLE (48)(49)(50).  ATF3 has been ascribed roles as diverse as tumor suppression, development, and cell survival, however, little is known regarding its role specifically in SLE.
We found a limited signal of type I interferon effects with diverse genes impacted across the three cell types. The type I interferon signature appears to be less robust in purified cell populations studied by other groups as well (9,10). In our study, this could be due to our focus on patients with low-disease activity, ascertainment of patients with autoantibody profiles that are less associated with type I interferon signatures (51,52), or an effect related the separation of cell populations and removal of low-density neutrophils.
We analyzed our data from the perspective of known variants associated with SLE because they have been postulated to act primarily at the level of gene regulation. Cross-referencing of SNP loci with DNase hypersensitive sites or other chromatin marks have generally utilized existing datasets from normal tissues (53,54). These data have implicated T cells (CD4 > CD8), B cells, and monocytes in SLE. One study utilizing neutrophil ChIP-seq data found enrichment of genes with potential ETS TF binding near SLE GWAS SNPs (27). To our knowledge, there are no previous studies that have used disease-related changes in RNA and H3K4me3 to analyze GWAS variants. This study offers a new perspective on the function of the GWAS variants using our unique lens of SLE-related changes to the transcriptome and epigenome. This paves the way for improved approaches to functionality.
Our study had several limitations including a limited sample size, a result of pairing RNA, and ChIP-seq in three cell types. We also chose to focus on a single histone modification. We had previously examined H4 acetylation and H3K4me3 (18,20,21). This study focused on H3K4me3 due to the greater understanding of the variables impacting this specific modification. Studies of other modifications, other cell types, and validating patient cohorts will improve our overall understanding of the epigenome of SLE. Our GWAS analysis focused on the region immediately surrounding the variant. SNP linkage disequilibrium blocks may be much larger and therefore the analysis may have under-estimated the effect. Nevertheless, our study provides unique new data on the pathogenesis of SLE using an innovative approach. This study is the first to cross-reference GWAS data with chromatin marks directly impacted by disease. This study presents the most robust analysis of gene expression effects in purified cell types to begin to dissect the concordant and discordant effects in disease. A pivotal and surprising finding is that RNAs cross-referenced with GWAS SNPs had markedly decreased expression in SLE. This study also identified potential drug targets in ETS family members and ATF3. ETS family TFs have been implicated in SLE (45,55,56) and represent a credible target for drug development.
In summary, this study provides key insights into molecular pathways relevant in SLE. Slightly more than half of the gene expression changes in SLE patients were concordant. This approach could facilitate efforts at drug discovery where pathologic cell behavior could be more effectively targeted by focusing on concordant pathways. Nevertheless, this study also provides a sobering view of the complexity of this disease. Nearly, half of the gene expression changes were cell specific and identifying the key mediators in distinct cells may prove to be complex. As a first analysis of combined gene expression and histone modifications in SLE across different cell types, this study provides both hope and reservation. This innovative approach to focus GWAS data and to define conserved pathways offers strategies for drug development and analysis of variants.

eThics sTaTeMenT
The female healthy donor blood samples were obtained from the Center for AIDS Research, under a protocol approved by reFerences the University of Pennsylvania Perelman School of Medicine Institutional Review Board (IRB) in accordance with their guidelines. All subjects gave informed consent in accordance with the Declaration of Helsinki. Samples were handled identically as the shipped SLE samples. The female SLE patient blood samples were from the Johns Hopkins Lupus Cohort (23,24) under a separate IRB protocol in accordance with the recommendations and guidelines of the Johns Hopkins IRB Committee. All patients gave informed consent in accordance with the Declaration of Helsinki.
aUThOr cOnTribUTiOns LHS, LS, and KM performed experiments and did data analysis. ZZ, KS, and MP did data analysis and wrote the manuscript.

acKnOWleDgMenTs
The authors would like to thank the patients and their families for participation. Nursing support and staff assistance and enthusiasm are gratefully acknowledged.

FUnDing
This work was supported by the Wallace Chair of Pediatrics and NIH grants ES017627 (KS), AR43727/AR69572 (MP).