Front. Immunol.Frontiers in ImmunologyFront. Immunol.1664-3224Frontiers Media S.A.10.3389/fimmu.2021.664024ImmunologyOriginal ResearchIdentification of Two Novel Candidate Genetic Variants Associated With the Responsiveness to Influenza VaccinationWenSimin1WeiHejiang2LiaoQijun1LiMao1ZhongShuyi1ChengYanhui2HuangWeijuan2WangDayan2*ShuYuelong1*1School of Public Health (Shenzhen), Sun Yat-sen University, Guangzhou, China2National Institute for Viral Disease Control and Prevention, Chinese Center for Disease Prevention and Control, Beijing, China
Edited by: Rong Hai, University of California, Riverside, United States
Reviewed by: Arun Iyer, Boehringer Ingelheim, United States; Duo Xu, University of California, Riverside, United States
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Background
Annual vaccination is the most effective prevention of influenza infection. Up to now, a series of studies have demonstrated the role of genetic variants in regulating the antibody response to influenza vaccine. However, among the Chinese population, the relationship between genetic factors and the responsiveness to influenza vaccination has not been clarified through genome-wide association study (GWAS).
Method
A total of 1,968 healthy volunteers of Chinese descent were recruited and 1,582 of them were available for the subsequent two-stage analysis. In the discovery stage, according to our inclusion criteria, 123 of 1,582 subjects were selected as group 1 and received whole-genome sequencing to identify potential variants and genes. In the verification stage, 29 candidate variants identified by GWAS were selected for further validation in 481 subjects in group 2. Besides, we also analyzed nine variants from previously published reports in our study.
Results
Multivariate logistic regression analysis showed that compared with the TT genotype of ZBTB46 rs2281929, the TC + CC genotype was associated with a lower risk of low responsiveness to influenza vaccination adjusted for gender and age (Group 2: P = 7.75E-05, OR = 0.466, 95%CI = 0.319–0.680; Combined group: P = 1.18E-06, OR = 0.423, 95%CI = 0.299–0.599). In the combined group, IQGAP2 rs2455230 GC + CC genotype was correlated with a lower risk of low responsiveness to influenza vaccination compared with the GG genotype (P = 8.90E-04, OR = 0.535, 95%CI = 0.370–0.774), but the difference was not statistically significant in group 2 (P = 0.008). The antibody fold rises of subjects with ZBTB46 rs2281929 TT genotype against H1N1, H3N2,and B were all significantly lower than that of subjects with TC + CC genotype (P < 0.001). Compared with IQGAP2 rs2455230 GC + CC carriers, GG carriers had lower antibody fold rises to H1N1 (P = 0.001) and B (P = 0.032). The GG genotype of rs2455230 tended to be correlated with lower antibody fold rises (P = 0.096) against H3N2, but the difference was not statistically significant. No correlation was found between nine SNPs from previously published reports and the serological response to influenza vaccine in our study.
Conclusion
Our study identified two novel candidate missense variants, ZBTB46 rs2281929 and IQGAP2 rs2455230, were associated with the immune response to influenza vaccination among the Chinese population. Identifying these variants will provide more evidence for future research and improve the individualized influenza vaccination program.
influenzavaccineimmune responsegenetic variantgenome-wide association studyZBTB46IQGAP2Shenzhen Science and Technology Innovation Program10.13039/501100017610National Natural Science Foundation of China10.13039/501100001809Introduction
Influenza is an acute infectious respiratory disease with an annual occurrence of one billion cases worldwide, including three to five million serious cases and 290,000–650,000 deaths (1). For some vulnerable populations such as infants, pregnant women, and elderly people, influenza poses a serious threat to their lives and health. Annual vaccination is the most effective prevention of influenza infection. However, according to the data from the World Health Organization (WHO) (2), the average effectiveness of seasonal influenza vaccine over the last 16 years was only 39.9%, suggesting that current vaccines cannot provide sufficient protection. There are at least two factors that may influence vaccine effectiveness (VE). One is the matching degree of circulating strains and vaccine strains, and the other is the immunogenicity of vaccine. Currently, although a well-developed global influenza surveillance network has been established, the prediction of annual vaccine strains remains a huge challenge. Therefore, it is of great significance to improve the immunogenicity of vaccines in case the vaccine strains do not match circulating strains. Many factors may affect the immune response to influenza vaccination, including gender, age, adjuvant use, delivery mode, and so on. For individuals who share the same features mentioned above, their immune responses to the same vaccine may also vary significantly. This indicates that genetic factors may play a role in the immunologic processes.
Single nucleotide polymorphism (SNP) is one of the most common types of heritable variation in humans. Currently, some studies have observed some SNPs that may affect the responsiveness to influenza vaccination. HLA is a polygenic and polymorphic complex involved in antigen presentation. The relationship between SNPs in HLA and the immune responses to influenza vaccination has been demonstrated (3–7). Besides, some studies also reported that genetic variants in some immune-related genes may also influence the immunological responses (6, 8).
To our knowledge, a genome-wide association study (GWAS) has not yet been performed to examine the genetic factors that were related to the responsiveness to influenza vaccination in the population of Chinese descent. Therefore, we conducted a GWAS in healthy volunteers of Han Chinese. A replication study was also conducted to validate the results of the GWAS. In addition, in our cohort, we also validated some SNPs previously reported to be associated with the immune responses to influenza vaccine. Identification of these genetic variants may provide more evidence for future research and improve the individualized influenza vaccination program.
Materials and MethodsSubjects and Study Design
A total of 1,968 healthy volunteers from September 2009 to September 2019 were enrolled from Yunnan Center for Disease Control and Prevention and Urumqi Center for Disease Control and Prevention. All participants were administered intramuscular injections of trivalent inactivated vaccine (TIV) which contains 15 μg hemagglutination (HA) before the flu season. The vaccine strains of TIVs used in our study were consistent with the northern hemisphere vaccine component recommended by WHO. Blood samples were collected before and 28 days after the vaccination. The exclusion criteria were as follows: not Han Chinese, lost to follow-up, inadequate blood samples, and repeated vaccination. In the end, 386 subjects were excluded, and 1,582 subjects in total were suitable for further research.
According to the age of subjects, we divided them into four groups: infants (≤5 years), children (6–17 years), adults (18–64 years) and elderly people (≥65 years). Seroprotection rate was defined as a percentage of subjects with post-vaccination titer ≥1:40. Seroconversion rate was defined as a percentage of subjects with either a pre-vaccination hemagglutination inhibition (HAI) titer <1:10 and post-vaccination HAI titer ≥1:40 or a pre-vaccination HAI titer ≥1:10, and a minimum four-fold increase in post-vaccination HAI titer. Considering the results of the HAI test, we selected some of the subjects according to our inclusion criteria and divided them into two groups. In group 1, low responders (LRs) were subjects from phenotypic extremes whose HAI titer reached neither the seroprotection rate nor the seroconversion rate to all three vaccine strains. Under this strict standard, only 41 of 1,582 subjects were eligible to be classified as LRs in group 1. According to the age and gender of these 41 LRs, 82 matched subjects were selected as responders from the cohort, the matching ratio was 1:2. Responders were defined as those whose HAI titers reached both the seroprotection rate and the seroconversion rate to all three vaccine strains. In this way, a total of 123 healthy volunteers with 41 LRs and 82 responders were included in group 1. In group 2, subjects included in group 1 were all removed. To validate the potential variants in larger sample size, we broadened the inclusion criteria appropriately. LRs were defined as HAI titer failed to reach the seroconversion rate to all three vaccine strains. And responders were defined as HAI titer reached the seroconversion rate to all three vaccine strains. Therefore, 481 subjects with 193 LRs and 288 responders were included in group 2. The combined group is a combination of group 1 and group 2. All the individuals in group 1 met the inclusion requirements of group 2.
In this study, we conducted two-stage sequencing. In the discovery stage, 123 subjects in group 1 received whole-genome sequencing (WGS) to identify potential variants and genes through GWAS and gene-based association study. In the verification stage, 29 potential variants from GWAS were selected for genotyping to validate the results in group 2. The flow chart of this study is shown in Figure 1.
Study flow chart.
Whole-Genome Sequencing and Quality Control
Genomic DNA was extracted from peripheral blood leukocytes using a TIANamp Blood DNA kit (TIANGEN BIOTECH (Beijing) Co., LTD, China) according to the manufacturer’s instructions. We performed WGS on the BGISEQ-500 sequencing platform with paired-end 100 base run. The average sequencing depth was 39.09×, and the coverage rate was 96.49% (Supplementary Table 1). Raw reads were filtered by removing adapters and low-quality sequences; the remaining clean reads of each sample were aligned to hg38 using the Burrows–Wheeler Aligner (BWA) software (version 0.75) (9). Duplicated reads were marked using Picard (https://broadinstitute.github.io/picard). Base quality score recalibration was then performed using the Genome Analysis Toolkit (GATK v3.7) (10). SNPs and small insertions/deletions (indels) were detected by GATK HaplotypeCaller and GenotypeGVCFS tools according to the best practices recommended. To obtain high-quality variants, we applied a series of exclusion criteria to remove variants that had an average depth <8, a minor allele mean depth <4, <90% samples covered at least eight reads, a low mapping quality score, strand bias or allelic imbalance, and deviations from Hardy–Weinberg equilibrium (HWE) (Supplementary Table 2). Retained variants were annotated using ANNOVAR (11).
For sample quality control, samples with a mean sequencing depth <30, <90% genome covered by at least 10×, inbreeding coefficient >0.1 or <−0.1, a 2nd-degree or closer relationship with other samples or a population outlier (smartpca) (12) were excluded. Therefore, a total of 12 samples including seven LRs and five responders were removed according to these exclusion criteria.
SNP Selection and Genotyping
According to the results of GWAS, six highly significant variants (P <1 × 10−5) were selected for further validation. To identify more potential variants, we also focus on SNPs (P <1 × 10−2) in coding region. Considering the function of genes, 23 exonic SNPs of immune-related genes were selected. Therefore, a total of 29 variants were included in the subsequent verification analysis. In addition, we also validate nine previously reported SNPs associated with the responsiveness to influenza vaccination. And the minor allele frequencies (MAFs) of these nine SNPs were no less than 0.05 in the Chinese Han population according to the 1000 Genomes Project data. Genotyping of each SNP was conducted using the MassARRAY technology platform (Sequenom, San Diego, California, USA) and determined by BioMiao Biological Technology (Beijing, China). Genotyping was performed by technicians who were blinded to the study design. According to the number of mutations carried by the subjects, their genotypes were classified into wild-type homozygote, mutant heterozygote, and mutant homozygote. SNPs with HWE-P <0.001 or call rate <95% were excluded.
HAI Test
Serological specimens from all subjects were separated and stored at −30°C. Red blood cells (RBCs) of turkey were used for influenza A/H1N1 and B HAI assays, while RBCs of guinea pig were used for influenza A/H3N2 HAI assays. Before the detection, receptor destroying enzyme (RDE) was mixed with the serum in the ratio of 3:1 and placed in a 37°C water bath for 16–18 h to remove the non-specific inhibitor. Then, the mixture was bathed in water at 56°C for 30 min to inactivate the RDE. The HAI titer was designated to be the highest serum dilution to inhibit hemagglutination. HAI assays were performed against the northern hemisphere influenza vaccine strains according to the standardized protocol by Technical Guidelines for National Influenza Surveillance (13). The recommended vaccine strains from 2009 to 2019 were summarized in Supplementary Table 3.
Prediction of SNP Function
The prediction of SNP function was performed using SNPinfo (https://snpinfo.niehs.nih.gov/) provided by the National Institute of Environmental Health Sciences (NIEHS).
Genetic Models
In this study, a total of four genetic models (additive genetic model, dominant genetic model, recessive genetic model, and overdominant genetic model) were applied in the analysis. The additive genetic model was used to compare the distribution frequencies of the three genotypes. The other three genetic models divided the subjects into two groups by different combinations of the three genotypes. The dominant genetic model compared the distribution frequencies of wild-type homozygotes with other subjects. The recessive genetic model compared the distribution frequencies of mutant homozygotes with other subjects. The overdominant genetic model divided the subjects into homozygotes and heterozygotes to analyze whether the mutation had heterozygous superiority. According to the results of data analysis, only the optimum genetic model with the highest impact was selected and presented.
Statistical AnalysisDiscovery Stage
For all the variants identified in WGS, we used Fisher’s exact test to assess the distribution frequencies of alleles and genotypes respectively. For variants with MAF >5%, we applied logistic regression adjusted for the top four principal components and age. Association analysis was performed by the PLINK (14) (v1.07). A gene-based association analysis was conducted for low-frequency variants with MAF <0.05. We defined two variant groups, including: 1. a “Deleterious (Broad)” set defined as non-sense, splice-site, indel frameshift, and missense, which is annotated as deleterious by at least one of the five protein prediction algorithms of likelihood ratio test (LRT) score, Mutation Taster, PolyPhen-2 HumDiv, PolyPhen-2 HumVar, and SIFT. 2. “Deleterious (Strict)” set comprising non-sense, splice-site, indel frameshift, and missense, which is annotated as deleterious by all five protein prediction algorithms. For variants that have different annotations from multiple transcripts of the gene, each variant with the highest impact was presented. We performed Fisher’s exact test and SKAT (15) (EPACTS v3.2.6, https://genome.sph.umich.edu/wiki/EPACTS) to calculate the gene-based collapsing on these two defined variants sets.
Verification Stage
The HWE test for assessing the SNP genotype frequency among subjects was conducted. Assessment of pairwise linkage disequilibrium (LD) was performed by the Haploview V.4.2 software. Categorical data were expressed in frequencies (percentages) and compared by χ2 test or Fisher’s exact test when appropriate. Associations between variants and responsiveness to influenza vaccination were calculated by multivariate logistic regression models adjusted for age and gender. When the fold rises of the antibodies are calculated, the HAI titers before and after vaccination were transformed to their logarithms. The antibody fold rises with unnormal distribution through Kolmogorov–Smirnov Test were described as the median [interquartile range (IQR)] and compared by Mann–Whitney U Test. The level of statistical significance was P <0.05. In the analysis of the relationship between variants and influenza vaccine response, the significance level should be turned to P <1.32E-03 (0.05/38 = 1.32E-03) according to Bonferroni correction. Analyses were performed by the IBM SPSS Statistics (version 25.0). Graphs were processed by GraphPad Prism v 6 (GraphPad Software, CA). The statistical power of the association analysis was estimated using the Quanto software (version 1.2.4).
ResultsCharacteristics of Subjects
In group 1, a total of 123 subjects including 41 LRs as cases and 82 matched subjects as controls received WGS. After the quality control, 12 samples were removed. Therefore, 34 cases and 77controls were finally included in the GWAS. In the verification study conducted in group 2, there were 193 cases (LRs) and 288 controls (responders). As shown in Table 1, between LRs and responders, no significant differences were observed in gender and age distribution in all three groups (P > 0.05).
Characteristics of subjects in three groups.
Variants
Group 1 (n = 111)
Group 2 (n = 481)
Combined (n = 592)
LR (%) (n = 34)
Responder (%) (n = 77)
P
LR (%) (n = 193)
Responder (%) (n = 288)
P
LR (%) (n = 227)
Responder (%) (n = 365)
P
Gender
Male
13 (32.5)
27 (67.5)
0.748
78 (40.4)
115 (59.6)
0.915
91 (39.1)
142 (60.9)
0.774
Female
21 (29.6)
50 (70.4)
115 (39.9)
173 (60.1)
136 (37.9)
223 (62.1)
Age
≤5
9 (32.1)
19 (67.9)
0.952a
45 (40.9)
65 (59.1)
0.077
54 (39.1)
84 (60.9)
0.102
6-17
3 (27.3)
8 (72.7)
13 (26.0)
37 (74.0)
16 (26.2)
45 (73.8)
18-64
13 (28.3)
33 (71.7)
94 (39.7)
143 (60.3)
107 (37.8)
176 (62.2)
≥65
9 (34.6)
17 (65.4)
41 (48.8)
43 (51.2)
50 (45.5)
60 (54.5)
Fisher exact test. LR, low responder.
GWAS and Gene-Based Association Analysis
Association P values from GWAS were presented in the Manhattan plot (Figure 2). No variant displayed genome-wide significant association (P < 1 × 10−8). The most significant SNP was LDLRAD4 rs2847111 (P = 8.82E-07) (Table 2). A total of 45 highly significant variants (P < 1 × 10−5) identified by GWAS were summarized in Supplementary Table 4, and none of them was in the coding region. Among these variants, only six of them were selected for subsequent analysis according to the function of genes and variants, the location of variants, and the LD between variants (Table 2). For variants in LD (r2 > 0.8), only one of them was selected. In order to identify more potential variants, we also considered SNPs (P < 1 × 10−2) in the coding region. Based on previous reports on the function of genes, 23 exonic SNPs of immune-related genes were selected. Table 2 summarizes the information of 29 SNPs of interest. The 29 associated SNPs indicated a total of 27 potential genes. No studies have reported that these 29 variants are associated with the immune response to influenza vaccine.
Manhattan plot of negative log of corresponding p‐values from genome‐wide association study.
Loci of interest from genome-wide association study.
For rare variants with low frequency, a gene-based association analysis was conducted. The top 27 genes (P < 0.01) identified through Deleterious (Broad) are shown in Supplementary Table 5. Two genes including AGBL1 and FAM47E were significant (P < 0.01) in the set of Deleterious (Strict). No studies have reported that these genes are associated with the immune response to influenza vaccine.
Verification Study for 29 Candidate Variants
To verify the findings in GWAS, a validation study was conducted in 193 LRs as cases and 288 responders as controls in group 2. All SNPs tested had a call rate ≥95% and conformed to HWE (P > 0.001). The MAF of SNPs in this study was consistent with the published data for the Han Chinese population in the 1000 Genomes Project data, indicating that our data set is reliable (Supplementary Table 6). In the univariate association analysis, four SNPs including SERPINA4 rs5510, IQGAP2 rs2455230, HCG22 rs2523855, and ZBTB46 rs2281929, had different genotype frequencies between the LR group and the responder group (P < 0.05). The frequency of HCG22 rs2523855 C allele (P = 0.003) and ZBTB46 rs2281929 C allele (P = 1.75E-04) in the responder group tended to be higher than that in the LR group. However, after the Bonferroni correction (P = 1.32E-03), only the distribution frequencies of ZBTB46 rs2281929 genotypes (P = 2.47E-04) and alleles (P = 1.75E-04) were still significantly different between the LR group and the responder group (Supplementary Table 6). For variants with P <0.10 in the univariate analysis, the multivariate logistic regression analyses adjusted for gender and age were conducted. All results were presented using the optimum genetic model (Table 3). Multivariate analysis results of all the 29 variants were shown in Supplementary Table 7. The results showed that compared with the TT genotype of ZBTB46 rs2281929, the TC + CC genotype was associated with a lower risk of low responsiveness to influenza vaccination adjusted for gender and age (Group 2: P = 7.75E-05, OR = 0.466, 95%CI=0.319–0.680; Combined group: P = 1.18E-06, OR = 0.423, 95%CI = 0.299–0.599). In the combined group, IQGAP2 rs2455230 GC + CC genotype was also correlated with a lower risk of low responsiveness to influenza vaccination compared with the GG genotype (P = 8.90E-04, OR = 0.535, 95%CI = 0.370–0.774), but the difference was not statistically significant in group 2 (P = 0.008).
Multivariate logistic regression analysis adjusted for age and gender.
Gene
SNPs
Genetic Model
Genotype
Group 2
Combined
P
OR (95%CI)
P
OR (95%CI)
CPED1, WNT16
rs2707520
Dominant
CC
1.00
1.00
CA+AA
0.067
0.698 (0.475-1.026)
1.39E-03
0.566 (0.400-0.803)
rs2707521
Dominant
CC
1.00
1.00
CT+TT
0.089
0.717 (0.488-1.052)
0.002
0.580 (0.410-0.821)
VPS53
rs61644407
Dominant
AA
1.00
1.00
AC+CC
0.062
1.435 (0.982-2.096)
1.42E-03
1.750 (1.241-2.468)
WWOX
rs75559202
Dominant
CC
1.00
1.00
CG+GG
0.079
1.510 (0.953-2.390)
0.004
1.846 (1.223-2.785)
SERPINA4
rs5510
Recessive
CC+CT
1.00
1.00
TT
0.007
0.289 (0.116-0.717)
0.002
0.277 (0.121-0.635)
KIDINS220
rs1044280
Recessive
CC+CA
1.00
1.00
AA
0.556
0.720 (0.241-2.153)
0.224
1.722 (0.716-4.140)
AOAH
rs2228410
Overdominant
CC+TT
1.00
1.00
CT
0.024
0.650 (0.448-0.944)
0.009
0.638 (0.455-0.896)
IQGAP2
rs2455230
Dominant
GG
1.00
1.00
GC+CC
0.008
0.577 (0.384-0.866)
8.90E-04
0.535 (0.370-0.774)
HCG22
rs2523855
Dominant
GG
1.00
1.00
GC+CC
0.009
0.606 (0.415-0.884)
0.157
0.781 (0.555-1.100)
ZBTB46
rs2281929
Dominant
TT
1.00
1.00
TC+CC
7.75E-05
0.466 (0.319-0.680)
1.18E-06
0.423 (0.299-0.599)
SNP, single nucleotide polymorphism; OR, odds ratio; CI, confidence interval.
rs2281929 TT Carriers and rs2455230 GG Carriers Had Lower Fold Rises of Antibodies
For ZBTB46 rs2281929 TT genotype carriers, the antibody fold rises against H1N1, H3N2, and B were all significantly lower than that of subjects with TC + CC genotype (P < 0.001). Compared with IQGAP2 rs2455230 GC + CC carriers, GG carriers had lower antibody fold rises to H1N1 (P = 0.001) and B (P = 0.032). The GG genotype of rs2455230 tended to be correlated with lower antibody fold rises (P = 0.096) against H3N2, but the difference was not statistically significant. (Figure 3 and Supplementary Table 8).
Comparison of antibody fold rises against different vaccine strains. (A) The comparison of antibody fold rises between ZBTB46 rs2281929 TT genotype and TC+CC genotype to H1N1, H3N2 and B. (B) The comparison of antibody fold rises between IQGAP2 rs2455230 GG genotype and GC+CC genotype to H1N1, H3N2 and B.
Verification Study for Nine Previously Reported SNPs
Nine SNPs from published studies were verified in our combined group. All SNPs tested had a call rate ≥95% and conformed to HWE (P > 0.001). The allelic and genotypic frequencies of the nine SNPs are listed in Supplementary Table 9. Supplementary Table 10 presents other information of these variants. We observed that the allelic and genotypic distributions of IL-12B rs3212227 and IL-1R1 rs3732131 were significantly different between LR group and responder group (P < 0.05). In multivariate logistic regression analysis adjusted for gender and age, compared with the TT genotype of IL-12B rs3212227, the TG + GG genotype tended to be correlated with a higher risk of low responsiveness to influenza vaccination (P = 0.016, OR = 1.614, 95%CI = 1.093–2.385). The AG + GG genotype of IL-1R1 rs3732131 tended to be associated with a better immune response to influenza vaccine compared with the AA genotype (P = 0.014, OR = 0.644, 95%CI = 0.454–0.913). However, after the Bonferroni correction, no statistical difference was found in any SNP (P < 1.32E-03) (Table 4).
Multivariate logistic regression analysis of 9 SNPs from published studies.
Reference
Gene
SNPs
Function
Genetic Model
Genotype
P
OR (95%CI)
Cummins NW et al. (16)
HMOX1
rs743811
2.8kb 3’ of HMOX1
Dominant
CC
1.00
CT+TT
0.505
1.122 (0.800-1.572)
HMOX2
rs2160567
intronic
Overdominant
TT+CC
1.00
TC
0.957
1.009 (0.722-1.412)
Poland GA et al. (6)
IL-6
rs1800796
5’-UTR
Dominant
CC
1.00
CG+GG
0.167
1.268 (0.905-1.776)
IL-12B
rs3212227
3’-UTR
Dominant
TT
1.00
TG+GG
0.016
1.614 (1.093-2.385)
IL-1R1
rs3732131
3’-UTR
Dominant
AA
1.00
AG+GG
0.014
0.644 (0.454-0.913)
IL-10RB
rs3171425
3’-UTR
Overdominant
GG+AA
1.00
GA
0.490
1.127 (0.803-1.581)
IL-10RA
rs4252249
synonymous
Overdominant
GG+AA
1.00
GA
0.179
1.449 (0.843-2.489)
IL-1RN
rs315952
synonymous
Dominant
CC
1.00
CT+TT
0.112
0.753 (0.531-1.068)
Egli A et al. (8)
IL-28B
rs8099917
7.5kb 5’ of IL28B
Overdominant
TT+GG
1.00
TG
0.323
1.334 (0.753-2.363)
SNP, single nucleotide polymorphism; OR, odds ratio; CI, confidence interval.
Discussion
Immune response to vaccine is a complicated process, which requires regulations from a series of molecules. A variant in any of these molecules may alter the level of the antibody response. Identifying these variants will make great contributions to improve the individualized influenza vaccination programs.
In this study, we focused on identifying immune-related genetic factors associated with the responsiveness to influenza vaccination. In the discovery stage, in order to identify the variants that are most likely to influence the antibody response, we selected subjects with extreme phenotypes. In the GWAS, regardless of the type of vaccine strain, only individuals whose immune responses to all vaccine strains reached neither the seroconversion rate nor the seroprotection rate can be included in the case group. Subjects whose HAI titers reached both the seroconversion rate and the seroprotection rate to all vaccine strains can be included in the control group. In the verification stage, in order to obtain a larger sample size and enhance the statistic power, we appropriately broadened the inclusion criteria of subjects and removed the criteria of seroprotection rate. Compared with the seroprotection rate, the seroconversion rate is a better index to reflect the level of immunity because the baseline antibody level before immunization was taken into consideration. Finally, two SNPs including ZBTB46 rs2281929 and IQGAP2 rs2455230 were found to be significantly correlated with the serological response to influenza vaccine. According to the MAF and OR of these two significant SNPs, our sample size was sufficient with a statistical power of 0.80 (Supplementary Table 11).
ZBTB46 is a zinc finger and BTB domain-containing transcription factor, which is selectively expressed by classical antigen presenting dendritic cells (cDCs) and their committed progenitors (17, 18). The N-terminal BTB domain is responsible for mediating protein–protein interactions, while the C-terminal is for DNA binding. Transcription of target genes can be initiated or inhibited (often inhibited) in collaboration with chromatin remodeling complexes recruited by the BTB domain (19). cDCs are immune accessory cells critical for adaptive responses against pathogens and play an important role in the process of antigen presentation (20). Researchers identified that ZBTB46 may inhibit the activation of cDCs during steady-state conditions and maintain cDCs in this quiescent state (21). When exploring the target genes of ZBTB46, Meredith et al. found that ZBTB46 binds a variety of genes in cDCs, including MHC II antigens, and behaves as a negative regulator of cDC gene expression (22). In contrast to ZBTB46, Creb1 is an activator of MHC and may compete with ZBTB46 for MHC promoter binding (23, 24). In mouse model, Meredith et al. also observed that ZBTB46-deficient cDCs expressed higher Creb1 transcript levels than wild-type controls. Compared with wild-type littermates, cDCs also expressed up to twofold higher MHC II levels in ZBTB46-deficient mice by flow cytometry. Besides, the deletion of this gene may also altered the proportions of CD4+ and CD8+ cDCs in the spleen (22). Therefore, considering the important regulatory effect of ZBTB46 on cDCs, the relationship between ZBTB46 genetic variants and vaccine immunity deserves further investigation.
ZBTB46 rs2281929 is a missense and the most significant variant in our validation study. Some studies showed that rs2281929 may be associated with the telomere length (25, 26) and glioblastoma (27), but the correlation is not strong. No previous studies have reported the relationship between rs2281929 and vaccine immunity. Our study showed that compared with the TT genotype, the TC + CC genotype of ZBTB46 rs2281929 was significantly associated with better responsiveness to influenza vaccination in both qualitative and quantitative analyses. The change in allele T to C results in a change in Threonine to Alanine. Therefore, we hypothesized that the change from Threonine to Alanine in this site may downregulate the activation of ZBTB46, reactivated the expression of cDCs as well as the antigen presentation, and finally induce a better vaccine immune response. Besides, rs2281929 is also a splicing site according to the prediction from SNPinfo (Supplementary Table 12). Further study should be performed to analyze the relationship between this variant and the function of the gene.
IQGAP2 (IQ Motif Containing GTPase Activating Protein 2) is a 180 kDa cytoplasmic multidomain scaffolding protein that belongs to the IQGAPs family. This family consists of three highly homologous members, IQGAP1, IQGAP2, and IQGAP3 (28). IQGAP proteins may play a role in regulating various cellular processes, including cytokinesis, cell migration and proliferation, intracellular signaling, vesicle trafficking, as well as cytoskeletal dynamics (29). Currently, many studies have been conducted to evaluate the relationship between IQGAP2 and cancers. And most of them reported that the mRNA expression of IQGAP2 was negatively associated with a variety of cancers, including hepatocellular carcinoma, gastric cancer, bladder cancer, breast cancer, kidney cancer, lung cancer, and so on, indicating the potential of IQGAP2 being a tumor suppressor (30–33). Results of these studies also suggested that the expression of this gene could be a biomarker for the prognosis of cancers. In addition to the role as a tumor suppressor, IQGAP2 is also found to be a novel interferon-alpha (IFN-α) antiviral effector gene and may play a role in regulating interferon stimulated genes (ISG) through the NF-κB pathway unconventionally. IFN is a regulatory of many biological processes, such as antiviral responses and immune responses. ISGs are the primary effectors of the IFN antiviral response. IQGAP2 may physically interact with RelA (a subunit of the NF-κB transcription factor, also known as p65) downstream of the IFN binding site and act cooperatively to mediate the antiviral response of IFN (34). IFITM3 is a kind of ISG. Lei et al. reported that the deletion of IFITM3 attenuated the antibody response to influenza vaccination (35). However, no association has been investigated between IQGAP2 and IFITM3. Ghaleb et al. found that knockout of IQGAP2 gene in mice has a protective effect on dextran sulfate sodium (DSS)-induced colitis in mice. In IQGAP2-deficient mice, the p65 subunit of NF-κB diminished, and the levels of macrophages and neutrophils decreased in comparison to those of wild-type mice (36). Collectively, all these results point to a potential role of IQGAP2 in the immune system. It is worthy to analyze the function of this gene in the antibody response.
As the second significant SNP identified from our study, IQGAP2 rs2455230 is a missense variant that may result in a change of Leucine to Phenylalanine from the change of G allele to C allele. In addition, the potential of rs2455230 acting as a splicing variant is also predicted from the SNPinfo (Supplementary Table 12). Our study observed that rs2455230 GC + CC genotype was associated with a lower risk of low responsiveness to influenza vaccine compared with the GG genotype. Thus, we presumed that the change from Leucine to Phenylalanine in this site may enhance the serological response to vaccine through the mediation of IQGAP2.
According to published studies on the relationship between genetic variants and the immune response to influenza vaccine, nine significant SNPs were validated in our study but no correlation was found (P < 1.32E-03). In Supplementary Table 10, we summarized the characteristics of three published studies mentioned in our verification analysis. We presumed that the reasons why our study results differ from previously published results are mainly caused by differences in ethnicity and study design. Whether these nine SNPs have any impact on the antibody response to influenza vaccine remains to be confirmed.
The limitations of our study were summarized as follows. Firstly, our GWA stage included only 111 subjects, which is a small sample size and may lead to false-negative results and the inability to identify infrequent variants. Secondly, only gender and age were included as confounding factors in the multivariate analysis. Further study should collect and contain more factors that may influence the immunogenicity of influenza vaccines. Finally, we observed that ZBTB46 rs2281929 and IQGAP2 rs2455230 are associated with the antibody response to influenza vaccine. However, the exact role of these two genes and SNPs in vaccine immunity remains unclear and requires further exploration by experimental animal models. To deal with the limitations, a series of processes were performed. In GWAS, we adopted strict inclusion criteria of subjects and a standardized procedure of quality control. In the stage of SNP selection, we carefully screened the variants considering the gene function. Then, a replication study with a larger sample size was performed to validate the results and filter out the false negative variants. We also corrected the P-value statistically to control the probability of making mistakes.
To the best of our knowledge, this is the first GWAS conducted to identify genetic variants associated with the immune response to influenza vaccination in the Chinese population. In conclusion, we have identified two novel candidate missense variants, rs2281929 and rs2455230, located in the exon region of ZBTB46 and IQGAP2 respectively. Further studies will be carried out to explore the mechanism of the regulatory function of these two genes in vaccine immunity.
Data Availability Statement
The raw sequence data reported in this paper have been deposited in the China National GeneBank DataBase (CNGBdb), under accession number CNP0001871 that are publicly accessible at https://db.cngb.org/search/project/CNP0001871/.
Ethics Statement
The studies involving human participants were reviewed and approved by the Ethics Review Committee of the National Institute for Viral Disease Control and Prevention (NIVDC, assurance number, 200916), and written informed consent was obtained from all volunteers. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.
Author Contributions
YS and DW designed the study. SW and HW managed the cohort. SW, ML and SZ performed the experiments. YC and WH contributed in collecting the samples and epidemiology data. QL and SW processed the data. SW and YS are responsible for all the analyses and write this paper. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by the Shenzhen science and technology program (Grant number: kqtd20180411143323605) and The National Natural Science Foundation of China (Grant number: 82041043).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We are grateful for the valuable contribution of Yunnan Center for Disease Control and Prevention, Urumqi Center for Disease Control and Prevention and Mr. Congda Wang.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.664024/full#supplementary-material
ReferencesHswenYBrownsteinJSLiuJHawkinsJB. Use of a Digital Health Application for Influenza Surveillance in China. Am J Public Health (2017) 107(7):1130–6. doi: 10.2105/ajph.2017.303767CDC. Seasonal Influenza Vaccine Effectiveness, 2005-2019 (2019). Available at: https://www.cdc.gov/flu/professionals/vaccination/effectiveness-studies.htm.MackenzieJSWetherallJDFimmelPJHawkinsBRDawkinsRL. Host Factors and Susceptibility to Influenza A Infection: The Effect of ABO Blood Groups and HL-A Antigens. Dev Biol Stand (1977) 39:355–62.MossAJGaughranFPKarasuAGilbertASMannAJGelderCM. Correlation Between Human Leukocyte Antigen Class II Alleles and HAI Titers Detected Post-Influenza Vaccination. PloS One (2013) 8(8):e71376. doi: 10.1371/journal.pone.0071376GelderCMLambkinRHartKWFlemingDWilliamsOMBunceM. Associations Between Human Leukocyte Antigens and Nonresponsiveness to Influenza Vaccine. J Infect Dis (2002) 185(1):114–7. doi: 10.1086/338014PolandGAOvsyannikovaIGJacobsonRM. Immunogenetics of Seasonal Influenza Vaccine Response. Vaccine (2008) 26 Suppl 4:D35–40. doi: 10.1016/j.vaccine.2008.07.065NarwaneyKJGlanzJMNorrisJMFingerlinTEHokansonJERewersM. Association of HLA Class II Genes With Clinical Hyporesponsiveness to Trivalent Inactivated Influenza Vaccine in Children. Vaccine (2013) 31(7):1123–8. doi: 10.1016/j.vaccine.2012.12.026EgliASanterDMO’SheaDBarakatKSyedbashaMVollmerM. Il-28B is a Key Regulator of B- and T-Cell Vaccine Responses Against Influenza. PloS Pathog (2014) 10(12):e1004556. doi: 10.1371/journal.ppat.1004556LiHDurbinR. Fast and Accurate Short Read Alignment With Burrows-Wheeler Transform. Bioinformatics (2009) 25(14):1754–60. doi: 10.1093/bioinformatics/btp324McKennaAHannaMBanksESivachenkoACibulskisKKernytskyA. The Genome Analysis Toolkit: A MapReduce Framework for Analyzing Next-Generation DNA Sequencing Data. Genome Res (2010) 20(9):1297–303. doi: 10.1101/gr.107524.110WangKLiMHakonarsonH. ANNOVAR: Functional Annotation of Genetic Variants From High-Throughput Sequencing Data. Nucleic Acids Res (2010) 38(16):e164. doi: 10.1093/nar/gkq603PriceALPattersonNJPlengeRMWeinblattMEShadickNAReichD. Principal Components Analysis Corrects for Stratification in Genome-Wide Association Studies. Nat Genet (2006) 38(8):904–9. doi: 10.1038/ng1847CNIC. Technical Guidelines for National Influenza Surveillance (2017). Available at: http://ivdc.chinacdc.cn/cnic/zyzx/jcfa/201709/P020170930331067634607.pdf.PurcellSNealeBTodd-BrownKThomasLFerreiraMABenderD. PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am J Hum Genet (2007) 81(3):559–75. doi: 10.1086/519795LeeSEmondMJBamshadMJBarnesKCRiederMJNickersonDA. Optimal Unified Approach for Rare-Variant Association Testing With Application to Small-Sample Case-Control Whole-Exome Sequencing Studies. Am J Hum Genet (2012) 91(2):224–37. doi: 10.1016/j.ajhg.2012.06.007CumminsNWWeaverEAMaySMCroattAJForemanOKennedyRB. Heme Oxygenase-1 Regulates the Immune Response to Influenza Virus Infection and Vaccination in Aged Mice. FASEB J: Off Publ Fed Am Soc Exp Biol (2012) 26(7):2911–8. doi: 10.1096/fj.11-190017SatpathyATKcWAlbringJCEdelsonBTKretzerNMBhattacharyaD. Zbtb46 Expression Distinguishes Classical Dendritic Cells and Their Committed Progenitors From Other Immune Lineages. J Exp Med (2012) 209(6):1135–52. doi: 10.1084/jem.20120030SatpathyATBrownRAGomuliaEBriseñoCGMumbachMRPanZ. Expression of the Transcription Factor ZBTB46 Distinguishes Human Histiocytic Disorders of Classical Dendritic Cell Origin. Mod Pathol (2018) 31(9):1479–86. doi: 10.1038/s41379-018-0052-4BeaulieuAMSant’AngeloDB. The BTB-ZF Family of Transcription Factors: Key Regulators of Lineage Commitment and Effector Function Development in the Immune System. J Immunol (2011) 187(6):2841–7. doi: 10.4049/jimmunol.1004006SteinmanRM. Decisions About Dendritic Cells: Past, Present, and Future. Annu Rev Immunol (2012) 30:1–22. doi: 10.1146/annurev-immunol-100311-102839AmonLLehmannCHKBaranskaASchoenJHegerLDudziakD. Transcriptional Control of Dendritic Cell Development and Functions. Int Rev Cell Mol Biol (2019) 349:55–151. doi: 10.1016/bs.ircmb.2019.10.001MeredithMMLiuKKamphorstAOIdoyagaJYamaneAGuermonprezP. Zinc Finger Transcription Factor zDC is a Negative Regulator Required to Prevent Activation of Classical Dendritic Cells in the Steady State. J Exp Med (2012) 209(9):1583–93. doi: 10.1084/jem.20121003KrawczykMSeguín-EstévezQLeimgruberESperisenPSchmidCBucherP. Identification of CIITA Regulated Genetic Module Dedicated for Antigen Presentation. PloS Genet (2008) 4(4):e1000058. doi: 10.1371/journal.pgen.1000058HandunnetthiLRamagopalanSVEbersGCKnightJC. Regulation of Major Histocompatibility Complex Class II Gene Expression, Genetic Variation and Disease. Genes Immun (2010) 11(2):99–112. doi: 10.1038/gene.2009.83SuraweeraNMouradovDLiSJorissenRNHampsonDGhoshA. Relative Telomere Lengths in Tumor and Normal Mucosa Are Related to Disease Progression and Chromosome Instability Profiles in Colorectal Cancer. Oncotarget (2016) 7(24):36474–88. doi: 10.18632/oncotarget.9015BerglundKReynoldsCAPlonerAGerritsenLHovattaIPedersenNL. Longitudinal Decline of Leukocyte Telomere Length in Old Age and the Association With Sex and Genetic Risk. Aging (Albany NY) (2016) 8(7):1398–415. doi: 10.18632/aging.100995ZhaoYYunDZouXJiangTLiGHuL. Whole Exome-Wide Association Study Identifies a Missense Variant in SLC2A4RG Associated With Glioblastoma Risk. Am J Cancer Res (2017) 7(9):1937–47.WhiteCDBrownMDSacksDB. Iqgaps in Cancer: A Family of Scaffold Proteins Underlying Tumorigenesis. FEBS Lett (2009) 583(12):1817–24. doi: 10.1016/j.febslet.2009.05.007HedmanACSmithJMSacksDB. The Biology of IQGAP Proteins: Beyond the Cytoskeleton. EMBO Rep (2015) 16(4):427–46. doi: 10.15252/embr.201439834XiaFDWangZLChenHXHuangYLiJDWangZM. Differential Expression of IQGAP1/2 in Hepatocellular Carcinoma and Its Relationship With Clinical Outcomes. Asian Pac J Cancer Prev (2014) 15(12):4951–6. doi: 10.7314/apjcp.2014.15.12.4951WhiteCDKhuranaHGnatenkoDVLiZOdzeRDSacksDB. IQGAP1 and IQGAP2 are Reciprocally Altered in Hepatocellular Carcinoma. BMC Gastroenterol (2010) 10:125. doi: 10.1186/1471-230x-10-125KumarDHassanMKPattnaikNMohapatraNDixitM. Reduced Expression of IQGAP2 and Higher Expression of IQGAP3 Correlates With Poor Prognosis in Cancers. PLoS One (2017) 12(10):e0186977. doi: 10.1371/journal.pone.0186977XuLShaoYRenLLiuXLiYXuJ. Iqgap2 Inhibits Migration and Invasion of Gastric Cancer Cells Via Elevating Ship2 Phosphatase Activity. Int J Mol Sci (2020) 21(6):1968. doi: 10.3390/ijms21061968BrisacCSalloumSYangVSchaeferEAHolmesJAChevaliezS. IQGAP2 is a Novel Interferon-Alpha Antiviral Effector Gene Acting Non-Conventionally Through the NF-κb Pathway. J Hepatol (2016) 65(5):972–9. doi: 10.1016/j.jhep.2016.06.028LeiNLiYSunQLuJZhouJLiZ. IFITM3 Affects the Level of Antibody Response After Influenza Vaccination. Emerg Microbes Infect (2020) 9(1):976–87. doi: 10.1080/22221751.2020.1756696GhalebAMBialkowskaABSniderAJGnatenkoDVHannunYAYangVW. Iq Motif-Containing Gtpase-Activating Protein 2 (Iqgap2) Is a Novel Regulator of Colonic Inflammation in Mice. PloS One (2015) 10(6):e0129314. doi: 10.1371/journal.pone.0129314