Importance of HBsAg recognition by HLA molecules as revealed by responsiveness to different hepatitis B vaccines

Hepatitis B (HB) vaccines (Heptavax-II and Bimmugen) designed based on HBV genotypes A and C are mainly used for vaccination against HB in Japan. To determine whether there are differences in the genetic background associated with vaccine responsiveness, genome-wide association studies were performed on 555 Heptavax-II and 1193 Bimmugen recipients. Further HLA imputation and detailed analysis of the association with HLA genes showed that two haplotypes, DRB1*13:02-DQB1*06:04 and DRB1*04:05-DQB1*04:01, were significantly associated in comparison with high-responders (HBsAb > 100 mIU/mL) for the two HB vaccines. In particular, HLA-DRB1*13:02-DQB1*06:04 haplotype is of great interest in the sense that it could only be detected by direct analysis of the high-responders in vaccination with Heptavax-II or Bimmugen. Compared with healthy controls, DRB1*13:02-DQB1*06:04 was significantly less frequent in high-responders when vaccinated with Heptavax-II, indicating that high antibody titers were less likely to be obtained with Heptavax-II. As Bimmugen and Heptavax-II tended to have high and low vaccine responses to DRB1*13:02, 15 residues were found in the Heptavax-II-derived antigenic peptide predicted to have the most unstable HLA-peptide binding. Further functional analysis of selected hepatitis B patients with HLA haplotypes identified in this study is expected to lead to an understanding of the mechanisms underlying liver disease.

continued until 2016 in Japan, the Japanese government started routine HBV vaccination for infants from October 2016. Several hepatitis B (HB) vaccines have been developed to prevent HBV infection corresponding to HBV genotypes [i.e. genotype A (Heptavax-II R ), genotype A2 (Engerix-B R ; Recombivax HB R ), or genotype C (Bimmugen R )]. In Japan, Bimmugen derived from HBV genotype C has been used in accordance with the background of hepatitis B patients, but in recent years, Heptavax-II derived from HBV genotype A has been increasingly used. Both HB vaccines are given three times, and antibody titers > 100 mIU/mL are considered to provide adequate protection against HBV infection over two years 3,4 . However, the immune response to HB vaccination differs among individuals, with 5-10% of healthy individuals failing to acquire protective levels of antibodies.
There have been several reports in which associations of SNPs in the HLA class II region with a response to HB vaccines have been identified, mainly in Asian populations [5][6][7][8] . In these previous reports, the studied individuals were vaccinated with HB vaccines designed based on the HBV genotype A2. In our previous report, a genomewide association study (GWAS) using a total of 1193 Japanese individuals who were vaccinated with Bimmugen revealed the importance of the HLA-DR-DQ and BTNL2 genes for response to an HB vaccine designed based on HBV genotype C 9 . To date, no studies have examined the response of different types of HB vaccines in the same population. Here, we conducted a genomic analysis, i.e. SNP-based GWAS and association tests of HLA class II alleles, of Japanese individuals inoculated with Heptavax-II vaccine and investigated the differences between the genetic factors related to Heptavax-II responsiveness and those related to Bimmugen responsiveness.

