Epistatic interactions between killer immunoglobulin-like receptors and human leukocyte antigen ligands are associated with ankylosing spondylitis

The killer immunoglobulin-like receptors (KIRs), found predominantly on the surface of natural killer (NK) cells and some T-cells, are a collection of highly polymorphic activating and inhibitory receptors with variable specificity for class I human leukocyte antigen (HLA) ligands. Fifteen KIR genes are inherited in haplotypes of diverse gene content across the human population, and the repertoire of independently inherited KIR and HLA alleles is known to alter risk for immune-mediated and infectious disease by shifting the threshold of lymphocyte activation. We have conducted the largest disease-association study of KIR-HLA epistasis to date, enabled by the imputation of KIR gene and HLA allele dosages from genotype data for 12,214 healthy controls and 8,107 individuals with the HLA-B*27-associated immune-mediated arthritis, ankylosing spondylitis (AS). We identified epistatic interactions between KIR genes and their ligands (at both HLA subtype and allele resolution) that increase risk of disease, replicating analyses in a semi-independent cohort of 3,497 cases and 14,844 controls. We further confirmed that the strong AS-association with a pathogenic variant in the endoplasmic reticulum aminopeptidase gene ERAP1, known to alter the HLA-B*27 presented peptidome, is not modified by carriage of the canonical HLA-B receptor KIR3DL1/S1. Overall, our data suggests that AS risk is modified by the complement of KIRs and HLA ligands inherited, beyond the influence of HLA-B*27 alone, which collectively alter the proinflammatory capacity of KIR-expressing lymphocytes to contribute to disease immunopathogenesis.


Introduction
Ankylosing spondylitis (AS) is an immune-mediated arthritis in which inflammation targeting particularly the pelvis and spine contributes to joint erosion and reactive bone deposition. The disease is strongly associated with inheritance of the human leukocyte antigen (HLA) class I allele HLA-B � 27 (carried by >80% of patients) [1,2], as well as polymorphisms in the endoplasmic reticulum aminopeptidases (ERAP1 and ERAP2) [3][4][5] involved in trimming endogenous peptides for HLA class I presentation. Despite a thoroughly characterised genetic architecture, implicating >100 loci in disease pathogenesis [3][4][5][6][7][8][9], the immunological mechanisms underlying AS are not fully understood. Of interest is the potential role of killer immunoglobulin-like receptors (KIRs) in disease, a diverse collection of paired signalling receptors expressed predominantly on natural killer (NK) cells (and some T-cell subpopulations) that exhibit variable specificity for class I HLA ligands [10][11][12]. Opposing inhibitory and activating signals transduced through surface KIRs buffer the threshold of activation of KIR-expressing cells, particularly NK cells, serving to quench or promote innate killing activity as required to maintain immune homeostasis and control inflammatory responses. The hypervariable nature of KIRs and their HLA ligands bestows upon the human population a large spectrum of proficiencies in lymphocyte responses to activating stimuli and self-tolerance. Accordingly, KIR-HLA co-inheritance has been associated with various immune-related phenotypes, including infection outcome [13][14][15], and susceptibility to autoimmunity [16][17][18] and cancer [19,20].
Fifteen KIR genes have been identified in humans, all of which share significant sequence homology (85-99% similarity) attributed to the duplication and progressive evolution of a common ancestral gene [21,22]. The KIR locus (chromosome 19q13.4) exhibits extreme copy number variation, with genes inherited in haplotypes of diverse content that often carry duplication, deletion and hybridisation events [22,23]. Classically, KIR haplotypes have been categorised into two groups; group A haplotypes containing the single activating receptor gene KIR2DS4 among six inhibitory receptors, and B haplotypes carrying a variable number of both activating and inhibitory receptor genes [24]. Stochastic expression of inherited KIRs is largely controlled by variable promoter methylation, with the vertical transmission of methylation 70], with the KIR3DL2 receptor upregulated on CD4+ T-cell populations in spondyloarthritis patients and shown to prompt differentiation in to pathogenic Th17 cells upon B27 2 ligation [71]. Intriguingly, loss of HLA-B � 27:05 recognition by KIR3DL1+ NK cells drives target cell killing in peptide-specific contexts [34], posing the hypothesis that features of a disease-specific peptidome may disrupt canonical inhibitory KIR interactions. In support of this, disease-associated polymorphisms in ERAP1 alter class I peptide production and destruction and exhibit strong genetic epistasis with HLA-B � 27 [4]. Alternatively, it is plausible that the complement of coinherited KIR and HLA alleles collectively modifies AS risk by shifting the activation threshold of lymphocyte populations, exacerbating autoinflammation.

