Evidence for a second ankylosing spondylitis-associated RUNX3 regulatory polymorphism

Objectives To explore the functions of RUNX3 single nucleotide polymorphisms (SNPs) associated with ankylosing spondylitis (AS). Methods Individual SNP associations were evaluated in 4230 UK cases. Their effects on transcription factor (TF) binding, transcription regulation, chromatin modifications, gene expression and gene interactions were tested by database interrogation, luciferase reporter assays, electrophoretic mobility gel shifts, chromatin immunoprecipitation and real-time PCR. Results We confirmed the independent association of AS with rs4265380, which was robust (P=4.7×10−6) to conditioning on another nearby AS-associated RUNX3 SNP (rs4648889). A RUNX3 haplotype incorporating both SNPs was strongly associated with AS (OR 6.2; 95% CI 3.1 to 13.2, P=1.4×10−8). In a large UK cohort, rs4265380 is associated with leucocyte counts (including monocytes). RUNX3 expression is lower in AS peripheral blood mononuclear cells than healthy controls (P<0.002), independent of rs4265380 genotype. Enhancer function for this RUNX3 region was suggested by increased luciferase activity (approximately tenfold; P=0.005) for reporter constructs containing rs4265380. In monocytes, there was differential allelic binding of nuclear protein extracts to a 50 bp DNA probe containing rs4265380 that was strongly augmented by lipopolysaccharide activation. TF binding also included the histone modifier p300. There was enrichment for histone modifications associated with active enhancer elements (H3K27Ac and H3K79Me2) that may be allele dependent. Hi-C database interrogation showed chromosome interactions of RUNX3 bait with the nearby RP4-799D16.1 lincRNA. Conclusions The association of AS with this RUNX3 regulatory region involves at least two SNPs apparently operating in different cell types. Monocytes may be potential therapeutic targets in AS.


AbstrAct
Objectives to explore the functions of RUNX3 single nucleotide polymorphisms (SnPs) associated with ankylosing spondylitis (aS). Methods individual SnP associations were evaluated in 4230 UK cases. their effects on transcription factor (tF) binding, transcription regulation, chromatin modifications, gene expression and gene interactions were tested by database interrogation, luciferase reporter assays, electrophoretic mobility gel shifts, chromatin immunoprecipitation and real-time Pcr. Results We confirmed the independent association of aS with rs4265380, which was robust (P=4.7×10 −6 ) to conditioning on another nearby aS-associated RUNX3 SnP (rs4648889). a RUNX3 haplotype incorporating both SnPs was strongly associated with aS (Or 6.2; 95% ci 3.1 to 13.2, P=1.4×10 −8 ). in a large UK cohort, rs4265380 is associated with leucocyte counts (including monocytes). RUNX3 expression is lower in aS peripheral blood mononuclear cells than healthy controls (P<0.002), independent of rs4265380 genotype. enhancer function for this RUNX3 region was suggested by increased luciferase activity (approximately tenfold; P=0.005) for reporter constructs containing rs4265380. in monocytes, there was differential allelic binding of nuclear protein extracts to a 50 bp Dna probe containing rs4265380 that was strongly augmented by lipopolysaccharide activation. tF binding also included the histone modifier p300. there was enrichment for histone modifications associated with active enhancer elements (H3K27ac and H3K79Me2) that may be allele dependent. Hi-c database interrogation showed chromosome interactions of rUnX3 bait with the nearby rP4-799D16.1 lincrna. Conclusions the association of aS with this RUNX3 regulatory region involves at least two SnPs apparently operating in different cell types. Monocytes may be potential therapeutic targets in aS.

InTROduCTIOn
Although genetic predisposition to ankylosing spondylitis (AS) is strongly associated with the human leukocyte antigen (HLA)-B27 immune response gene genome-wide association studies (GWAS) have clearly shown that it is polygenic. 1 2 Several robust associations have been identified with single nucleotide polymorphisms (SNPs) in coding sequences, including other genes in the major histocompatibility complex, ERAP1 and IL23R. [1][2][3][4][5][6] However, as in other complex polygenic disorders, AS-associated SNPs are more typically located in non-coding regions of the genome, such as the IL23R-IL12RB2 intergenic region, where they may influence gene expression. 7-10

Key messages
What is already known about this subject? ► From genome-wide association studies, it appears that more than 100 genes are involved in susceptibility to ankylosing spondylitis (aS). One of these is RUNX3, a transcription factor involved in the maturation of several immune cell types. ► Several single nucleotide polymorphisms (SnPs) upstream the promoter of RUNX3, which are strongly associated with aS lie in a putative genomic regulatory region. One of these SnPs, rs4648889, appears to be operative particularly in cD8+ t cells.
What does this study adds?
► We report potential functional effects from a second independently associated SnP in this RUNX3 region (rs4265380) that may be restricted to monocytes rather than t cells. the independent associations of aS with these two SnPS and the likelihood that they have functional effects in two different cell types suggest that both cD8+ t cells and monocytes could be therapeutic targets in aS. Understanding the functions of this RUNX3 could provide new insights into the pathogenesis of aS.
How might this impact on clinical practice? ► By investigating the aS-associated mechanisms controlling rUnX3 transcription and the regulatory effects of rUnX3 itself on downstream gene expression, we expect to identify pathways that include credible therapeutic targets for future drug development.

RMD Open
Convincing associations between AS and SNPs upstream of RUNX3 (Runt-related transcription factor 3) have previously been demonstrated by GWAS. 4 11 12 Subsequently, the International Genetics of AS (IGAS) Immunochip study demonstrated that 22 SNPs in a ~15 kb linkage disequilibrium (LD) block upstream of the RUNX3 promoter are strongly associated with AS (lead SNP rs6600247: P=1.3×10 −14 ). 1 We have previously shown that this LD block contains a putative regulatory region; further, we demonstrated that the risk allele at rs4648889 within this block reduces recruitment of the interferon regulatory factor 4 (IRF4) transcription factor (TF), leading to reduced RUNX3 expression in CD8+ T cells in an allele-dependent fashion. 11 RUNX3 not only plays a key role in the development of CD8+ T cells 13 14 but also has plausible roles in the development of many other immune cell types, including the differentiation, activity and function of monocytes. [15][16][17] Both CD8+ T cells and monocytes have been shown to have a role in AS pathogenesis. 18 19 Here, we present evidence for the complexity of both the molecular and cellular associations of RUNX3 with AS. Having previously established the strong association of AS with rs4648889 in this region (P=1.3×10 −14 ), 11 we now confirm the presence of a second independent association with a second SNP-rs4265380. We characterise some of the functional effects of this SNP and conclude that these two adjacent but independent signals probably exert their effects in different cell types.

MaTeRIals and MeTHOds Genotyping
Historical genotype data from the UK subset of cases included in the AS Immunochip study (n=4230) and ethnically matched healthy controls (n=9200) were used in the association studies. 1 Where necessary for cases used in the functional studies, DNA was extracted using the Qiagen AllPrep DNA/RNA Mini Kit (Qiagen Ltd, Manchester, UK) and genotyped for rs4265380 using TaqMan SNP assay (Life Technologies, Paisley, UK).

electrophoretic mobility gel shift assay (eMsa)
A 50-bp DNA probe including rs4265380 was incubated with nuclear extracts obtained either from primary CD14+ monocytes stimulated with LPS for 24 hours or a monocyte cell line from histiocytic lymphoma (U937) stimulated with phorbol-12-myristate-13-acetate (PMA). Nuclear extracts from monocytes and U937 cells were prepared using the Thermo Scientific NE-PER Nuclear and Cytoplasmic Extraction kit (Thermo Scientific, Waltham, USA). EMSA was performed using LightShift Chemiluminescent EMSA Kit (Thermo Scientific) with no modifications. Single-stranded biotinylated oligonucleotides (prepared with biotin 3′ end DNA labelling kit) were mixed and annealed at room temperature for 1 hour. Five micrograms of nuclear extract and 0.5 ng biotin labelled doublestranded oligonucleotides (50 bp fragment; Eurofins, Wolverhampton, UK) were used in each experiment. Unlabelled probes in 100-fold excess were used as competitor.
Quantitative real-time RT-PCR RNA was isolated with TRIzol (Invitrogen, Paisley, UK) and reverse transcribed with Superscript III (Invitrogen)

RMD Open
to synthesise cDNA as previously described. 11 The specific primers were: RUNX3 sense PCRs were performed in triplicate and the 2−ΔCt method was used to calculate the expression of each gene. All values were normalised to β-actin (ID Assay qHsaCED0036269, Bio-Rad Laboratories, UK). miRNA quantitative real-time RT-PCR miRNeasy Mini Kit (Qiagen Ltd) was used to isolate total RNA, according to the manufacturer's instructions. Ten nanograms of RNA was reverse transcribed and quantified by TaqMan real-time PCR using the TaqMan microRNA Assay (hsa-miR-4425, assay ID: 462954_mat, Applied Biosystems). RNU6 was used as endogenous control (U6 snRNA, assay ID: 001973). For detection, we used a vii7 PCR machine (Applied Biosystems).

Bioinformatics
We used publicly available data from the ENCODE (https:// genome. ucsc. edu/ ENCODE/) and Roadmap Epigenomics (http:// epigenomegateway. wustl. edu/ browser/ roadmap/) projects and previously published GWAS association data 1 to identify regulatory elements, histone modifications and TF binding sites relevant to AS in the RUNX3 region. 23 24 Potential chromatin interactions were evaluated using data generated by the International Human Epigenome Consortium (https://www. chicp. org/ chicp). 25 RUNX3 transcription in AS cases and controls was evaluated from historical data derived from RNA-seq in PBMCs from 72 AS cases and 62 healthy controls. 19 statistical analysis Association analysis was performed using the logistic regression function in PLINK, 26 V.1.07, accounting for population structure with 10 principal components as including H3K4Me1, H3K27Ac and H3K79Me2 for CD14+ monocytes are shown. ChIP-seq for p300 transcription factor is shown for the available lymphoblastoid cell line GM12878. 23 24 The position of rs4648889 SNP and the location of the 250-bp construct used in luciferase assay are also shown. Spondyloarthritis covariates in the regression analysis and conditioning on the RUNX3 SNP rs4648889. 11 Haplotypes were inferred by phasing using the SHAPEIT tool 27 and all SNPs genotyped in the locus (chr1: 23038732-27660734). Haplotype counts in cases and controls were compared using the χ 2 test and measured relative to the most common haplotype found in the control samples. One-way analysis of variance and two-tailed Student's t-test were used to determine statistical significance using the GraphPad Prism software (V.7.03) package.

ResulTs
Identification of rs4265380 as a putative candidate causal variant in monocytes The region upstream of RUNX3 contains potentially important regulatory sequences. However, the strongest association with AS in this region is with rs6600247 (P=4.23×10 −9 ) that does not overlap with obvious functionally active sites. We therefore investigated SNPs in strong LD with rs6600247 that lie in a region associated with open chromatin for evidence of regulatory activity; these include rs4648889 (unconditional analysis P=9.5×10 −9 , figure 1A) that is functionally active in CD8+ T cells.
Of the 22 RUNX3 SNPs previously associated with AS in the IGAS Immunochip study at P≤5×10 −7 , rs4265380 (P=4.58×10 −7 ) was the most strongly associated (with the 'C' allele) after conditioning on rs6600247 (lead IGAS SNP) or rs4648889 ( figure 1B). Interrogation of the publicly available ENCODE and Roadmap databases (figure 2) shows that rs4265380 is close to potential regulatory sequences active in CD14+ monocytes that include: (1) DNase I hypersensitivity sites; (2) H3K4Me1, H3K27Ac and H3K79Me2 histone modifications that are associated with regulatory regions and replication initiation sites and (3) ChIP sequencing peaks for TF binding, particularly p300 (a transcriptional activator and histone reader). This region, therefore, appears to contain a potential regulatory element upstream of RUNX3 specifically active in monocytes, in which the C allele at rs4265380 increases the risk of AS.
Haplotype analysis showed the presence of a specific asassociated haplotype We included the two SNPs-rs4265380 and rs4648889that were both independently associated with AS in a haplotype analysis (table 1). Both SNPs were formally genotyped rather than imputed. The CA haplotype was associated with a particularly high risk for AS (OR 6.2; 95% CI 3.1 to 13.2, P=1.4×10 −8 ) compared with the most common CG haplotype. The OR of CA relative to TA is 5.35 (P=2.7×10 −7 ).

rs4265380 affects monocyte cell counts
We interrogated the Roslin Gene Atlas (http:// geneatlas. roslin. ed. ac. uk/), a public database including more than 400 000 individuals (from the UK Biobank), to investigate a plausible role for rs4265380 associations with leucocyte cell counts. The results for rs4265380 showed a significant reduction in monocyte (P=0.006), neutrophil (P=6.4×10 −7 ) and eosinophil (P=8.3×10 −7 ) counts associated with the C allele. No association was apparent for lymphocyte counts (P=0.2), consistent with our hypothesis that rs4265380 is acting specifically in monocytes.

RUNX3 expression is downregulated in as cases compared with healthy controls
The historical data on RUNX3 mRNA in PBMCs measured by RNA-seq showed significantly reduced expression (P<0.002) in AS cases compared with healthy controls. When these results were stratified on rs4265380, there was no apparent influence of this SNP on expression levels (online Supplementary figure 1A).
The region spanning rs4265380 shows regulatory activity in a reporter gene assay We performed luciferase reporter assays to evaluate possible enhancer activity of the region spanning rs4265380. Reporter constructs including the 250-bp DNA sequence spanning rs4265380 increased luciferase activity approximately 10-fold (P=0.0005) compared with the minimal promoter alone, which was consistent with its having enhancer activity. Since there was no apparent difference between constructs containing either the rs4265380 C (risk) or T (protective) allele in this system Figure 3 The genomic region spanning rs4265380 shows enhancer activity. The reporter activity of rs4265380 compared with minimal promoter construct (minP) and normalised to Renilla internal control was measured by luciferase assays on HEK293T (n=3). The values of relative luciferase activity are expressed as mean±SEM. (P=0.0005, Student's t-test). AS, ankylosing spondylitis.

RMD Open
(figure 3), we chose to investigate possible allelic differences by other methods. rs4265380 shows differential allelic binding of monocyte nuclear extract that involves the p300 transcription activator We performed EMSAs to evaluate the impact of the rs4264380 polymorphism on DNA/TF binding. A major protein/DNA complex was seen with both monocyte (figure 4A arrowed) and U937 nuclear extracts (figure 4B arrowed). Little difference in binding was seen between the two rs4265380 alleles in unstimulated cells (online Supplementary figure 2A arrowed), but when stimulated monocyte and U937 cell nuclear extracts were used, a much stronger band was seen for the T (protective) allele. Successful competition with 100-fold excess Spondyloarthritis of unlabelled probe confirmed the specificity of these experiments ( figure 4A and B and online Supplementary figure 2A). Quantification of the DNA/protein complex is shown in figure 4C (n=6 experiments, monocytes+LPS: P=0.02; U937+PMA: P=0.03). Preincubation with p300 antibody also greatly reduced the binding of nuclear extracts from both stimulated monocyte and U937 cells, suggesting that p300 is included in the complex ( figure 4A and B). We conducted ChIP-qPCR (n=3) on LPS-stimulated CD14+ monocytes to assess the enrichment of p300, which was found to be similar for both alleles ( figure 4D). Additional experiments performed with CD8+ T cells and stimulated Jurkat nuclear extracts (online Supplementary figure  2A) showed no differential binding thereby confirming that the effect of rs4265380 appears to be restricted to a monocyte-like milieu.
Other cis-interactions with the RUNX3 locus The RUNX3 locus lies close to two non-coding RNAs (miR-4425 and RP4-799D16.1 lincRNA) and also the SYF2 locus, which is 240 kb upstream of the RUNX3 promoter. By interrogating the Capture Hi-C plotter database (https:// www. chicp. org/ chicp), 25 we observed interactions between RUNX3 bait (overlapping rs4265380) and other regions of chromosome 1, specifically with RP4-799D16.1 and miR4425, but not with the more distant SYF2 locus ( figure 5A). Expression of RP4-799D16.1, a lincRNA with unknown function, showed a significant (P=0.048) increase in LPS-stimulated CD14+ monocytes from AS cases with the CC (AS-risk genotype) compared with the TT genotype ( figure 5B). miR-4425 was not detected (data not shown). There was no discernible effect on the expression of SYF2 (data not shown). rs4265380 AS-risk C allele alters enrichment for specific histone marks in stimulated monocytes.
We assessed the enrichment of specific histone marks associated with active chromatin in this region by ChIP-qPCR (n=3) on LPS-stimulated CD14+ monocytes. There was relative enrichment for H3K79Me2 (average ~3.1fold, P=0.04) and H3K27Ac (average ~4.6-fold, P=0.03) in the presence of the C (risk) allele (online Supplementary figure 3A), while H3K4Me1 was unchanged. No significant allelic difference was seen in unstimulated monocytes (data not shown).

dIsCussIOn
Finding mechanistic explanations for the multiple genetic associations of common diseases, such as AS, is challenging. A small minority of such SNP associations alter protein function that have yielded insights with potential therapeutic implications. Thus, the AS association with rs11209026, encoding an R381Q loss of signalling function dimorphism in the cytoplasmic tail of the IL23 receptor, highlights the importance of the IL23/ IL17 axis in AS. 6 28 This influenced the successful introduction of anti-IL17A biologics for the treatment of AS. 29 However, most GWAS 'hits' lie in non-coding regions of the genome with little or no immediately obvious functional significance. 9 30 In a cross-disease analysis, Ellinghaus et al identified 113 possible non-HLA-B27 AS associated loci, of which only 13 were exonic. 31 Furthermore, of the 39 AS associations in this study, only one key SNP was in an exon. 19 It seems likely that at least some of these associated SNPs exert regulatory functions either on individual genes or gene networks but unravelling these relatively subtle effects will be difficult, 8 32 and the RUNX3 locus is unlikely to be an exception.
Previously, we showed that rs4648889, an AS-associated SNP in this regulatory region upstream of RUNX3, influenced TF binding and gene transcription in CD8+ T cells. 11 Here, we have demonstrated additional complexity in the RUNX3 association with AS. Although rs4265380 is only ~500 bp from rs4648889, not only is its association with AS independent of other AS-associated SNPs in this region but also it likely acts in a different cell type. DNA sequences in the vicinity of rs4265380 show regulatory activity (demonstrated by luciferase reporter activity in vitro) and allele-specific effects on TF binding and histone modifications that could potentially alter gene expression. The relevant TF/DNA complexes include p300, a histone acetyltransferase associated with open chromatin that is involved in transcriptional activation. 33 34 Such coactivators may not bind directly to DNA in a sequence-specific manner, but can function as histone modifiers (eg, histone acetyltransferases) or chromatin remodellers as part of the TF/DNA complex. 35 Our data suggest that the rs4265380 association with AS is probably explained by operational differences in monocytes rather than CD8+ T cells. First, interrogation of epigenetic databases shows relevant areas of open chromatin, TF binding and enrichment for chromatin markers of active transcription limited to monocytes. Second, the allele-specific effects on TF binding appear to be restricted to nuclear extracts from monocytes (particularly after stimulation with LPS).
In the EMSA studies, binding of p300 to rs4265380 was allele specific, but ChIP in heterozygous C/T monocytes from AS cases indicates that p300 actually binds to both rs4265380 alleles in the context of native chromatin. The allelic difference we observed in EMSA is unlikely to be driven by p300 alone but more likely by a multicomponent TF complex recruiting the relevant transcriptional machinery. It is not yet clear whether these findings translate to influences on the expression of nearby genes, including RUNX3. The complexity of the RUNX3 associations and their multiple potential physical chromosomal interactions is further highlighted by the possible interaction of RUNX3 bait (including rs4265380) with a neighbouring non-coding RNA (RP4-799D16.1 lincRNA) of unknown function (see Hi-C plot in figure 5).
The observation that PBMCs from AS cases appear to have reduced RUNX3 expression may be relevant to the pathogenesis of AS. However, it does not correlate with the RUNX3 genotypes that we have described here.

RMD Open
We cannot exclude a disease-specific effect on RUNX3 expression or an effect arising elsewhere in the genome. In contrast, analysis of the UK Biobank data on the Roslin Gene Atlas database showed a significant association of rs4265380 with monocyte, neutrophil and eosinophil counts that might be relevant to an immune-mediated disease like AS. Further evaluation of these observations in well-defined cell types from AS cases seems justified.
We have also identified an extended RUNX3 AS-risk haplotype that includes the independent risk alleles for rs4265380 and rs4648889. 11 The finding of two independent AS-associated SNPs in such close proximity possibly operating in different cell types highlights a role for

Spondyloarthritis
RUNX3 in AS and perhaps suggests that both CD8+ T cells and monocytes may be involved in the pathogenesis of AS. At least in principle, both these cell types could be targeted in the treatment of AS, analogous to the use of anti-CD20 monoclonal antibody therapy targeting B cells in the treatment of rheumatoid arthritis. 36 In summary, we have identified a regulatory region upstream of RUNX3 that could be relevant to the pathogenesis of AS. The association of this region with AS is complex, including at least two neighbouring SNPs. There is some evidence that these have distinct cellular functional profiles, possibly indicative of pathogenic roles for both CD8+ T cells (shown previously) and monocytes in AS. However, substantial further work will be necessary to confirm this.
translational genomics group, institute of Health and Biomedical innovation, School of Biomedical Sciences, Queensland University of technology at translational research institute, Queensland, australia acknowledgements We wish to thank all the study participants who generously donated their Dna to this study. the authors wish to thank Maria cristina carena for technical support.