Results
GWAS revealed that HLA genes significantly contribute to the responsiveness to both Heptavax-II and Bimmugen. Genome-wide multiple regression analysis was carried out to compare poorresponders (i.e. Group_0) and responders (i.e. Group_1 + Group_2) ( Supplementary Fig. 1A), and poor responders and high-responders (i.e. Group_2) ( Supplementary Fig. 1B), in which age, sex, and the number of vaccinations were used as covariates. No significant association was detected in genome-wide association analysis after genotype imputation (Fig. 1). The SNP rs4248166 present on the BTNL2 gene, which showed a significant association in the Bimmugen GWAS 9 , was not significantly associated in the Heptavax GWAS (P = 0.018). When the OR in a GWAS (Group_0 vs. Group_1 + Group_2) of Heptavax-II was compared to the one in Bimmugen, the top 70 SNPs with P values lower than 0.001 in GWAS of Bimmugen were highly correlated (r 2 = 0.886) (Fig. 2). Sixty-six SNPs out of 70 were located in the HLA region (Chr.6: 32,570,311-33,080,360, dbSNP b151), and only one SNP, rs12595417, which is located in an intron of long non-coding RNA (LINC02250, Chr. 15: 25,793,489), showed an opposite OR in the GWAS between Heptavax-II (OR = 0.52) and Bimmugen vaccines (OR = 1.87) (Supplementary Table 1). A combined association analysis of the cases of Heptavax-II and Bimmugen was subsequently performed, because the genetic background for their responsiveness was found to be similar. Significant associations were detected from the HLA region in genome-wide association analysis after genotype imputation in both comparisons between poor-responders and responders (Supplementary Fig. 2A), and between poor-responders and high-responders ( Supplementary Fig. 2B). Finally, we carried out genomewide association analysis after genotype imputation to compare each of the three groups between Heptavax-II and Bimmugen. Although a significant association of rs57190718, located 50 kb downstream from the PCDH10 gene, was detected when high-responders were compared between Heptavax-II and Bimmugen ( Supplementary  Fig. 3C), no significant association was found in the validation test by TaqMan assay. The other three SNPs (rs4287780, Chr.2: 166860126 on GRCh38; rs10904760, Chr.10: 16507427 on GRCh38; rs17770803, Chr.15: 24694843 on GRCh38) that met the genome-wide significance level were also genotyped by TaqMan assay, but the associations were not replicated and were therefore revealed to be false positives.
Two HLA-DRB1-DQB1 haplotypes showed different responses to Heptavax-II and Bimmugen. Since the GWAS suggested that HLA genes play a very important role in responsiveness to either Heptavax-II or Bimmugen vaccines, HLA imputation was performed, and HLA allele frequencies were compared in detail. We performed statistical imputation of classical HLA alleles for three HLA loci including HLA-DRB1, DQB1, and DPB1 using 555 genome-wide SNP typing data as established in our previous report 10,11 . A total of 515 individuals, all of whom showed high-quality callings in imputation (CT > 0.5) for three HLA alleles, were used for the association tests in a comparison between poor-responders and high-responders (Table 1), and between poor-responders and responders (Supplementary Table 2). When poor-responders were compared with high-responders, i.e. excluding intermediate-responders, both HLA-DRB1*04:05-DQB1*04:01 haplotype (OR = 2.09, P = 5.84 × 10 -4 ) and DPB1*05:01 allele (OR = 2.03, P = 3.70 × 10 -4 ) were found to have higher ORs than the same haplotype and allele found in a comparison between poor-responders and responders, and these met the significance level. Subsequently, ORs from HLA association tests (i.e. poor responders vs. high responders) for Heptavax-II were compared with those for Bimmugen. All 13 HLA class II alleles showing P < 0.05 in the association tests of Bimmugen showed the same OR trend in both Heptavax-II and Bimmugen, and were highly correlated (r 2 = 0.842) (Supplementary Table 3). The frequencies of HLA-DRB1-DQB1 haplotypes and DPB1 alleles were compared between Heptavax-II and Bimmugen recipients in high-responders (Table 2), responders (Supplementary Table 4), and poor-responders (Supplementary Table 5 Table 6). In contrast, no DRB1-DQB1 haplotype showed a significant association in Bimmugen recipients.  Table 7). Binding coefficients were estimated for a total of 212 peptides for both Heptavax-II and Bimmugen, 120 of which had different binding coefficients (IC 50 ) between Heptavax-II and Bimmugen. The predicted bindings to DRB1*04:05 and DRB1*13:02 were performed using the NetMHCII 1.1 method, and the top 10 peptides with small IC 50 (that is, HLA-peptide binding is strong) are summarized in Supplementary  Tables 8 and 9, respectively. The frequency of the DRB1*13:02 allele was significantly lower in the group with