Ethics and sample acquisition
The IGAS test cohort comprised 8,107 AS cases and 12,214 healthy controls of Caucasian decent, originally recruited by the IGAS Consortium as reported in Cortes et al. 2013 [5]. Written informed consent was obtained from all participants with research ethics approval granted by the relevant ethics committee at each participating centre [5]. All cases had an AS diagnosis as defined by the modified New York criteria [74]. A semi-independent replication cohort comprised 14,844 controls sourced from the UK Biobank (UKB) public resource (project 21024), and the subset of 3,497 IGAS AS cases recruited from the UK, also included in the test cohort. UKB controls, aged between 40-69 years, were selected through exclusion of individuals coding as having one of the following disease states: ankylosing spondylitis (code 1313), inflammatory bowel disease (code 1461), Crohn's disease (code 1462), ulcerative colitis (code 1463), psoriasis (code 1453) or spine arthritis/spondylitis (code 1311); a kinship coefficient > = 0.0442 and non-Caucasian ancestry, with 20,000 of the remaining identifiers selected by random number generation for study inclusion prior to quality control filtering as described below and in Genotyping IGAS test cohort samples were genotyped using the Illumina Immunochip array on the Illumina Infinium platform as previously described [5]. UKB replication cohort controls were genotyped on the UKB Affymetrix Axiom array. Identity by decent was calculated using thegenome command in PLINK [75], with exclusion of one individual from each pair with a PI_HAT score >0.05. Principal components (PCs) for ethnicity confirmation and population stratification correction were calculated for test and replication cohorts independently based on 20,783 and 12,485 autosomal SNPs respectively outside of long-range LD regions [76]. Homogeneity of ethnic background was confirmed by visualisation of the first two PCs with exclusion of individuals falling beyond plus or minus three standard deviations from the mean of the European sample cluster. PCs were recalculated for the filtered European test and replication cohorts (S1 Fig shows stratification of the European cluster along the first two principal components) and the first ten principal components fitted as covariates in all regression models.

KIR imputation
The KIR � IMP software [72] Table. HLA imputation HLA class I alleles were imputed at four-digit subtype resolution using HLA � IMP:03 [80]. Prior to allele imputation, genotypes were converted to Hg19 positions and Illumina plus strand orientation and the Michigan Imputation Server [81] was used to SNP impute HLA locus genotypes (Chr6: 20Mb-40Mb, Hg19) against the 1000 Genomes Phase 3v5 reference panel, with phasing using SHAPEIT [77]. The returned imputed SNP haplotype and sample files were passed to HLA � IMP:03 for HLA imputation, which outputs the presence or absence of each allele in each individual with a corresponding posterior probability of imputation accuracy. The HLA subclass of each allele (i.e. HLA-Bw4, Bw6, C1 or C2) was assigned based on known allele classifications ( Table 1). The number of HLA alleles satisfying the following classes were summed per individual for use in statistical analyses: HLA-Bw4, HLA-Bw4A, HLA-Bw4B, HLA-Bw4B(I80), HLA-Bw4B(T80), HLA-C1 and HLA-C2. KIR3DL1/S1 functional group typing DNA previously extracted from whole blood or saliva was available for 236 HLA-B � 27+ AS patients and 99 HLA-B � 27+ control individuals from the test cohort. DNA concentrations were checked via Qubit fluorometric quantitation and normalised to 10ng/μL in UltraPure H 2 O. Five separate PCR reactions designed to segregate KIR3DL1 alleles into six distinct functional groups were performed on each sample as per the protocol published by Boudreau et al. [73]. In brief an array of primers were designed to target polymorphic sites within three exons of the KIR3DL1 gene (exon 3, 4 and 7), combinations of which stratify alleles into two highly expressed (High-1 and High-2), two lowly expressed (Low-1 and Low-2), null and activating (KIR3DS1) groups (S3 Table). 50ng of input genomic DNA was amplified across each 25uL PCR reaction with reagent concentrations and cycling temperatures as described [73]. PCR products were run on a 1.5% agarose gel with 0.025μL/mL ethidium bromide at 125V for 40 minutes and visualised under UV light. Individuals were assigned a dosage between 0 and 2 for each of the six KIR3DL1/S1 allele groups. Nine individuals who typed positive for three KIR3DL1 alleles, likely due to locus duplication, were excluded from further analysis.

Statistical analysis
All statistical analyses were conducted using custom scripts in R [82]. LD across the KIR locus was calculated using the 'LD' function from the package genetics [83] based on test cohort imputed haplotypes. A generalised linear model (GLM) with logarithmic link (logistic regression) was used to test for an association between imputed haplotype calls or KIR gene dosages and disease status, with inclusion of the first 10 principal components calculated per cohort for population stratification correction. When testing for an interaction between KIR and HLA subclasses or alleles an interaction term was included in the model. KIR dosage was assessed under a dominant (0 = gene absent, 1 = one or more gene copy) and recessive (0 = dosage less than two, 1 = homozygosity at locus) inheritance model. Given that some KIR genes harbour alternate alleles that differ in their activating or inhibitory classification, and genes are strongly linked in both positive and negative LD, dominant and recessive inheritance models were considered the most relevant when assessing the biological consequences of differing KIR dosages. HLA status was treated as dominant. All analyses were conducted in both test and replication cohorts, however replication cohort results are shown only for the top gene associations and epistatic interactions returned from the test cohort. Where applied, correction for multiple testing was performed using the Bonferroni method with statistical significance defined as analyses achieving P<0.05. For interaction P-values, multiple testing correction was applied across all tested KIR-HLA interactions, for KIR association P-values, multiple testing correction was applied across all tested KIRs within the specified HLA subtype or allele group. It is acknowledged that there is presently no established convention for a standard of evidence when reporting genetic interactions (whereas the genome-wide significance threshold is used for single-variant associations). We have thus used the conventional P-value threshold of P<0.05 following multiple testing correction (which controls the family-wise type 1 error rate) to highlight genetic epistatic interactions with greater evidence of disease association.
The association of ERAP1 SNP rs30187 with disease risk was assessed under an additive model using a logistic regression in KIR3DL1+ or-and KIR3DS1+ or-cohorts, defined based on genotype at the top KIR3DL1/S1 tag SNP rs592645 [72] (A allele = KIR3DL1 present, T allele = KIR3DL1 absent). An additional 26,607 HLA-B � 27+ and 299,267 HLA-B � 27-Caucasian controls from the UKB were added to the test cohort to boost sample size. Disease associations with KIR3DL1 functional groups were assessed using a GLM as previously, under both dominant and recessive modes of inheritance.

Study cohorts
Immunochip genotyping array data for 10,619 AS cases and 15,145 controls of mixed ethnicity from the IGAS Consortium, and UK Biobank (UKB) Axiom chip genotyping array data for 20,000 unrelated UKB controls of Caucasian ancestry without AS or associated inflammatory disease, were available prior to sample filtering and imputation quality control (Fig 1). We identified and excluded samples with cryptic relatedness, non-Caucasian ancestry, poor KIR locus SNP clustering or poor KIR � IMP imputed KIR haplotype posterior probability scores. This left 8,107 AS cases and 12,214 healthy controls in the IGAS test cohort, and 3,497 IGAS AS cases and 14,844 UKB controls in the replication cohort. As the 3,497 AS cases in the replication cohort were selected from the British AS cases in the IGAS test cohort, to match the ethnicity of the UKB controls, the replication dataset is considered only semi-independent. KIR haplotypes for these individuals were re-imputed using the reduced SNP set available for the UKB control samples (Fig 1). KIR gene and haplotype frequencies taken from imputed test and replication cohorts were largely in accordance with those reported for published European populations with dosages typed by laboratory-based methods. However, implementation of imputation quality control led to a slight bias in both cases and controls toward A group haplotypes, due to their lower copy-number variability and higher imputation accuracy relative to B group haplotypes (S2 Table, S3 Fig). A direct comparison of KIR dosages for 32 Centre d'Etude du Polymorphisme Humain (CEPH) and 20 control samples for which both imputed data and laboratory KIR typing was available revealed >90% concordance across loci with a per gene imputation posterior probability threshold of > = 0.4, and 100% at a threshold of > = 0.6 (S4 Table). All sample data, including disease status, imputed KIR gene and HLA allele dosages, HLA subtype counts and principal components are included in S5  S4 Fig) and are inherited together with KIR2DL5B as part of the telomeric tB01 motif in a number of common B haplotypes (B3, B7 and B8, S1 Table). Under a recessive inheritance model, homozygosity for KIR3DL1 and KIR2DS4TOTAL (R 2 = 1) was nominally associated with disease protection, and homozygosity for KIR2DS5 with disease risk ( Table 2). Haplotypes B3 and B11 were seen at increased frequency in HLA-B � 27+ AS cases ( Table 2). Multiple testing correction ablated KIR gene and haplotype associations in the test cohort and associations were not evident in the replication cohort.

KIR interactions with HLA subclass ligands
The association of co-inherited KIRs and HLA subclass ligands with AS risk was assessed using an interaction term in a generalised linear model. A significant KIR-HLA interaction term indicates that the association of a given KIR with disease risk is seen particularly in the context of a specific HLA subclass. HLA subclass ligand count was treated as dominant and KIRs assessed under both dominant and recessive modes of inheritance. Only interactions between KIRs and HLA subclass ligands known to biologically interact, and showing evidence of a disease association, are reported in Table 3, full results (showing KIRs in LD with reported genes and interactions between receptors and ligands without present functional evidence) are reported in S7 Table. There was evidence for a statistical interaction between KIR3DL1 and HLA-Bw4A (HLA-A alleles with the Bw4 motif) in both test and replication cohorts. Carriage of KIR3DL1 in those who had inherited at least one HLA-Bw4A subclass allele was associated with increased disease risk, whereas no association was observed in HLA-Bw4A negative subjects ( Table 3). Nominal interactions were also observed between KIR3DL1 and HLA-Bw4B(I80), whereby in the absence of an HLA-Bw4B(I80) ligand KIR3DL1 homozygosity was associated with increased disease risk ( Table 3). There was no evidence of an interaction between KIR3DL1 and HLA-Bw4B(T80) group ligands, of which the most common AS-associated allele in European populations, HLA-B � 27:05, is a member. There was also evidence of an interaction between the activating receptor gene KIR2DS1 and alleles of the HLA-C2 subclass, encoding canonical KIR2DS1 ligands. Carriage of KIR2DS1 in the absence of an HLA-C2 allele (i.e. HLA-C1 homozygosity) was associated with reduced disease risk in the test cohort only ( Table 3).

KIR interactions with HLA class I alleles
Several KIR-HLA subclass interactions appeared to be driven by specific HLA class I alleles ( Table 4). KIR3DL1 exhibited a statistical interaction with the HLA-Bw4A allele HLA-A � 32, with odds of disease increased >3 fold in those co-inheriting both the receptor and ligand. The association of KIR3DL1 in HLA-A � 32+ individuals remained evident following multiple testing correction in the test cohort and was independently detected in the replication cohort. The absence of HLA-C � 05 (HLA-C2) specifically was associated with disease protection in KIR2DS1 carriers. Alternatively, in the HLA-C � 05+ cohort, carrying the deletion isoform of KIR2DS4 (KIR2DS4D) was also associated with disease protection, although evidence was lost with multiple testing correction. KIR2DS1/4-HLA-C � 05 interactions were not evident in the replication cohort. KIR2DS2/L2/L3 (KIR2DS2 and KIR2DL2 are in strong LD, R 2 = 1 S4 Fig, and KIR2DL2 and KIR2DL3 are alternate alleles of the same gene) showed evidence of an interaction with the HLA-C1 alleles HLA-C � 12, HLA-C � 07 and HLA-C � 08 in the test cohort, though it was not possible to determine which of the tightly linked KIRs were driving the interaction effect. Interactions with HLA-C � 12 were the strongest observed and of the nature that HLA-C � 12+ individuals homozygous for KIR2DS2 and KIR2DL2 were at a two times increased risk of disease, whereas those carrying at least one copy of KIR2DL3 were protected ( Table 4). KIR associations tested in HLA-C � 12+ individuals remained significant at P<0.05 upon multiple testing correction in the test cohort however trends were less evident in the smaller replication cohort. There was a modest interaction observed between KIR3DL1 homozygosity and HLA-B � 27 carriage in the test cohort alone, through statistical significance was lost with multiple testing correction. HLA-B � 27+ individuals who were homozygous for KIR3DL1 showed protection from disease, whereas there was no evidence of a protective association in the HLA-B � 27cohort. Interactions between receptors and HLA alleles of the canonical ligand subclass are reported in Table 4, with full results (showing KIRs in LD with reported genes and interactions between receptors and ligands without present functional evidence) reported in S8 Table. KIR

interactions with ERAP polymorphisms in HLA-B � 27+ AS
Given that the HLA-B � 27-bound peptide can disrupt the affinity of KIR3DL1 for HLA-B � 27 ligands, the lead AS associated SNP in ERAP1 (rs30187), known to alter peptide trimming by  Fig 2), providing no support for a three-way interaction between AS-associated ERAP1 genetic variation, KIR3DL1/S1 and HLA-B � 27 influencing AS susceptibility. The reduction in the significance of the rs30187 association in the HLA-B � 27+ KIR3DL1-cohort can likely be ascribed to the low number of AS cases in this group. As reported previously, there was no disease-association of rs30187 in HLA-B � 27-individuals, irrespective of KIR3DL1/S1 dosage (Fig 2).

KIR3DL1/S1 allele group associations with HLA-B � 27+ AS
The PCR protocol designed by Boudreau et al. (2014) was used to segregate KIR3DL1 alleles from 236 HLA-B � 27+ AS patients and 99 HLA-B � 27+ controls into six distinct functional groups (Null, High-1, High-2, Low-1, Low-2 and KIR3DS1) based on a combination of polymorphic sites that segregate with surface expression phenotype [73]. High and low frequency alleles captured within each functional group are listed in S3 Table. Functional group frequencies in the control cohort were in accordance with those reported in the published method (S3 Table). When assessed under a dominant or recessive mode of inheritance, there was no significance difference in the frequency of functional cell-surface expressed (KIR3DL1 � f: High-1, High-2, Low-1 and Low-2 alleles) or unexpressed (KIR3DL1 � n: Null) KIR3DL1 alleles between HLA-B � 27+ve AS patients and controls. Similarly, there was no statistically significant difference in allele group frequency between patients and controls ( Table 5). All data, including disease status and KIR3DL1/S1 functional group typing is included in S9 Table.

Discussion
Assessing the contribution of KIR and HLA co-inheritance to an immune-related phenotype is inherently challenging, with immense genomic diversity at both loci converging to influence the dynamics of an individual's lymphocyte responses. Very large cohorts are required to identify KIR associations with human disease amid the genetic 'noise' ascribed to copy-number variable haplotypes, and account for the modifying effect of independently inherited HLA ligands on KIR signalling. Here we have used imputation methods to type HLA class 1 alleles [80] and gene dosage level KIR haplotypes [72] from SNP genotype data in order to perform  79]. However, removal of haplotypes with low (<0.4) posterior probability of imputation accuracy largely excluded individuals with poorly imputed rare gene content B group haplotypes. This led to a slight A group haplotype bias in both test and replication cohorts, hindering the ability to dissect genetic associations attributed to closely linked genes by conditional analysis (which are more likely to be found in rare arrangements on uncommon B group haplotypes). Albeit, given our sizeable dataset of high-confidence imputed haplotypes, we were able to detect suggestive interactions between KIR genes and HLA ligands at both subclass and allele resolution showing association with AS risk. It is possible that such interactions modify the activation threshold of lymphocyte populations in a disease setting, contributing to or sustaining a damaging autoinflammatory environment. The profound association of the HLA-Bw4 allele HLA-B � 27 with AS has encouraged investigation into a role for KIR3DL1/S1 in disease; the inhibitory KIR3DL1 receptor being the only KIR known to engage canonical HLA-Bw4 ligands [40]. Its activating allele, KIR3DS1, has been demonstrated to mediate NK killing of HLA-B � 57 (HLA-Bw4) target cells in vitro in a peptide dependent fashion [47], of potential additional interest in AS given the modest protective association of HLA-B � 57:01 with the disease [7], and may recognise other HLA-Bw4 ligands with similar peptide specificity. We identified a modest increase in the frequency of KIR3DS1 in patients (29.2% in controls, 32.5% in AS) and genes KIR2DL5, KIR2DS5 and KIR2DS1 in the HLA-B � 27+ test cohort, all of which occur together on the telomeric half of haplotype B3 [23], which also showed a modest association with disease risk. Association of KIR3DS1 with HLA-B � 27+ AS has been detected in numerous studies [57-59, 61, 64], with some studies reporting large KIR3DS1 frequency differences between AS cases and healthy controls (e.g. 26.6% in controls and 64.3% in AS in a HLA-B � 27+ Chinese cohort of 72 individuals [58], and 40.4% in controls and 59.2% in AS in a HLA-B � 27+ Spanish cohort of 705 individuals [61]). Despite studying a much larger case-control dataset, we did not find a markedly increased frequency of KIR3DS1 in AS patients in either the test or replication cohorts. Although differences in ethnic background may account for disparate results across studies, our data suggests that any association of KIR3DS1 with AS must have quite a small overall effect size.
Conversely, homozygosity for KIR3DL1 was nominally associated with disease protection in the HLA-B � 27+ test cohort, replicating previous studies demonstrating that the recessive genotype reduces disease risk [57,58,61,64]. Despite being an HLA-Bw4T80 ligand with weak KIR3DL1 binding affinity relative to HLA-Bw4I80 alleles, the high surface density of HLA-B � 27:05 (the predominant subtype in European HLA-B27 carriers) makes it a strong inhibitor of KIR3DL1+ NK cytotoxicity [27]. It could be hypothesised that the increased proportion of KIR3DL1+ NK cells in KIR3DL1 homozygotes minimises the risk of undue lymphocyte activation and autoreactivity in healthy HLA-B � 27 carriers. No KIR3DL1 association was seen with HLA-B � 27-disease, resulting in a detectable epistatic interaction between the two genetic factors in the test cohort, however the interaction was weaker than that seen between a number of other KIR-HLA pairs and significance was ablated with multiple testing correction.
It is possible that a protective effect of KIR3DL1 in HLA-B � 27+ AS is confined to a specific allelic variant of the gene (for instance, engagement of HLA-B � 27:05 inhibits NK cytotoxicity through KIR3DL1 � 002 but not KIR3DL1 � 007 [29], and a point mutation in KIR3DL1 � 004 disrupts protein folding and abolishes cell surface expression of this allotype altogether [84]), which would be obscured in gene-dosage resolution analyses. To assess this, we used a PCR approach in a subset of the HLA-B � 27+ test cohort to type and stratify KIR3DL1 allotypes into six independent functional groups differing in surface expression and binding affinity [73]. We observed no disease-association with KIR3DL1 genotype irrespective of whether functional (cell surface expressed) KIR3DL1 � f or non-functional (KIR3DL1 � 004) KIR3DL1 alleles were carried. Such finding differs from data presented by Zvyagin et al. (2010), showing disease protection afforded by carriage of the KIR3DL1 � f allele and the KIR3DL1 � f homozygous genotype, but no association with KIR3DL1 � 004 in a HLA-B � 27+ Caucasian cohort of 83 patients and 107 controls [64]. We also observed no significant disease-association with any specific KIR3DL1 functional group, noting the power limitations given the modest sample size of this component of the study. This and previous studies have been underpowered to deconvolute KIR3DL1 allelic associations at this locus [61], largely given the hyperpolymorphic nature of the receptor (for which 147 independent alleles are currently known [67]). Evidently, highthroughput KIR typing and allele imputation approaches need to advance in line with those available for the HLA locus if appropriately powered, high-resolution allotype association tests are to be feasible.
Additional to both the KIR allele and HLA subtype present [85], the outcome of KIR3DL1 signalling upon recognition of an HLA-B � 27 ligand is known to be modified by the HLA bound peptide [34]. There is substantial genetic and functional evidence supporting a role for altered peptide processing in AS, particularly given a non-synonymous coding variant in ERAP1, rs30187, is uniquely associated with HLA-B � 27+ disease [4,5] and modifies the HLA-B � 27 peptidome [86,87]. We tested the hypothesis that the rs30187 risk-association with HLA-B � 27+ disease is dependent on the KIR3DL1/S1 background, with altered peptides either perturbing inhibitory signalling through KIR3DL1 or promoting NK activation through KIR3DS1 to enhance pathogenic lymphocyte responses. Disproving this, the rs30187 association with HLA-B � 27+ AS was clearly detectable in both KIR3DL1-and KIR3DS1-individuals. Although it cannot be disregarded that the altered HLA-B � 27 peptidome shaped by diseaseassociated ERAP1 alleles may influence KIR3DL1/S1 signalling dynamics and subsequently lymphocyte activity, this does not appear to be the immunological mechanism by which the strong ERAP1-HLA-B � 27 epistasis in AS arises.
Dissecting genetic interactions between inherited KIRs and HLA class I alleles beyond HLA-B � 27 requires a large enough cohort to retain a degree of statistical power when stratifying by both KIR and HLA carriage (some combinations of which are found at low frequency in the population). This study presents the first attempt to test for associations of KIR-HLA epistasis with an immune-mediated disease across all represented gene-gene combinations. Despite the substantial sample size of our test cohort, lack of an equally well powered replication cohort of independent AS patients made validation of findings difficult. Further, given strong LD across both the HLA and KIR locus, members of a statistically interacting receptorligand pair may not biologically interact but rather tag closely linked genes for proteins that do engage to alter immune cell function. Rare KIR haplotypes, in which unique gene combinations are observed, were likely overrepresented in those discarded during imputation quality filtering, hindering our ability dissect strong LD structures within this data. Here, only those interactions that can be interpreted in terms of experimentally validated KIR-HLA engagements are discussed, though to ascertain the relevance of these findings to AS pathogenesis on a functional level will require further experimentation. For the most part, this data supports a model whereby the contribution of many co-inherited KIR and HLA pairs shifts an individual's predisposition to AS, as is likely the case with many immune-mediated diseases, some combinations imparting more dominant effects than others.
In both test and replication cohorts, carriage of KIR3DL1 in the presence of an HLA-Bw4A allele was associated with increased disease risk, the association remaining statistically significant upon multiple testing correction in the former. At an allelic level, the KIR3DL1 interaction could be detected most strongly with HLA-A � 32, encoding a HLA-Bw4 ligand shown to protect target cells from lysis by KIR3DL1+ NK cells [88]. Why coinheritance of an inhibitory receptor-ligand pair may be associated with risk of, not protection from, a condition associated with chronic inflammation is unclear; HLA-A � 32:01 binds with stronger affinity to some KIR3DL1 alleles than most HLA-Bw4B molecules [27,28]. Increased inhibitory KIR-HLA binding strength has been shown to positively correlate with NK cytotoxicity against HLA-negative target cell lines in vitro [27], and this may exacerbate autoreactivity in instances that HLA expression is downregulated or altered peptides subvert binding in an inflammatory-disease setting. Numerous inflammatory diseases, including AS, have been associated with endoplasmic reticulum stress [89], which contributes to interrupted assembly and surface expression of appropriately folded HLA glycoproteins and increases target cell susceptibility to NK cell lysis [90].
In agreement with previous studies [57,65], coinheritance of KIR3DL1 and HLA-B ligands of the I80 subtype was decreased in AS patients, KIR3DL1 homozygotes lacking an inhibitory HLA-Bw4B(I80) ligand significantly more often than controls in the test cohort. We did not detect a specific allele of the HLA-Bw4B(I80) subclass that appeared to drive this interaction. Strong KIR3DL1 binding to a number of HLA-Bw4B(I80) ligands [27,28] may assist in reducing the activation threshold of NK and T-cell populations under homeostatic conditions in those who inherit both genetic factors. Absence of an HLA-C2 ligand in individuals who had inherited the activating receptor KIR2DS1 was also associated with disease protection in the test cohort. Unlike other activating receptors for which true ligands are unknown, KIR2DS1 recognition of HLA-C2 has been previously shown to activate KIR2DS1+ NK cells (with signalling sufficient to override KIR2DL1 inhibition on KIR2DL1+KIR2DS1+ cells), inducing IFNγ secretion and degranulation in a peptide dependent fashion [37,91,92]. Disease risk attributed to co-inheritance of KIR2DS1 and HLA-C2 has been previously observed in two independent AS cohorts [59,65], as well as in psoriasis vulgaris [93,94] and psoriatic arthritis [18], whereas it is associated with protection from Hodgkin's lymphoma [95] and the anti-leukemic activity of alloreactive NK cells in a HLA-C2 context [96]. Evidently this is a proinflammatory combination of factors. The KIR2DS1-HLA-C2 interaction appeared to be driven most strongly by HLA-C � 05 in this study, perhaps suggesting that this ligand is a particularly potent activator of KIR2DS1+ NK or T-cells. Finally, we detected a number of KIR2D/-HLA-C1 allele interactions associated with AS risk that are yet to be reported, however tight positive and negative LD between receptor genes KIR2DS2, KIR2DL2 and KIR2DL3 make it difficult to discern which receptor may be driving the biological effect. Strong disease-associations with activating and inhibitory KIR2D genes in HLA-C � 12 and HLA-C � 07 carriers in the test cohort, and suggestive associations in HLA-C � 01 and HLA-C � 08 carriers, emphasises that the HLA background likely imparts an array of modifying effects on disease risk beyond those ascribed to HLA-B � 27 alone.
In silico prediction of KIR dosages from genotype data has opened avenues for KIR association testing in cohorts too large for conventional laboratory-based typing methods. Albeit, prediction algorithms inherently operate with some degree of error, and the findings of computational studies should be validated with laboratory methods where possible. Although KIR locus imputation is presently erroneous for rare gene content haplotypes, as large human reference datasets begin to capture the genomic diversity of the KIR locus these approaches will likely gain accuracy in both gene dosage and allelic discrimination. The targeted coverage of the KIR locus by the Illumina Immunochip has enabled us to impute high-confidence KIR dosages in a large test cohort of AS cases and controls for disease-association testing. We have endeavoured to replicate findings in a semi-independent cohort, incorporating UKB controls, however the lower imputation accuracy (ascribed to reduced SNP coverage) and reduced size of this sample set has hindered this somewhat. Further replication of the results presented here are required to ensure accuracy and their biological interpretation will inevitably require functional investigation.
In conclusion, here we report the largest analysis of KIR and KIR-HLA co-associations with any immunological phenotype to date, addressing the contribution of these complex receptors to AS immunopathogenesis. We identified multiple nominally significant epistatic interactions between the genes encoding KIRs and their HLA ligands on both the subtype and allele level, suggesting that AS risk, as likely the case for many immune-mediated diseases, is modified by a mosaic of genetic effects that converge to influence the proinflammatory capacity of KIR expressing lymphocytes. Notably, although we replicated the direction of effect for previously reported KIR3DL1/S1 associations with HLA-B � 27+ disease, we demonstrated that the KIR3D background does not modify the profound HLA-B � 27-ERAP1 epistasis observed in AS. Thus a three-way interaction between these factors is unlikely to be the primary driver of pathogenesis. As the resolution and throughput of KIR typing improves, more studies of this nature will assist in defining patterns in KIR-HLA coinheritance that contribute to complex immune phenotypes, improving understanding of the dynamic role of these receptors in health and disease.  Table. Percent concordance between KIR � IMP imputed gene dosages and direct typing for a subset of control samples (including 32 CEPH samples). Concordance was calculated on an individual rather than a haplotype basis as phasing was undetermined for laboratory typed samples. Percent concordance is reported either for all individuals (n = 52), or just those with all gene dosages imputed with posterior probability above 0.4 (n = 51) or 0.6 (n = 37). KIR2DS4D = KIR2DS4 deletion allele, KIR2DS4W = KIR2DS4 wild-type.  Table. Interactions between KIR genes and HLA class I subclasses (all interactions with interaction P-value = < 0.05). Interactions in blue are those between receptors and HLA subclasses ligands known to biologically interact. Alternate shading indicates groups of KIRs in strong LD that demonstrate statistical interactions with the same HLA subclass ligand. P = P-value denoting significance of the KIR association with AS when assessed in the specified HLA allele group, Int.P = P-value for the KIR-HLA allele interaction term, � = P-values that remain significant (P<0.05) after multiple testing correction, + = dominant inheritance, ++ = recessive inheritance (homozygosity), AS = ankylosing spondylitis, CTRL = control, OR = odds ratio, SE = standard error, NS = not significant, KIR2DS4D = KIR2DS4 deletion allele, KIR2DS4W = KIR2DS4 wild-type. Proportions are the proportion of individuals with the specified KIR genotype in cohorts split by HLA allele carriage. (DOCX) S8 Table. Interactions between KIR genes and HLA class 1 alleles (all interactions with interaction P-value = < 0.05). Interactions in blue are those between receptors and HLA ligands of a subclass known to biologically interact. Alternate shading indicates groups of KIRs (mostly in strong LD) that demonstrate statistical interactions with the same HLA allele. P = P-value denoting significance of the KIR association with AS when assessed in the specified HLA allele group, Int.P = P-value for the KIR-HLA allele interaction term, � = P-values that remain significant (P<0.05) after multiple testing correction, + = dominant inheritance, + + = recessive inheritance (homozygosity), AS = ankylosing spondylitis, CTRL = control, OR = odds ratio, SE = standard error, NS = not significant, KIR2DS4D = KIR2DS4 deletion allele, KIR2DS4W = KIR2DS4 wild-type. Proportions are the proportion of individuals with the specified KIR genotype in cohorts split by HLA allele carriage. (DOCX) S9 Table.  and KIR2DS5 genes are duplicated on some B haplotypes and can occur on both the centromeric and telomeric halves of the haplotype. At present KIR � IMP is unable to distinguish centromeric from telomeric copies of these genes so they have been grouped together as a single locus. The database did not include frequencies for wild type and deletion variants of KIR2DS4. (TIF) S4 Fig. Pairwise LD between KIR genes calculated using KIR � IMP imputed haplotypes from test cohort controls. Strength of LD is coloured according to R 2 , with gene pairs in perfect positive or negative linkage (always or never occurring together in a haplotype) coloured red with an R 2 value of 1. Genes are ordered according to genomic position from KIR2DS2 (centromeric) to KIR2DS4 (telomeric), with exclusion of framework genes KIR3DL3, KIR3DP1, KIR2DL4 and KIR3DL2. Distinction could not be made between centromeric and telomeric copies of KIR2DS3/5 or KIR2DL5.