Classical HLA-DRB1 and DPB1 Alleles Account for HLA Associations with Primary Biliary Cirrhosis

Susceptibility to primary biliary cirrhosis (PBC) is strongly associated with HLA region polymorphisms. To determine if associations can be explained by classical HLA determinants we studied Italian 676 cases and 1440 controls with genotyped with dense single nucleotide polymorphisms (SNPs) for which classical HLA alleles and amino acids were imputed. Although previous genome-wide association studies and our results show stronger SNP associations near DQB1, we demonstrate that the HLA signals can be attributed to classical DRB1 and DPB1 genes. Strong support for the predominant role of DRB1 is provided by our conditional analyses. We also demonstrate an independent association of DPB1. Specific HLA-DRB1 genes (*08, *11 and *14) account for most of the DRB1 association signal. Consistent with previous studies, DRB1*08 (p = 1.59 × 10−11) was the strongest predisposing allele where as DRB1*11 (p = 1.42 × 10−10) was protective. Additionally DRB1*14 and the DPB1 association (DPB1*03:01) (p = 9.18 × 10−7) were predisposing risk alleles. No signal was observed in the HLA class 1 or class 3 regions. These findings better define the association of PBC with HLA and specifically support the role of classical HLA-DRB1 and DPB1 genes and alleles in susceptibility to PBC.


Introduction
The human major histocompatibility complex, HLA, has been implicated in the etiopathogenesis of primary biliary cirrhosis (PBC), similar to many other autoimmune diseases. Genome-wide association studies (GWAS) of PBC including our own find the strongest association with single nucleotide polymorphisms (SNPs) within the HLA region [1][2][3] . In these studies the peak association signal is between HLA-DQA1 and HLA-DQB1. Multiple studies of PBC also show association with particular classical HLA alleles in PBC (reviewed in Invernizzi 4 ). These studies have variably implicated different DRB1 alleles in European populations with most studies including all larger cohorts showing association of DRB1*08 5,6 . Our previous studies in an Italian cohort with PBC showed the association of DRB1*08 as predisposing, and DRB1*11 and DRB1*13 as protective alleles 6 . A study using a small cohort (32 German PBC cases and 47 controls) suggest that DPB1 associations may also be present in Europeans 7 . However, a comprehensive study of HLA region associations has not been performed and like other autoimmune diseases it is unclear which determinants are actually causally related to pathogenesis.
To further study HLA associations in PBC, in the current study we used the most recent advances in imputation algorithms and sequence information resources including the 1000 genome database to accurately impute missing SNPs, and importantly HLA classical alleles. Specifically, our investigation rests on recent development and resources for imputing HLA classical alleles including a reference set of European subjects 8 . For our study we used an inference set of SNP genotypes from both GWAS and a designed chip array, the Immunochip 9 that contains a set of SNPs that have been used in multiple studies of HLA 10 . We perform a series of conditioning analyses that clarify which HLA genes and alleles underlie the major component of the genetic associations of PBC.
As expected from the analysis of classical HLA alleles, PBC also showed strong association with specific AAs in these genes. Most of the HLA AA showing association signals corresponded to the key residues that distinguish the specific classical alleles which for DRB1, included lysine (L) at AA74 in DRB1 (DRB*08), glutamate (E) at AA58 (DRB*11), alanine (A) at AA57 (DRB1*14), and histidine (H) at AA60 (DRB1*14). Similar results were observed for specific DQB1, and DQA1 AAs that are in strong LD with specific DRB1 alleles and AAs (Supporting Table 2).

Conditioning Studies Using Classical HLA Genes
To examine whether these associations could be explained for by known coding differences in genes we next performed a series of conditional analyses. These studies were done by conditioning on a combination of the alleles from an HLA gene (e.g. DRB1) to control for the association that might be attributable to each gene, albeit some of the effect may not be directly attributed to that gene due to extensive LD across this region. The residual signals after controlling for the effect of various combinations of classical alleles and AA residues in these HLA genes show that both DRB1 and DQB1 could account for most of the association signal ( Fig. 1 B-E, Table 2 and Supporting Table 3). In addition, the signal in the DPB1 region was only marginally decreased conditioning on DRB1, DQB1 or DQA1. Conditioning on DPB1 eliminated the signal in the DPB1 region and showed a modest increase in the signal in the across DRB1, DQA1 and DQB1.
To further assess these conditional analyses, we also examined the relative difference of the conditioning by different HLA genes by examining beta estimates and their differences. The beta estimate is the measure of the increase in log-odds that can be attributed to each copy of a given minor allele. The largest effect is from the composite of DRB1 alleles as shown by the residual beta estimates (and odds ratios) after conditioning and the mean change in the beta estimates (Table 3). This is most evident examining the SNPs with the strongest signals from association (original signal < 5 × 10 −8 ). For example, the DRB1 conditioning had a much larger mean change in beta estimate (−0.424) compared with DQB1 (−0.236) (p value < 10 −10 , paired T test). Additional conditional analyses using combinations with DQA1 demonstrated that DQA1 could not substitute for DRB1 or DPB1 in any of the combinations tested (data not shown).

Conditioning on Specific HLA alleles
We next examined the effect of conditioning on specific DRB1, DQA1, DQB1 and DPB1 classical genes and AAs. A clear pattern emerged showing that the association of groups of SNPs was specifically controlled by different alleles. These results are highlighted in Table  4 and in a more complete version Supporting Table 3. Individually, the specific SNPs conditioned a part of the association signal largely corresponding to those SNPs in moderate or strong LD (r 2 > 0.5) with the particular classical specificity or AA. Similar effects were observed for specific alleles in one gene in strong LD with another gene. This is particularly evident for the DRB1*14 and DQB1*05:03 in which the effects of controlling for either of these alleles was virtually indistinguishable (Supporting Table 3). Other pairs (e.g. DRB1*08 and DQB1*04:02) showed small but consistent differences in which the DRB1 allele diminished more signal than when the DQB1 allele was used in the conditioning analyses (Supporting Table 3). Notably the strongest AA association was at position 74 that is in the antigen binding pocket of DRB1 11, 12 . Similar to when we conditioned on genes ("all" alleles at a particular gene), the DPB1 region SNP signal was only substantially reduced when DPB1*03:01 (Table 4 and  Supporting Table 3) or when specific DPB1 AA's were used in conditioning. The strongest effects were observed for the lysine at AA position 11 and the methionine at AA position 76 that are both members of the 16 AA in the putative antigen binding pocket of this gene 11 .
Thus, the vast majority of the HLA region association signal can be accounted individually by conditioning on one of four specific alleles, three in DRB1 (*08, *11, and *14), and one in DPB1 (*03:01). Combinations of these specific alleles accounted for most of the remaining signal (Supporting Table 3 and Fig. 1F) and are also reflected in the strong reduction in beta estimates (Table 3). However, there are signals from several SNPs in the DRB1-DQB1 region that are not accounted for by these conditioning studies. None of these SNPs with signals p <10 −5 after conditioning on DRB*08, *11 and *14, and DPB1*03:01 were among the stronger associated SNPs prior to conditioning (all with original association p values >10 −6 ). In particular, the strongest associated SNP after conditioning (rs9268668, p = 1.67 × 10 −7 ), showed no signal prior to conditioning (p = 0.40). Whether these residual or new signals are also due to other specific classical HLA genes is not clear, however, conditioning on "all" DRB1 and DPB1 alleles ablated all signals with resulting p-values > 10 −5 (Fig. 1F) suggesting that additional sequence differences (e.g. putative regulatory SNPs) do not have to be postulated.
Most of the signal observed for specific AAs was also specifically eliminated when conditioning on the DRB1 or DPB1 classical alleles. However, there were several exceptions in which the association signal was not readily decreased by controlling for single classical HLA alleles. These AAs included DRB1-AA47F, DRB1 AA74A, DQB1-AA26G, and DQB1-AA74S. For these AAs the signal was ablated when conditioning on two DRB1 alleles (DRB1*08 and DRB1*11) (Supporting Table 3). Conversely, conditioning on these AA could not account for the association of the most of the other SNPs that were not ablated by single classical HLA alleles (Supporting Table 3). Therefore, it may be less likely that these particular AAs are critical to explaining the association patterns we observed. However, we cannot exclude a specific functional role for these AAs and it is notable that DRB1-AA47 and DRB1-AA74 are both in the antigen binding pocket of DRB1 11,12 and that conditioning on DRB1-AA47F did ablate several of the association signals that were not controlled by individual DRB1 classical alleles.