Discussion
In SNP-based association analysis, no SNP significantly associated with a response to the Heptavax-II vaccine. However, the genetic factors involved in the responsiveness to both Heptavax-II and Bimmugen vaccines were found to be very similar. Seventy SNPs with P values less than 0.0001 in the GWAS of 1,193 Bimmugen recipients were highly correlated, with r 2 = 0.8856, when compared with the odds ratios obtained in the GWAS of 555 Heptavax-II recipients. Of the 70 SNPs, 66 were located in the HLA class II region, and only 1 (rs12595417 located on the LINC02250 gene) had the OR reversed. Although rs12595417 did not reach a genome-wide significance level in either Heptavax-II or Bimmugen GWAS (P = 1.82 × 10 −5 and P = 1.87 × 10 −3 , respectively), it is present in a multigenic region containing the SNHG14 gene, which has been reported to be associated with the progression of cervical 12 and lung cancer 13 , UBE3A, which is known to be critical for the development of human papillomavirus-associated cancers 14 , and ATP10A (also known as ATP10C), which has been reported to encode www.nature.com/scientificreports/ a putative aminophospholipid translocase associated with Angelman syndrome 15 . Because the ORs of SNPs with low P values in the GWAS of Heptavax-II and Bimmugen were highly correlated, the combined GWAS including both Heptavax-II and Bimmugen and comparisons of each of the three groups (poor-responders, responders, and high-responders) between Heptavax-II and Bimmugen was performed. The combined GWAS showed that HLA class II genes contribute significantly to hepatitis B vaccine responsiveness. In order to understand the contribution of HLA class II genes to vaccine responsiveness in more detail, we carried out HLA imputation to determine HLA alleles, and performed HLA association tests. As a result of HLA association tests, HLA-DRB1*04:05-DQB1*04:01 haplotype and DPB1*05:01 allele showed significant association with Heptavax-II responsiveness. These HLA class II genes were the same as those that we previously reported showing a significant association in Bimmugen responsiveness 9 . Therefore, we compared the ORs of HLA class II alleles that were significant in HLA association analysis of Bimmugen, and found a very high correlation between Heptavax-II and Bimmugen (r 2 = 0.8419). In addition, HLA allele/haplotype frequencies were compared in each of the three groups. Unexpectedly, two DRB1-DQB1 haplotypes, DRB1*04:05-DQB1*04:01 and DRB1*13:02-DQB1*06:04, were significantly associated in a comparison of high-responders (P = 2.22 × 10 -3 , P = 3.40 × 10 -4 , respectively).
HLA-DRB1*04:05-DQB1*04:01 was significantly associated with vaccine hypo-responsiveness in both Heptavax-II and Bimmugen recipients (OR = 2.09 and OR = 2.85, respectively). From the greater OR for Bimmugen, it can be understood that individuals with HLA-DRB1*04:05-DQB1*04:01 are more likely to have a significantly lower vaccine response with Bimmugen than with Heptavax-II. As a result of the binding prediction of HBsAg peptide to DRB1*04:05 by IEDB, there were many peptides predicted only for Heptavax-II that bind more stably than those predicted only for Bimmugen (i.e. low IC 50 ). It can be understood that due to the presence of many HBsAg peptides derived from Heptavax-II that bind stably to DRB1*04:05, DRB1*04:05 is an allele that contributes to vaccine hypo-responsiveness but to a lesser extent than to Bimmugen.
HLA-DRB1*13:02-DQB1*06:04 was a haplotype that was not significantly associated with vaccine responsiveness in either Heptavax-II or Bimmugen recipients (P = 2.91 × 10 -1 , P = 1.09 × 10 -1 , respectively). However, when high-responders in both vaccinations were compared, the frequency of the HLA-DRB1*13:02-DQB1*06:04 haplotype was significantly lower in the Heptavax-II group. Although the difference was not significant, focusing on the ORs in Heptavax-II-vaccinated and Bimmugen-vaccinated individuals, the haplotype was associated with a tendency for low vaccine responses in Heptavax-II (OR = 1.66) and high vaccine responses in Bimmugen (OR = 0.53). From this result, it can be understood that individuals with the HLA-DRB1*13:02-DQB1*06:04 haplotype can obtain higher antibody titers (HBsAb > 100 mIU/mL) by being vaccinated with Bimmugen rather than with Heptavax-II. In the binding prediction of vaccine-derived antigen peptides to the HLA-DRB1*13:02 allele by IEDB, a peptide that is unstable for Heptavax-II but is stable for Bimmugen was found. The 15 residues from position 134-148 (FPSCCCTKP [T/S] DGNCT) are the most suitable for binding.
A previous report in Japan has examined whether Bimmugen recipients are protected against HBV genotype A infection and whether Heptavax-II recipients are protected against HBV genotype C infection 16 . The results showed that sufficient antibody titers (i.e. HBsAb > 100 mIU/mL) could protect against infection from different HBV genotypes. Our findings showed that HBsAb > 100 mIU/mL is less likely to be obtained in individuals with the HLA-DRB1*13:02-DQB1*06:04 haplotype when vaccinated with Heptavax-II, indicating that Bimmugen vaccination is preferred for these individuals.
As we reported previously, HLA-DRB1*13:02-DQB1*06:04 is a significantly less frequent haplotype in patients with chronic hepatitis B 17 (Table 3). Taken together with the reverse contribution of Heptavax-II and Bimmugen to vaccine responsiveness, it can be understood that HBsAg peptides common to Heptavax-II and Bimmugen do not contribute to disease susceptibility. Approximately 80% of hepatitis patients in Japan are infected with HBV genotype C, but the mechanism by which HLA molecules contribute to disease resistance remains unclear. If HLA-DRB1*13:02-DQB1*06:04 is significantly associated with a high vaccine response in Bimmugen recipients from the same HBV genotype C, it may be demonstrated that recognition of Bimmugen-derived HBsAg peptide by HLA-DRB1*13:02-DQB1*06:04 plays an important role in disease resistance.
In conclusion, genomic analyses of a Japanese population vaccinated with different hepatitis vaccines (Heptavax-II and Bimmugen) revealed that there are specific HLA-DRB1-DQB1 haplotypes with different responses, although the genetic factors involved in vaccine responsiveness are very similar. Further functional analysis of selected hepatitis B patients with HLA haplotypes identified in this study is expected to lead to an understanding of the mechanisms underlying liver disease. Table 3. Haplotype frequencies and genotype counts in Japanese vaccines (high-responders), Japanese chronic hepatitis B patients, and healthy Japanese individuals. *Significant association was identified in comparison with healthy Japanese individuals (OR = 0.34, P = 4.30 × 10 −4 ). **Significant association was identified in comparison with healthy Japanese individuals (OR = 0.42, P = 2.05 × 10 −10 ).

