Identification of Two Novel Candidate Genetic Variants Associated With the Responsiveness to Influenza Vaccination

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.

Introduction

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 Methods

Subjects 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.

FIGURE 1

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⁻⁵) were selected for further validation. To identify more potential variants, we also focus on SNPs (P <1 × 10⁻²) 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 Analysis

Discovery 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 χ² 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).

Results

Characteristics 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).

TABLE 1

Table 1 Characteristics of subjects in three groups.

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⁻⁸). The most significant SNP was LDLRAD4 rs2847111 (P = 8.82E-07) (Table 2). A total of 45 highly significant variants (P < 1 × 10⁻⁵) 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 (r² > 0.8), only one of them was selected. In order to identify more potential variants, we also considered SNPs (P < 1 × 10⁻²) 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.

FIGURE 2

Figure 2 Manhattan plot of negative log of corresponding p‐values from genome‐wide association study.

TABLE 2

Table 2 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).

TABLE 3

Table 3 Multivariate logistic regression analysis adjusted for age and gender.

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).

FIGURE 3

Figure 3 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).

TABLE 4

Table 4 Multivariate logistic regression analysis of 9 SNPs from published studies.

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

References

1. Hswen Y, Brownstein JS, Liu J, Hawkins JB. 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.303767

PubMed Abstract | CrossRef Full Text | Google Scholar

2. CDC. Seasonal Influenza Vaccine Effectiveness, 2005-2019 (2019). Available at: https://www.cdc.gov/flu/professionals/vaccination/effectiveness-studies.htm.

Google Scholar

3. Mackenzie JS, Wetherall JD, Fimmel PJ, Hawkins BR, Dawkins RL. 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.

PubMed Abstract | Google Scholar

4. Moss AJ, Gaughran FP, Karasu A, Gilbert AS, Mann AJ, Gelder CM, et al. 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.0071376

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gelder CM, Lambkin R, Hart KW, Fleming D, Williams OM, Bunce M, et al. Associations Between Human Leukocyte Antigens and Nonresponsiveness to Influenza Vaccine. J Infect Dis (2002) 185(1):114–7. doi: 10.1086/338014

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Poland GA, Ovsyannikova IG, Jacobson RM. Immunogenetics of Seasonal Influenza Vaccine Response. Vaccine (2008) 26 Suppl 4:D35–40. doi: 10.1016/j.vaccine.2008.07.065

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Narwaney KJ, Glanz JM, Norris JM, Fingerlin TE, Hokanson JE, Rewers M, et al. 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.026

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Egli A, Santer DM, O’Shea D, Barakat K, Syedbasha M, Vollmer M, et al. 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.1004556

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Li H, Durbin R. Fast and Accurate Short Read Alignment With Burrows-Wheeler Transform. Bioinformatics (2009) 25(14):1754–60. doi: 10.1093/bioinformatics/btp324

PubMed Abstract | CrossRef Full Text | Google Scholar

10. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. 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.110

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wang K, Li M, Hakonarson H. ANNOVAR: Functional Annotation of Genetic Variants From High-Throughput Sequencing Data. Nucleic Acids Res (2010) 38(16):e164. doi: 10.1093/nar/gkq603

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal Components Analysis Corrects for Stratification in Genome-Wide Association Studies. Nat Genet (2006) 38(8):904–9. doi: 10.1038/ng1847

PubMed Abstract | CrossRef Full Text | Google Scholar

13. CNIC. Technical Guidelines for National Influenza Surveillance (2017). Available at: http://ivdc.chinacdc.cn/cnic/zyzx/jcfa/201709/P020170930331067634607.pdf.

Google Scholar

14. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. 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/519795

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lee S, Emond MJ, Bamshad MJ, Barnes KC, Rieder MJ, Nickerson DA, et al. 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.007

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Cummins NW, Weaver EA, May SM, Croatt AJ, Foreman O, Kennedy RB, et al. 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-190017

CrossRef Full Text | Google Scholar

17. Satpathy AT, Kc W, Albring JC, Edelson BT, Kretzer NM, Bhattacharya D, et al. 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.20120030

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Satpathy AT, Brown RA, Gomulia E, Briseño CG, Mumbach MR, Pan Z, et al. 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-4

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Beaulieu AM, Sant’Angelo DB. 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.1004006

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Steinman RM. Decisions About Dendritic Cells: Past, Present, and Future. Annu Rev Immunol (2012) 30:1–22. doi: 10.1146/annurev-immunol-100311-102839

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Amon L, Lehmann CHK, Baranska A, Schoen J, Heger L, Dudziak D. Transcriptional Control of Dendritic Cell Development and Functions. Int Rev Cell Mol Biol (2019) 349:55–151. doi: 10.1016/bs.ircmb.2019.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Meredith MM, Liu K, Kamphorst AO, Idoyaga J, Yamane A, Guermonprez P, et al. 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.20121003

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Krawczyk M, Seguín-Estévez Q, Leimgruber E, Sperisen P, Schmid C, Bucher P, et al. Identification of CIITA Regulated Genetic Module Dedicated for Antigen Presentation. PloS Genet (2008) 4(4):e1000058. doi: 10.1371/journal.pgen.1000058

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Handunnetthi L, Ramagopalan SV, Ebers GC, Knight JC. 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.83

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Suraweera N, Mouradov D, Li S, Jorissen RN, Hampson D, Ghosh A, et al. 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.9015

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Berglund K, Reynolds CA, Ploner A, Gerritsen L, Hovatta I, Pedersen NL, et al. 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.100995

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhao Y, Yun D, Zou X, Jiang T, Li G, Hu L, et al. Whole Exome-Wide Association Study Identifies a Missense Variant in SLC2A4RG Associated With Glioblastoma Risk. Am J Cancer Res (2017) 7(9):1937–47.

PubMed Abstract | Google Scholar

28. White CD, Brown MD, Sacks DB. Iqgaps in Cancer: A Family of Scaffold Proteins Underlying Tumorigenesis. FEBS Lett (2009) 583(12):1817–24. doi: 10.1016/j.febslet.2009.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Hedman AC, Smith JM, Sacks DB. The Biology of IQGAP Proteins: Beyond the Cytoskeleton. EMBO Rep (2015) 16(4):427–46. doi: 10.15252/embr.201439834

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Xia FD, Wang ZL, Chen HX, Huang Y, Li JD, Wang ZM, et al. 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.4951

PubMed Abstract | CrossRef Full Text | Google Scholar

31. White CD, Khurana H, Gnatenko DV, Li Z, Odze RD, Sacks DB, et al. IQGAP1 and IQGAP2 are Reciprocally Altered in Hepatocellular Carcinoma. BMC Gastroenterol (2010) 10:125. doi: 10.1186/1471-230x-10-125

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Kumar D, Hassan MK, Pattnaik N, Mohapatra N, Dixit M. 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.0186977

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Xu L, Shao Y, Ren L, Liu X, Li Y, Xu J, et al. 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/ijms21061968

CrossRef Full Text | Google Scholar

34. Brisac C, Salloum S, Yang V, Schaefer EA, Holmes JA, Chevaliez S, et al. 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.028

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Lei N, Li Y, Sun Q, Lu J, Zhou J, Li Z, et al. IFITM3 Affects the Level of Antibody Response After Influenza Vaccination. Emerg Microbes Infect (2020) 9(1):976–87. doi: 10.1080/22221751.2020.1756696

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ghaleb AM, Bialkowska AB, Snider AJ, Gnatenko DV, Hannun YA, Yang VW, et al. 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

PubMed Abstract | CrossRef Full Text | Google Scholar