Comprehensive epigenomic profiling reveals the extent of disease-specific chromatin states and informs target discovery in ankylosing spondylitis

Summary Ankylosing spondylitis (AS) is a common, highly heritable inflammatory arthritis characterized by enthesitis of the spine and sacroiliac joints. Genome-wide association studies (GWASs) have revealed more than 100 genetic associations whose functional effects remain largely unresolved. Here, we present a comprehensive transcriptomic and epigenomic map of disease-relevant blood immune cell subsets from AS patients and healthy controls. We find that, while CD14+ monocytes and CD4+ and CD8+ T cells show disease-specific differences at the RNA level, epigenomic differences are only apparent upon multi-omics integration. The latter reveals enrichment at disease-associated loci in monocytes. We link putative functional SNPs to genes using high-resolution Capture-C at 10 loci, including PTGER4 and ETS1, and show how disease-specific functional genomic data can be integrated with GWASs to enhance therapeutic target discovery. This study combines epigenetic and transcriptional analysis with GWASs to identify disease-relevant cell types and gene regulation of likely pathogenic relevance and prioritize drug targets.


In brief
Brown et al. performed comprehensive transcriptomic and epigenomic mapping in ankylosing spondylitis. Comparison of blood immune cells from patients and healthy controls revealed diseasespecific chromatin patterns with potential pathogenic effects via mechanisms such as chromosome looping. The authors used these results to prioritize genetic variants, genes, and pathways for drug discovery.

INTRODUCTION
Ankylosing spondylitis (AS) is a common inflammatory arthritis characterized by inflammation of the sacroiliac joints and spinal entheses, which causes extensive new bone formation and vertebral fusion, resulting in pain, loss of movement, and disability. 1,2 Combined with other systemic manifestations of the disease, such as inflammation of the gut, skin, and eyes, this leads to significant morbidity and disease burden. 3,4 Twin and other family studies indicate that AS is highly heritable (ls $50) with broad-sense heritability greater than 90%. 5,6 This involves strong association with the major histocompatibility complex (MHC) allele HLA-B27 7,8 and more than 100 other loci identified through genome-wide association studies (GWASs). 9-13 Several of these associations implicate genes involved in interleukin-23 (IL-23)-driven inflammation and Th17 responses; these include IL23R (encoding the IL-23 receptor), IL6R (IL-6 receptor), TYK2 (tyrosine kinase 2 receptor), and IL27R. 9 In a few cases, such as ERAP1 (endoplasmic reticulum aminopeptidase 1) and IL23R, functional non-synonymous single-nucleotide Article ll OPEN ACCESS polymorphisms (SNPs) have been described. 9,13,14 However, most AS associations involve non-coding SNPs, which may be regulatory in nature and act in a cell-specific manner to modulate a variety of epigenetic, transcriptional, and post-transcriptional mechanisms. [15][16][17][18] Recently we demonstrated how AS-associated SNPs at RUNX3 modulate the binding of transcription factors (TFs) and regulatory complexes in T cells and monocytes, 15-18 but for other associated loci, the causal genes and pathways remain largely unresolved. 19 The expression and coordination of regulatory mechanisms for genes involved in the disease pathophysiology are likely to be cell type specific.
Here we focus on CD4 + and CD8 + T cells and monocytes from patients with active AS and healthy controls (HC) because these cell types have been implicated previously in the pathogenesis of AS, 20-26 and previous studies have largely sampled whole peripheral blood mononuclear cells (PBMCs). [27][28][29][30] The outlook for patients with more severe forms of AS has been greatly improved in recent years by the introduction of new biologic treatments inhibiting the inflammatory cytokines tumor necrosis factor alpha (TNF-a) and IL-17A. Nevertheless, fewer than half are likely to achieve sustained remission even with these targeted therapies, 31 highlighting the need for patient stratification of potential responders and new therapeutic targets in AS. Human genetic evidence supporting the identification of therapeutic targets strongly increases the likelihood of success in late-stage clinical trials. 32 We and others have shown that GWASs, in combination with functional genomic evidence and knowledge of network connectivity, can be used to prioritize target genes and pathways through, for example, the priority index (Pi) algorithm. 33,34 Logistical and technical challenges have limited the number of studies generating omics data from patient samples to date, and it remains unresolved how best to maximize the value of such data through integration across assay modalities, including genetics. 35 Here, we present a comprehensive map of the epigenomic landscape of AS defining the global transcriptome, chromatin accessibility, and enhancer-and promoter-associated histone modifications in disease-relevant subsets of immune cells from patients and HCs. We identify global changes in chromatin architecture in the AS disease state in monocytes and characterize specific GWAS loci to identify interactions between lead SNPs in enhancers and cognate genes, including prostaglandin E receptor 4 (PTGER4) and ETS proto-oncogene 1 (ETS1). Furthermore, we show how functional genomic evidence can be integrated with GWAS data through Pi to identify candidate therapeutic targets for future study.