Materials and methods
Ethics approval. This study was approved by the Ethics Committee of the National Center for Global Health and Medicine and of all of the following Institutes and Hospitals throughout Japan that participated in this collaborative study: the University of Tokyo; Kawasaki Medical School; National Hospital Organization Nagasaki Medical Center; Osaka City University; Saga University; the University of Tsukuba; Iwate Medical University; and Chiba University. All participants provided written informed consent for participation in this study, and the methods were carried out in accordance with the approved guidelines.
Samples and clinical data. New cases vaccinated with a recombinant absorbed HB vaccine (Heptavax-II R , MSD KK, Tokyo, Japan) were collected under conditions similar to those used in our previous study using individuals vaccinated with the Bimmugen R vaccine (Kaketsuken, Kumamoto, Japan) 9 . All 555 Japanese genomic DNA samples were obtained from healthy adult volunteers (≥ 18 years) who were vaccinated in three doses (0.5 mL) at 0, 1, and 6 months with the Heptavax-II vaccine. Only those who received all three doses of the Heptavax-II vaccine were included in this study. Serum anti-HBV surface antibody (HBsAb) and serum anti-HBV core antibody (HBcAb) were tested before vaccination and at 1 month after final inoculation using the ARCHI-TECT anti-HBs kit (Abbott Japan, Tokyo, Japan) and the ARCHITECT anti-HBc II kit (Abbott Japan, Tokyo, Japan), respectively, with a fully automated chemiluminescent enzyme immunoassay system using the Architect i2000SR analyzer (Abbott Japan, Tokyo, Japan). Individuals who were HBcAb-positive (> 1.0 S/CO) were not included in this study. Although there has been a report of an HBV vaccinated donor showing an HBsAb value of 96 IU/L who was found to be positive for HBV DNA, we set levels of antibodies against HBV surface antigen as 10 IU/L or more to be considered as having immunity against HBV, following the World Health Organization recommendation 4  Genome-wide SNP genotyping in individuals vaccinated with Heptavax-II. For the GWAS, the 555 Japanese genomic DNA samples from individuals vaccinated with Heptavax-II were genotyped using the Axiom Genome-Wide ASI 1 array, according to the manufacturer's instructions. The following analysis was performed according to the same way as the analysis in our previous GWAS on Bimmugen 9 . All samples had an overall call rate of more than 95%; the average overall call rate was 99.41% (95.01-99.93) and passed a heterozygosity check. No related individual (PI ≥ 0.1) was identified in identity-by-descent testing. Principal component analysis was carried out to check the genetic background in the studied 555 samples together with HapMap samples (43 JPT, 40 CHB, 91 YRI, and 91 CEU samples). This analysis showed that all 555 samples formed the same cluster using the first and second components, indicating that the effect of population stratification was negligible in the studied samples ( Supplementary Fig. 4).

SNP genotype imputation and association test.
For high-density mapping following the GWAS, we carried out SNP genotype imputation from 555 Heptavax-II and 1193 Bimmugen recipients by Beagle version 4.1 (URL: https ://facul ty.washi ngton .edu/brown ing/beagl e/b4_1.html) using the 1000 genomes as the imputation reference 18 . The following thresholds were then applied for SNP quality control during the data cleaning: SNP call rate ≥ 95%; minor allele frequency ≥ 5%; and Hardy-Weinberg Equilibrium P value ≥ 0.001. Of the SNPs on autosomal chromosomes, 5,835,420 SNPs finally passed the quality control filters and were used for genomewide multiple regression analysis using age, sex, and the number of vaccinations as covariates. Figure 1 shows association analysis using 555 Japanese HB vaccinated individuals with Heptavax-II in a comparison between poor-responders (Group_0) and responders (Group_1 + Group_2) (Fig. 1A), and between poor-responders and high-responders (Group_2) (Fig. 1B). A combined analysis of Heptavax-II (n = 555) and Bimmugen (n = 1193) ( Supplementary Fig. 2), and a comparison of each group (i.e. Group_0, Group_1, and Group_2) to two different vaccines ( Supplementary Fig. 3), were carried out. To avoid false positives due to multiple testing, genome-wide significance was set at P < 5.0 × 10 −8 . Four SNPs satisfying the significance level in a comparison of Group_2 between Heptavax-II and Bimmugen were validated to see if the association could be reproduced by the TaqMan assay (Thermo Fisher Scientific Inc., MA, USA).
HLA imputation. SNP data from 555 samples were extracted from an extended MHC (xMHC) region ranging from 25,759,242 to 33,534,827 bp based on the hg19 position. We conducted 2-field HLA genotype imputation for three class II HLA genes using the HIBAG R package as the same way in our previous report 9 . For HLA-DRB1, DQB1, and DPB1, our in-house Japanese imputation reference was used for HLA genotype imputation 10,11 . We applied post-imputation quality control using call-threshold (CT > 0.5). A total of 515 samples consisting of 63 in Group_0, 174 in Group_1, and 278 in Group_2 showed three estimated HLA genotypes, Haplotype estimation. Haplotype phasing was performed in a manner similar to that performed in the previous analysis of Bimmugen 9 . The phased haplotypes consisting of three HLA class II loci (HLA-DRB1, DQB1, and DPB1) were estimated using the PHASE program version 2.1 (URL: http://steph ensla b.uchic ago. edu/index .html) 19 . The estimated 3-locus haplotypes were further used for the estimation of haplotypes of HLA-DRB1 and DQB1 loci (i.e., the collapsing method was applied to the phased data for three HLA loci). See our previous paper for more details 9 .
HLA association test. To assess the association of HLA allele or haplotype with a response to HB vaccination, P values and odds ratios (OR) were calculated by Pearson's chi-square test in the presence vs. absence of each allele or haplotype. In Table 1 and Supplementary Table 2, HLA association analysis using 555 Japanese HB vaccinated individuals with Heptavax-II was carried out in a comparison between poor-responders (Group_0) and high-responders (Group_2), and between poor-responders (Group_0) and responders (Group_1 + Group_2). HLA association analysis between Heptavax-II and Bimmugen recipients was carried out in high-responders, responders, and poor-responders ( Table 2, Supplementary Tables 4 and 5, respectively). HLA frequencies of DRB1-DQB1 haplotypes were compared between healthy individuals and high-responders in Heptavax-II and Bimmugen inoculation (Supplementary Table 6). The healthy individuals used here were the same as those used in our previous study 17 .
To avoid false positives due to multiple tests, the family-wise error rate was set to be less than 0.05 in the association analysis. The threshold of the P value for each allele frequency table or haplotype frequency table was determined by the permutation test.
Web resources. GWAS data in this study will be submitted to the NBDC Human Database, a public database in Japan (NBDC Human Database, https ://human dbs.biosc ience dbc.jp/en/).