Genotypic Associations
We also examined genotypic associations including combinations of susceptibility alleles and combinations of risk susceptibility and protective alleles. Examining individuals with combinations of risk alleles DRB1*08 combined with DPB1*03:01 or DRB1*14 combined with DPB1*03:01 we found higher odds ratios for disease association than when examining only individuals with single susceptibility alleles (Supporting Table 4). There were insufficient numbers of DRB1*08/DRB1*14 (frequency <1%) to evaluate these heterozygote genotypes. The increased risk of the combining DPB1*03:01 with DRB1*08 or DRB1*14 was observed whether DPB1*0301 was on the same or different haplotype as that of the DRB1 risk allele (Supporting Table 4). When risk alleles were combined with the DRB1*11 protective allele the odds ratios were near 1 and there was no significant association with disease.
Finally, we also examined the cumulative combination of risk predisposing and protective alleles (Table 5 and Supporting Table 4). The count of predisposing alleles minus protective alleles showed a strong correspondence with the odds ratio for PBC. Individuals with an excess of one or two or more risk alleles showed an odds ratio (OR) of 3.05 and 5.25 between cases and controls, and conversely individuals with an excess of one or two or more protective alleles had OR of 0.5 or 0.38, respectively ( Table 5). All results were compatible with an additive model of action between each of the alleles similar to our previous observations 6 .

Discussion
The current study of PBC association with HLA differs from previous investigations by providing the most comprehensive analysis of the entire HLA region while correcting for multiple confounding factors. Our results are consistent with a predominant role for class II genes and we believe exclude any substantial effect from either HLA class I or class III genes (there were no residual signals for these genes with p <0.0005 after accounting for class II genes). This contrasts other autoimmune diseases in which HLA class I or class III plays a predominant role (e.g. myasthenia gravis 13 ) or strong class I gene effects are observed independent of class II associations (e.g. type 1 diabetes 14 and multiple sclerosis 15,16 ).
Our study strongly suggests that the major gene in HLA that underlies susceptibility to PBC is DRB1. Overall DRB1 alleles show the strongest associations and conditioning studies show that DRB1 could account for almost all (except DPB1 region) of the association signal. HLA-DQB1 shows association that is only marginally less than that observed for DRB1. However, several points suggest that these associations are secondary to the strong LD between DRB1 and DQB1: 1) the overall strength of association of particular DRB1 alleles is stronger than the corresponding DQB1 allele; 2) conditioning on DRB1 could account for all DQB1 associations; and 3) residual beta-estimates after conditioning showed a substantially stronger DRB1 than DQB1 effect.
In addition, our study provides strong evidence for an independent effect of DPB1. Although previous studies have as indicated in the introduction suggested DPB1 associations in PBC, these were based on small subject sets and were difficult to evaluate. Our study demonstrates that the association of DPB1 cannot be accounted for by controlling for other HLA region genes. These results are also consistent with findings in some, but not other, autoimmune diseases in which an independent effect of DPB1 has been reported. These include juvenile idiopathic arthritis 17 , type 1 diabetes 18 , multiple sclerosis 16,19 and particular autoantibodies in systemic lupus erythematosus 20 .
We note that our study does not directly address whether DRB3, DRB4, DRB5 or structural variations might have additional independent associations. At present such studies are challenging due to absence of reference sets for imputation and/or difficulty in assessing these polymorphisms including whether missing genotypes (excluded SNPs with call rates <0. 95) may have excluded analysis or inclusion of SNPs within these genes in available arrays.
This study has also addressed the association of specific HLA-gene alleles. Most of the HLA association with PBC can be attributed to specific associations with DRB1*08, DRB1*11, DRB1*14 and DPB1*03:01. DRB1*08 has the strongest association, followed by DRB1*11 consistent with several previous studies 6,7 . PBC associations with DRB1*14 has not been previously demonstrated, however, this weaker effect is supported by our conditioning studies that show that this classical allele can control for a set of associated SNPs and AAs that are not strongly influenced by other classical alleles (Table 3 and Supporting Table 3). The DPB1*0301 association is consistent with a previous study of a small German cohort 7 .
In a previous study we observed that DRB1*13 was a protective allele 6 . In the current study the association of DRB1*13 was weak (p = 4.9 × 10 −3 , OR= 0.69, 95% CL = 0.53 -0.89) compared to the previous study (p = 3.6 × 10 −6 ). This may be due to several factors: i) the previous study did not explicitly control for population substructure; ii) the overlap of subjects with the previous study is < 25% and the difference may reflect statistical noise; and iii) the previous study used DNA typing rather than the imputation used in the current study. It may be worth noting that DRB1*13 like DRB1*11 has an alanine at AA position 74 and thus contributes to the protective effect observed for this AA (p = 1.33 × 10 −11 , Supporting Table 3). Similarly, our previous study 6 showed only a marginal association of DRB1*14 (uncorrected p value = 0.004) compared with a strong association (p = 6.9 × 10 −7 ) observed in our current study. Here, the conditioning study results including the effect of controlling SNPs with very strong associations (see group 3, Table 4) provide additional support for the role of DRB1*14.
Finally, we have also considered specific HLA AAs. Most of the associated AAs are both nearly unique to the specific HLA classical alleles discussed above and also correspond to critical residues for the antigen binding pocket. Thus, associated AAs in DRB1 at AA positions 37, 47, 57, 60, 67, 70 and 74; and DPB1 at AA positions 9, 11, 76, 84, 87 are antigen pocket AAs 11,12 . Consistent with our results strong associations have been recently observed with serine at position 57 and leucine at position 74 in a Japanese PBC cohort 21 . We also note that many of the associated DQB1 AAs are also in critical residues for antigen binding (DQB1-AA13, 26, 70, 71, 74) 11,12 . Of the DRB1 associated only AA58 is not among the AAs in this functional class, whereas for DQB1 several are not in this functional class (DQB1-AA45, 56, 75, 167, 185).
In conclusion, the most parsimonious explanation consistent with the current study is that classical HLA genes and the coding variations within these genes are responsible for the HLA associations with PBC. Although we cannot exclude the possibility that other sequence variations affecting for example gene regulation could be important, our data indicates that a limited set of classical DRB1 and DPB1 alleles are sufficient to explain HLA associations with this disease. We believe the current data provides cogent information for understanding HLA-associations in PBC. Studies in other ethnic groups both within Europe and in other continental groups will also be important in further definition of the role of particular HLA genes and alleles. Lastly, our results provide additional rational for functional studies examining specific HLA genes and their relative binding to the putative disease associated epitopes of the PDC-E2 the immunodominant autoantigen epitopes of PBC 22 23 .

Study Population and design
The Italian PBC cases were obtained through a multi-center study and met internationally accepted criteria for the diagnosis of PBC as detailed in a previous study 6 . Each of the included cases also met ancestry criteria as defined below (see Ancestry). Controls were derived from several sources and this sample set information is detailed in Supporting Table  1. After data filtering and ancestry analyses contained 676 Italian PBC cases and 1440 Italian controls. All subjects enrolled in the study provided written informed consent and the study followed ethical guidelines of the most recent revision of the Declaration of Helsinki (Edinburgh, 2000).
All samples were genotyped with either Illumina (San Diego, CA 92121) genome-wide and/or Immunochip SNP platforms and the participants included the dataset from our previous GWAS as well as new samples (see Supporting Table 1). With the exception of ancestry information and assessment of relatedness the current study was restricted to genotypes in an ~4 megabase segment of human chromosome 6 (bps 28911802 -33813043, HG19 map). This dataset comprised a minimum of 1548 and a maximum of 5489 genotyped SNPs in each individual (Supporting Table 1) and was used for the SNP and HLA imputations (see Imputation).

Data Filtering
We used stringent quality control criteria to ensure that high-quality data were included in the analyses. We excluded individuals who had >5% missing data and all individuals with cryptic relatedness and duplicate samples based on identity-by-descent status for genomewide SNPs (PI^ > 0.15) using PLINK 24 .
We included only SNPs with <5% missing data, Hardy-Weinberg (H-W) equilibrium p values >10 −4 in controls and >10 −5 in combined cases and controls (to exclude most genotyping errors) applying these procedures in a stepwise approach separately for each dataset. For each of the separately derived control genotyping sets (Supporting Table 1), SNPs were excluded if they failed the above criteria within the individual control set or in combination with any of the other control groups, or in the complete data set. The H-W criteria were applied after exclusion of non-European individuals (see Ancestry). Finally, SNPs were excluded if allele frequency differed by > 10% between different control subject groups.

Ancestry
European ancestry was determined using 883 genome-wide SNPs with minimal or no linkage disequilibrium (LD) (r 2 < 0.1). SNPs analyzed using the STRUCTURE v2.1 program 25 and subjects of known European, Amerindian, East Asian and West African origin as previously described 26 . We used STRUCTURE to exclude non-European and admixed study participants since this method allows exclusion/inclusion criteria to be set using reference populations. Subjects with >15% non-European ancestry were excluded from further analysis.
Italian ancestry was defined using principal components analyses (PCA). For subjects with GWAS data, we used the same methods and criteria applied in a previous study with largely the same dataset. Briefly, PCA was performed using the EIGENSOFT statistical package 27 utilizing 34 thousand SNPs distributed throughout the genome (r 2 < 0.1) that we have previously used to define population genetic substructure 2 . These analyses used an independent set of Italian subjects for establishing membership [+/− 2 standard deviations (sd) in first 4 principal components (PCs)]. In the current study, a substantial portion of the samples did not have GWAS data (Supporting Table 1). For these samples we used a set of 12,579 SNPs from the Immunochip for which our empiric analyses demonstrated the ability of this set to discern Italian ancestry and exclude both other European ethnicities including Sardinian Italians (Supporting Fig. 1). Using the GWAS defined individuals, the Immunochip only genotyped samples were included using 2 sd in the first 4 PCs. In addition to the subject selection, we used the eigenvalues from the first four PCs (only the first four PCs were significant based on Tracy-Widom statistics) as covariates in our association analyses.

Imputation
We imputed SNPs, HLA classical alleles, and HLA gene amino acids (AAs) using phased reference genotypes from both the 1000 genome sequencing project (interim release June 2011) (http://www.1000genomes.org/) and an HLA defined reference set 8 . For the 1000 genome imputation we used IMPUTE version 2 28 under default parameters. The reference haplotypes for this imputation were from 1094 subjects including 381 European subjects and 98 Tuscan Italians. The number of genotyped (inference) SNPs that overlapped with the 1000 genome reference set ranged from 1435 SNPs (samples typed by GWAS), 4386 (samples typed by Immunochip) to 4981 SNPs (samples typed by GWAS plus Immunochip) (Supporting Table 1). For subsequent data analyses we utilized only imputed genotypes with maximum posterior probability scores of > 0.90. Using this parameter our empiric testing (leave one-out analyses) indicated that the maximum error rate for genotype assignment was < 0.05 and the mean error rate was < 0.01. To impute classical HLA alleles and corresponding amino acids determinants we utilized a reference separate dataset of collected by the Type 1 Diabetes Genetics Consortium (T1DGC). This reference data contains genotype data for 2,537 SNPs, selected to tag the entire MHC, and classical types for HLA-A, B, C, DRB1, DQA1, DQB1, DPA1 and DPB1 at 4-digit resolution in 2767 unrelated individuals of European descent 29 . The Beagle software package 30 was used for this imputation under default parameters. The number of inference SNPs that overlapped with this reference dataset ranged from 648 (samples typed by GWAS), 1444 SNPs (samples typed by Immunochip) to 1610 SNPs (GWAS plus Immunochip) (Supporting Table 1). Similar to the imputation using 1000 Genome data, only SNPs with posterior probabilities of >0.90 were included in our final analyses. For imputed SNPs that overlapped between the two imputation sets and algorithms used (Impute V2.0 and Beagle) there was a nearly complete concordance of the association testing results indicating similar performance of these algorithms for this dataset. After imputation and selecting only those markers meeting posterior probability criterion this region contained a total of 49,885 markers including the genotyped SNPs that were included in association test analyses.

Association and Conditional Association Tests of Imputed SNPs and HLA determinants
SNPTEST V2.0 28 (web) was used for the primary association analyses for the imputed genotypes. This software uses the genotype probabilities for the imputed SNPs or determinants and accounts for genotype uncertainty. The first four PC eigenvalue scores were used as continuous variables in the association test together with the gender covariate. Analyses were performed using the SNPTEST v2 Score test algorithm that enabled both the inclusion of the covariates and conditioning tests and all of our reported results used an additive model. To minimize potential spurious results we limited our main and conditioning analyses to markers with information scores (Inf) > 0.85. This parameter is a measure of the observed statistical information for the estimate of SNP allele frequency (for additional information see https://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.v2.pdf).
Conditioning on multiple markers either separately or together was performed using an additive model. For the HLA region, over 150 conditional analyses were performed using the SNPs and HLA determinants including all HLA determinants with p values < 10 −6 .
Nominal p values after correction for covariates and conditioning are provided throughout the manuscript. The p-values <10 −6 would remain significant after conservative (Bonferroni) correction for the number of markers (<50,000) tested after imputation.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.   all classical HLA determinants with p-values < 10-6 (rows 9-16); and 4) the top four HLA amino acid (AA) associations (rows [17][18][19][20]. P values after conditioning on particular HLA genes (see methods) are also shown. The p values were determined using SNPTEST V2 (see Patients and Methods) and used the genotype probabilities. A complete summary of strongly association markers is provided in Supplementary Table 1. b The SNP name, classical HLA determinant or AA is shown. c The position corresponding to the marker is shown (chromosomal base pair positions from HG19).

Invernizzi et al. Page 17
d Allele 1 is the minor allele used for the analysis. For sequences designated with amino acid (AA) numbers the symbol refers to the allele symbol corresponds to standard AA nomenclature. Where the AA is compared with multiple other AAs, "Other" is used to represent all alternative AAs. The AA position includes the leader sequence (e.g. DRB1* position 104 is the same as AA position 74 without the leader sequence). Where HLA determinants are provided, the specificity is compared with not having the specificity (OTHER).
f p values before conditioning using population substructure and gender as covariates. g p value after conditioning on the allele(s) listed in the heading line.