Design
Patient and control cohorts and experimental overview To generate a comprehensive functional genomic and epigenomic atlas of the immune response in peripheral blood of AS patients, we recruited 20 adult patients with active disease who were naive to biologic therapy and fulfilled diagnostic criteria for AS and 35 HCs recruited locally or from the Oxford Biobank (Table S1; Table S2; STAR Methods). Cell populations of interest (CD4 + and CD8 + T cells and CD14 + monocytes) were freshly isolated from peripheral blood by positive selection using immunomagnetic cell separation with more than 98% purity ( Figure S1; STAR Methods). Each cell type was then processed immediately for total RNA sequencing (RNA-seq), chromatin accessibility (assay for transposase-accessible chromatin with next-generation sequencing [ATAC-seq]), informative histone modifications for promoter (histone H3 lysine 4 trimethylation [H3K4me3]) and enhancer (histone H3 lysine 4 trimethylation [H3K27ac]) activity (ChIPmentation), and high-resolution chromosomal conformation capture (Capture-C) (Figure 1A; Table S2; STAR Methods).

Differential gene expression in active AS is cell type specific
We first investigated the nature of differential gene expression between AS patients and HCs for specific immune cell types. We focused on three major immune cell types previously implicated in AS (CD4 + and CD8 + T cells and CD14 + monocytes). 20-26 Analysis of gene expression by RNA-seq comparing the three cell types for each individual showed that gene expression segregated by cell type more strongly than by disease state (Figures 1B and S1B). However, for each cell type, we found hundreds of differentially expressed genes between AS patients and HC (CD4 + T cells, 122 genes; CD8 + T cells, 299 genes; monocytes, 300 genes; padj (adjusted p value) < 0.05, fold change [FC] > 1.5) (STAR Methods; Figures 1C and S1C; Table S3). The majority of differentially expressed genes were cell type specific ( Figure 1D), and where differentially expressed genes involved more than one cell type, the direction of effect was the same in the majority of cases. For example, CD83 (encoding (B) PCA of RNA-seq data in CD4 + T cells, CD8 + T cells, and monocytes from AS patients and HCs. (C) Volcano plot showing differentially expressed genes calculated using DEseq2 (padj < 0.05, FC > 1.5) between AS patients and HCs in CD4 + T cells, CD8 + T cells, and monocytes. Genes in AS-associated GWAS regions are purple. Red genes are upregulated and blue downregulated in AS patients. CD4 + T cells had 122 differentially expressed genes, CD8 + T cells 299 genes, and monocytes 300 genes. (D) Cell type specificity of differentially expressed genes; numbers of differentially expressed genes are given. (E) Examples of differential gene expression at CD83 and TNFSF14; **padj < 0.01, ***padj < 10 À7 (from DEseq2). (F) Enriched pathways in the Reactome database (FDR < 0.01 from XGR output) from significant differentially expressed genes in each cell type. Dot size represents percentages of genes represented in that pathway, and colors represent cell types. (G) CXC subfamily of the Kyoto Encyclopedia of Genes and Genomes (KEGG) ''cytokine-cytokine reception interaction'' pathway colored by gene expression log 2 FC. Significantly differentially expressed genes are marked by asterisks. The p value of CXC family subset over-representation is shown below for each cell type, calculated by chi-squared test with Yates' correction. See also Figure S1 and Table S3. CD83, a cell-surface glycoprotein involved in regulation of antigen presentation) has significantly higher expression only in AS patient monocytes, TNFSF14 (encoding TNF superfamily member 14) has significantly lower expression in AS patients in all three cell types ( Figure 1E), and SCAMP5 is an example of a gene upregulated in CD8 + T cells and downregulated in CD4 + T cells. We defined disease-enriched pathways from differentially expressed genes using eXploring Genomic Relations (XGR) 36 (STAR Methods). All cell types showed enrichment of immune-related pathways and G protein-coupled receptor (GPCR) signaling pathways in AS ( Figure 1F). We identified significant upregulation of the CXC subfamily of chemokine receptors in CD4 + and CD8 + T cells ( Figures 1G and S1D), which links to the important role of IL-17-producing cells in AS pathogenesis. We investigated cell subset composition by deconvolution of RNA-seq data using CIBERSORTx 37 and found no difference in abundance of the three cell types or major cell subsets within these when comparing AS patients and HCs ( Figure S1E).
Cell-type-specific epigenomic marks show limited differences between AS patients and controls We then resolved genomic regulatory features in CD4 + T cells, CD8 + T cells, and monocytes from AS patients. To do this, we assayed open chromatin with ATAC-seq and enrichment of H3K4me3 and H3K27ac histone modifications using ChIPmentation (ChIPm) and identified non-coding enhancer RNAs (eRNAs) within ATAC peaks outside of coding genes (STAR Methods; Figures S2A and S2B). Each omics modality showed cell type specificity on principal-component analysis (PCA) that outweighed the effect of disease state (Figure 2A). Only a very small number of significant differential signals were observed between AS patients and HCs (Figure 2B;  Tables S4-S6), and disease state could not be clearly distinguished on PCA ( Figure S2C). The transcription start site (TSS) score for ATAC and chromatin immunoprecipitation (ChIP) correlated with expression levels of their corresponding gene (Figure S2D). We assigned differential ATAC, ChIPm, and eRNA signals to genes by proximity or overlap with promoter capture Hi-C (PCHi-C) looping interactions identified in relevant cell types 38 (STAR Methods). Pathway enrichment analysis of genes linked to the top 200 differential regions between AS patients and HCs for each modality implicated immunological pathways and GPCR signaling along with transcriptional pathways and NOTCH signaling across cell types and modalities, consistent with our analysis of differentially expressed genes ( Figure 2B).
Disease-specific regulatory chromatin states are found in monocytes from patients with AS We further investigated whether there were disease-specific differences in the chromatin landscape by maximizing the informativeness of different sources of omics information using ChromHMM 39 ( Figure 3A Figure 3B) in four functional categories: promoter (states 1-4), transcribed (states 5-7), enhancer (states 8-11), and quiescent (states 12-14). Multiple correspondence analysis (MCA) on global ChromHMM data revealed major differences between cell types, albeit with considerable variation between individuals ( Figure S3A). Within each cell type, we identified thousands of 200-bp segments that were assigned to different chromatin states in AS patients and HCs (Table S7). We found that differential segments in two promoter states (1 and 3), three enhancer states (9, 10, and 11), and two quiescent states (12 and 14) were significantly over-represented in AS monocytes (Figure 3C), as determined by permutation analysis ( Figure S3B). Consistent with this, MCA comparison of AS patients and HCs showed separation on dimension (dim) 1 in monocytes but not CD4 + and CD8 + T cells ( Figure 3D). Within monocytes, all states except state 14 (quiescent) showed separation of AS patients and HCs on dim 1 or dim 2 ( Figure S3C). In contrast, only three states showed separation in CD8 + T cells (states 7, 8, and 11), and none in CD4 + T cells. Comparative state transitions between patients and HCs revealed a complex pattern of differential states, including enrichment of enhancer state 10 (EnhA) in AS patients corresponding to promoter state 1 (TssA) in HCs ( Figure S3D). Overall, this analysis demonstrates significant changes in the epigenomic landscape of active AS disease, notably in monocytes.
We next sought to determine which molecular pathways might be altered by these global changes in monocyte chromatin architecture. We assigned the differential ChromHMM fragments to genes based on proximity or PCHi-C looping events 38 (STAR Methods; Table S7) and performed pathway enrichment analysis ( Figure 3E). The highest number of enriched pathways contained genes linked to promoter state, followed by those linked to enhancer state. Only two pathways were enriched in genes linked to transcribed states and none with quiescent states. Six of 15 enriched pathways from promoter and enhancer states related to NOTCH signaling ( Figure 3E), further implicating this pathway in monocytes in AS. To aid interpretation of these findings, we further defined by flow cytometry which monocyte subpopulations were represented in our sorted monocyte population, showing that these are 80% CD14 + CD16 + classical monocytes with the remainder intermediate and non-classical monocytes ( Figure S3E).
To illustrate our results, we show two examples where differential ChromHMM segments correlate with alterations in gene expression ( Figure 3F). ELOVL3 (encoding elongation of verylong-chain fatty acids) is involved in fatty acid metabolism and downregulated in psoriasis, 40 a common extra-articular manifestation of AS. 41 The ELOVL3 promoter is marked by state 1 (TssA, promoter-like) in AS patients and HCs, but this mark spans a shorter genomic interval in AS patients ( Figure 3F). Consistent with this, expression of ELOVL3 is lower in AS patients (padj = 0.0003). TLE3 (encoding transducin-like enhancer family member 3) is a transcriptional co-repressor involved in the NOTCH signaling pathway. We identified up-and downstream enhancers that form looping interactions with the TLE3 promoter ( Figure 3F). The enhancers are defined by ChromHMM state 10 (enhancer) containing punctate regions of state 1 (promoter) that correspond with non-coding eRNA transcription. In two segments these promoter elements are narrower in AS patients compared with HC, which may contribute to the observed reduction in TLE3 expression in AS patients (padj = 0.045).

Regulatory chromatin signatures are enriched at AS GWAS regions
Having generated comprehensive epigenomic maps in three cell types from AS patients and HCs, we next addressed whether this could inform the functional basis of observed genetic associations in AS from GWASs, specifically seeking evidence to implicate/delineate the gene(s) responsible for the genetic association. We identified 35 differentially expressed genes located within GWAS regions (<500 kb from the lead SNP), of which 31 were cell type specific ( Figure 1C; STAR Methods). We found that regulatory ChromHMM states are over-represented in or near AS GWAS regions, and this is independent of the HLA-B27 association ( Figures 4A and S4). Enhancer state 10 is enriched near GWAS regions in all three cell types, while in monocytes, states 2, 6, 10, 11, 12, 13, and 14 are all enriched in GWAS regions. We also found significant enrichment of ATAC, H3K4me3, and H3K27ac signals and eRNAs at GWAS loci in all cell types ( Figure 4B).
The chr3p21.31 locus illustrates the intersection between an AS genetic association (rs1001007 10 ) and disease-specific chromatin state in monocytes ( Figure 4C). This intergenic SNP is (F) Visualization of differential ChromHMM regions (vertical arrows) at (i) the promoter of ELOVL3 and (ii) enhancer of TLE3, with (iii) magnified region of the TLE3 differential enhancer. ChromHMM states in AS patients and HCs are colored as in (B). PCHiC looping interactions 38 are shown, with loops intersecting differential ChromHMM segments in red. See also Figure S3 and Table S8. Article ll OPEN ACCESS located between CCR5 and CCRL2, but this locus has not been fine mapped to establish a functional variant and modulated gene in AS. Consequently, we interrogated the whole monocyte-specific topologically associating domain (TAD) and discovered several regions with disease-specific chromatin signatures. We found evidence of an intergenic region between XCR1 and CCR1 that has an enhancer-like profile (ChromHMM state 10, EnhA) in AS patients with a promoter-like profile (ChromHMM state 3, FlnkU) in HCs, consistent with observed differences in eRNA and H3K27ac between cases and controls.
Analysis of chromatin interaction data shows that this region is involved in DNA looping to CCR1, CCR2, CCR5, and CCRL2 but not with CCR3. Looping events correlate with expression of the cognate genes in monocytes. These findings indicate complex regulation of the CCR gene cluster in monocytes with epigenetic changes specific to AS.

Chromosome looping interactions from AS-associated SNPs in enhancers implicate disease-relevant genes
To further substantiate the relationship between disease-associated SNPs and likely modulated genes, we mapped chromosome looping events at GWAS loci in patient samples at high resolution (compared with PCHi-C). We performed Capture-C on CD4 + T cells, CD8 + T cells, and monocytes from AS patients and controls (STAR Methods). Baits were designed at promoters of genes with known immune roles within AS GWAS regions. We found that 18 of 44 promoter viewpoints assayed demonstrated chromatin interactions in at least one cell type, although no differences were found between AS patients and HCs. Follow-up Capture-C experiments were performed with baits at AS-associated SNPs to demonstrate reciprocal interactions between SNPs and promoters. Overall, nine reciprocal interactions were found between promoters and regions containing GWAS SNPs, nine of which were marked as enhancers (state 10 or 11) in our ChromHMM analysis in the same cell type (Table 1; Figures 5 and S5-S14). These loci also contained differential ATAC, H3K4me3, H3K27ac, or eRNA peaks, and five overlapped expression quantitative trait loci (eQTLs) in the same cell types 43 ( Figure 5A).
Two loci illustrate SNP-promoter interactions. ETS1 encodes ETS proto-oncogene 1, a TF with numerous roles in immune cells, including regulation of cytokine and chemokine gene expression. 46 Our ChromHMM analysis showed enhancer regions flanking the ETS1 gene, and Capture-C analysis showed looping events between the enhancer overlapping the AS-associated lead SNP rs7933433 10 and ETS1 promoter in CD4 + and CD8 + T cells but not monocytes ( Figures 5B and S11). PTGER4 encodes prostaglandin receptor E4, a GPCR whose expression is associ-ated with disease severity in AS. 47 We observed an interaction between the promoter of PTGER4 and the enhancer overlapping rs1992661 (AS GWAS lead SNP 10 ) in all cell types ( Figures 5C  and S9). We detected an interaction between an enhancer encompassing the known functional SNP rs9283753 48 and the PTGER4 promoter specifically in monocytes and only with the promoter bait ( Figures 5C and S9). There was no detectable interaction between the lead SNP rs12186979 and PTGER4 gene in any of the cell types. These results support a model where associated SNPs lying in an enhancer region interact with a distal gene promoter via a chromatin looping event. This, combined with additional local disease-context-specific epigenomic modifications, may lead to alterations in cognate gene expression.
We were interested to explore whether the genotypes of individual SNPs were associated with alterations in chromatin structure. The study was underpowered to perform such an analysis genome wide, but we were able to analyze the effect of individual lead AS-associated SNPs. In doing so, we identified an association between rs4672505 and ATAC-seq peak chr2:62559366-62561099, where the A risk allele correlates with a reduced ATAC-seq signal in CD8 + T cells ( Figure 5D). Analysis of nextgeneration sequencing (NGS) reads from the 16 heterozygous individuals with at least 5 mapping reads showed that 99.4% of ATAC peak reads encoded the G allele, strongly suggesting that the risk A allele prevents chromatin opening. This finding was not driven by mapping bias because A was the reference allele in the hg19 build used in this analysis. This SNP lies at chr2p15 and overlaps a putative enhancer that forms a looping interaction with B3GNT2 identified by publicly available PCHi-C ( Figure 5D). rs4672505 is associated with AS, Crohn disease, and psoriasis 10 and is an eQTL for B3GNT2. 45 Use of disease-specific functional genomic datasets enhances drug target discovery in AS The final aim of our study was to prioritize new therapeutic targets in AS. We previously developed Pi, a genetics-led approach that annotates GWASs with functional genomic data to prioritize therapeutic targets across a range of immune-mediated diseases. 33, 34 We modified the underlying algorithm of Pi to include our new AS-specific functional genomic datasets ( Figure 6A; STAR Methods), the algorithm previously having been limited to non-disease-context functional genomics data, and assessed whether this inclusion increased the power to identify potential therapeutic targets. In the original algorithm, we used diseasespecific genetic associations to define seed (core) genes, including (1) nearby genes (nGene) using genomic proximity and organization, (2) expression-associated genes (eGene) Figure 4. Differential chromatin regions are enriched at GWAS loci (A) Enrichment of differential ChromHMM differential regions at AS-associated GWAS loci 10 with association p value thresholds as indicated in CD4 + T cells, CD8 + T cells, and monocytes. ***p < 0.001, **p < 0.01, *p < 0.05. OR, odds ratio. (B) Enrichment of ATAC, H3K4me3, H3K27ac, and eRNA peaks at AS-associated GWAS loci 9 in CD4 + T cells, CD8 + T cells, and monocytes; GWAS association p value thresholds are indicated. Error bars represent 95% confidence interval of the OR (from GARFIELD 42 ). ***p < 0.001, **p < 0.01, *p < 0.05. (C) Visualization of multiple epigenomic datasets at chr3p21 (top, chr3:45890000-46600000; bottom, chr3:46090000-46200000). PCHiC looping events 38 that intersect regions of differential chromatin are indicated in red. Selected gene transcripts from the Ensembl database are shown. ChromHMM is shown for 4 AS patients and 4 HCs with colors as in Figure 3B. The lead AS-associated SNP rs1001007 10 is shown. Representative ATAC, H3K4me3, and H3K27ac tracks (RPKM) and total RNA tracks (log2 count) are shown for AS patients (orange) and HCs (blue), with called peaks marked by gray boxes and differential peaks shown in red (upregulated in AS) or blue (downregulated in AS). See also Figure S4.  using PCHi-C datasets. Here, we added five types of AS-specific functional genomic predictors using data from this study (denoted RNA, eRNA, H3K27ac, H3K4me3, and ATAC). The AS-specific expression predictor (RNA) was generated based on differential gene expression, and AS-specific epigenomic predictors (ATAC, H3K4me3, H3K27ac, and eRNA) were prepared on the basis of differential peaks linked to genes as above.
We benchmarked the performance of the Pi algorithm with and without the AS-specific functional genomic predictors to prioritize currently approved drug targets for AS versus simulated negative targets (STAR Methods). This showed that inclusion of disease-specific data improved the predictive power (area under the curve [AUC] = 0.869) compared with the original prediction (AUC = 0.822) and the state-of-the-art approach (Open Targets, 49 including text mining [AUC = 0.761] and genetic associations [AUC = 0.582] from the Open Targets Genetics Portal 50 ) ( Figure 6B). As expected, combined use of predictors performed much better than each predictor alone ( Figures 6B and S15A). Among the top 1% of prioritized genes (of >17,000 ranked genes) were the known AS drug targets IL23R, JAK2, and TYK2 ( Figure S15B). Pathway enrichment analysis of the top 1% of genes identified pathogenic AS pathways, such as Janus kinase-signal transducer and activator of transcription (JAK-STAT) signaling (false discovery rate [FDR] = 3.3 3 10 À32 ), TNF signaling (FDR = 1 3 10 À16 ) and Th17 cell differentiation (FDR = 3.1 3 10 À24 ) ( Figure 6C). The significantly enriched pathways also included cytokine-cytokine receptor interaction, Tolllike receptor signaling, ErbB signaling, chemokine signaling, and T cell receptor (TCR) signaling, which were consistent with our findings from the RNA-seq data ( Figure S2C). Pathway crosstalk analysis for potential therapeutic intervention identified a network of 53 interconnecting genes ( Figure 6D), including genes that are already therapeutic targets in AS and other autoimmune conditions, such as JAK1 and TYK2, alongside yet unexplored genes. The list of genes within this network ( Figure 6E) are highly prioritized as candidates for therapeutic intervention, providing an input for the drug development pipelines and further study. 51

Summary of findings
We performed transcriptomic and epigenomic profiling in specific primary immune cell populations isolated from carefully phenotyped patients with active AS. This revealed disease-specific changes in the chromatin landscape and differential regulatory signatures significantly enriched near AS-associated loci. Pathway analysis indicated the importance of NOTCH and chemokine receptor signaling in AS in addition to known diseaseassociated pathways. Capture-C identified physical interactions between associated SNPs lying in enhancer regions and nearby gene promoters. Taken together, these results provide a map of the epigenetic landscape in AS and evidence of the mechanisms by which genetic associations can alter immune cell function in AS. Earlier studies have compared gene expression between AS patients and HCs in PBMCs. [27][28][29][30] By looking at individual cell types, in this study we demonstrated a key role of monocytes in the pathogenesis of AS. This work highlights the importance of looking for functional effects of SNPs in the appropriate cellular and disease contexts in AS and other immune-mediated conditions.
Identification of regulatory disease-associated SNPs and cognate genes When analyzing the effects of a genetic variant with putative enhancer-modifying activity, a key component is to understand which gene(s) are regulated by that enhancer. Here, we described 10 enhancers that overlap AS-associated SNPs and that interact with a gene promoter within the same TAD (Table 1; Figures 5 and S5-S14). These enhancers share common features such as eRNA expression, activating histone marks, and open chromatin marks.
We discovered enhancer-gene interactions with two members of the ETS proto-oncogene (ETS) family of TFs, ETS1 and ETS2 (Figures 5, S11, and S13). At the ETS1 locus, we observed chromatin looping interactions with dual enhancers, up-and downstream of the gene, specifically in T cells, of which the downstream enhancer overlaps an AS-associated lead SNP (rs7933433 10 ). A monocyte-specific interaction was observed between an enhancer and the ETS2 promoter, which overlaps the 99% credible set comprising 5 AS-associated SNPs 52 that are also eQTLs in monocytes. 43 ETS1 and ETS2 enhancers exhibit non-coding eRNA transcription. ETS1 and ETS2 are TFs expressed across various immune cell types and have roles including regulation of T cell subset differentiation. 53 ETS1 regulates the expression of IL-7R 54 (encoding the IL-7Ra subunit), which is also associated with AS, 20 and RUNX3, which encodes Runt-related TF 3, a TF involved in T cell function with known AS-associated functional variants. 15-18 ETS1, ETS2, and IL7R are in the Pi network output ( Figure 6D), indicating that they form part of an important functional pathway with strong possibility for therapeutic intervention. Gene-SNP interactions were also observed at the PTGER4 locus ( Figures 5 and S9). PTGER4 is widely expressed throughout the immune system and is involved in the IL-23 and TNF-a pathways. 55 It is also expressed in osteoclasts and could potentially have a role in new bone formation in AS. 56 GWASs have found two independent associations at the PTGER4 locus. 9 Tewhey et al. 44 described a functional SNP (rs9283753) that lies in an enhancer and alters PTGER4 expression in lymphoblastoid cell lines. Our data show a chromatin looping interaction between AS-associated SNP rs1992661 and the PTGER4 promoter in all three cell types and a T cell-specific interaction between rs9283753 and the PTGER4 promoter. Our ChromHMM analysis shows that these SNPs overlap enhancer marks in the same cell types. Taken together, these data strongly support a functional role of rs1992661 and rs9283753 in regulation of PTGER4.
We found allele-specific differences in the ATAC-seq signal for an AS GWAS-associated SNP, rs4672505. This variant has been associated previously with differential abundance of B3GNT2, a poly-N-acetyllactosamine synthase in whole blood. 45 B3GNT2 is upregulated in T cells on activation, and a recent CRISPR screen showed evidence that this enzyme is important in modulating T cell activation in the setting of cancer. 57 Our findings suggest that reduced expression of B3GNT2 in individuals with the AS risk allele is likely caused by reduced chromatin openness at this locus. This may impact higher-order chromatin structures, such as the looping event identified between a distal enhancer and B3GNT2 ( Figure 5). The JASPAR 2022 58 database of TF binding profiles predicts that STAT1 binds to this site, a TF that plays an important role in transcriptional activation in the immune system. 59 Further work, such as genomic editing and functional assays, will be needed to identify the function of these or other unknown TFs at this locus.
Cytokine and NOTCH signaling pathways GPCR-related and cytokine signaling pathways were consistently enriched across modalities and cell types. In particular, expression of CXC cytokine subfamily genes was upregulated in CD4 + and CD8 + T cells, consistent with the inflammatory environment of AS. Chemokine levels have been shown previously to be increased at the gene and protein level in AS 60 and psoriatic arthritis, a related spondyloarthropathy. 61,62 This may contribute to the differentiation of pro-inflammatory Th1 and Th17 cells, which are expanded in AS, or trafficking of leukocytes and osteoclast precursors to the sites of inflammation at the joint. 63 NOTCH signaling has a wide range of functions in the innate and adaptive immune systems. 64 We found that NOTCH signaling was linked to changes in chromatin signatures in monocytes. NOTCH signaling has been shown to be important for monocyte-macrophage differentiation, with increased NOTCH signaling favoring inflammatory M1 macrophage development in atherosclerosis, systemic lupus erythematosus, and cancer. 65 NOTCH signaling has already been implicated in inflammatory states 66 and rheumatoid arthritis 67 and may be important for regulating monocyte differentiation to osteoclasts, 68 thus influencing the pathogenic ossification that is a key feature of latestage AS. 69 Wang et al. 70 showed that NOTCH1 expression is reduced in AS patients who have had anti-TNF biologic therapy. Further study of this pathway in AS is warranted and could lead to repurposing of existing NOTCH pathway inhibitors in AS. 71

Identification of novel therapeutic targets
We showed how disease-context-specific functional genomic data could be used to identify novel therapeutic targets in AS. We modified Pi, a previously published algorithm designed to identify the network of therapeutic targets in autoimmune disease from GWAS. The multi-omics approach used here supported the importance of known candidate pathways, such as Th17/IL-23 and TNF, and identifies new pathways and potential drug targets, including PTGER4, ErbB, phosphatidylinositol 3-kinase (PI3K)/AKT, NOTCH, and GPCR ( Figure 6). Existing inhibitors of these pathways could be repurposed in AS, such as PI3K inhibitors used in lymphoma treatment. 72 Our results show exciting promise for development of new therapeutics in AS. The exact roles these TFs and other network genes play in AS remain to be elucidated and will be investigated in the future using, for example, genomic editing and functional assays.

Conclusions
We demonstrated that the epigenomic landscape of immune cells is altered in AS. We used these results, together with evidence of chromosomal interactions, to inform the interpretation of GWASs for AS in terms of likely functional variants and modulated genes and prioritize potential drug targets and networks. This is important because existing therapies are only effective in a subset of AS patients and ultimately do not cure the disease.

Limitations of the study
This study used disease-relevant immune cell subsets isolated from PBMCs of AS patients and HCs. One limitation of our study was the use of bulk CD4 + T cell, CD8 + T cell, and CD14 + monocyte populations rather than investigating smaller subsets of these cell types. This was necessary because of the large cell numbers required to perform a panel of multi-omics experiments on the same sample, in particular mapping of chromosomal interactions. We showed that major cell subsets were present at expected frequencies and did not differ between AS patients and HCs (Figure S1), although disease-specific cell type frequencies have been reported for minor cell subtypes. 20,25,26,73 Future work sampling cells from sites of inflammation, such as sacroiliac joints, and utilizing single-cell-based methods will further unravel the cell types and pathogenic mechanisms of AS. 74 For individual modalities, signals from differential analysis were modest, especially in CD4 + T cells, which could be due to the heterogeneity of this cell type. Context specificity of regulatory regions is key, so in addition to direct ex vivo analysis of cells from patients with the active disease state, analysis of such cells subjected to immune challenges in vitro (such as lipopolysaccharide stimulation of monocytes, anti-CD3/28 stimulation of T cells) may identify additional functional SNPs specifically in those activation states. The study was underpowered to perform expression and chromatin quantitative trait mapping on a genome-wide scale, and this is an important area for future work in a disease context. We were unable to study the effect of other covariates, including sex and drug regimens (although all patients were biologic therapy naive), because of the small sample size. We have shown previously that eRNAs have a role in innate immune activation. 75 The observed widespread bidirectional eRNA expression at genomic enhancers can be further investigated using more sensitive methods (such as global run-on sequencing, small capped RNA sequencing, and precision run-on sequencing 76 ) to detect more subtle alterations in eRNA expression and specifically identify the role of SNPs therein. Findings from ChromHMM were limited by small sample size. We focused on the presence of activating chromatin modifications (chromatin accessibility, promoter-and enhancerassociated histone modifications), so future studies could investigate the role of repressive histone marks such as H3K9me3 and H3K27me3. Future studies will be required to further characterize the genes and pathways highlighted by this study to assess their effect at the protein level and on cellular phenotype and function; for example, through genome editing and smallmolecule inhibitors.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: