Abstract
Objective The efficacy of nucleic acid–based vaccines against SARS-CoV-2 varies across individuals, partly due to genetic factors influencing neutralizing antibody production. In patients with systemic autoimmune diseases (SADs), this response may be further altered by immune dysregulation.
Methods We conducted a genome-wide association study (GWAS) to identify genetic variants associated with postvaccination anti–SARS-CoV-2 IgG antibody levels and to assess whether these associations differ between patients with SAD and healthy individuals.
Results The study included 165 participants (138 with SADs, 27 healthy controls), all of whom received nucleic acid–based vaccines. Antibody levels targeting the spike protein receptor-binding domain (RBD) and nucleocapsid were measured between 1 and 12 months after vaccination. GWAS results were metaanalyzed with data from a previously published GWAS with 1076 healthy individuals. We identified a novel association near RACGAP1 (rs706785; βmeta = −0.30, Pmeta = 3.85 × 10−8) and replicated a known association at HLA-DRB1 position 71 (βmeta = −0.23, Pmeta = 1.94 × 10−11). No significant interactions were observed between genotype and disease status.
Conclusion This study highlights both MHC and non-MHC genetic contributions to SARS-CoV-2 vaccine responses and suggests these effects are consistent across patients with SADs and healthy individuals, supporting standard vaccination strategies for individuals with systemic autoimmune conditions.
PLAIN LANGUAGE SUMMARY
People with autoimmune diseases often have altered immune responses, which may affect how well vaccines work for them. In this study, we looked at genetic differences that might help explain why some people, including those with systemic autoimmune diseases (like systemic lupus erythematosus or rheumatoid arthritis), produce more or fewer protective antibodies after receiving coronavirus disease 2019 vaccines. We found that certain genetic markers—both in and outside the immune system—are linked to how strong an individual’s response is to the vaccine. Importantly, these genetic influences were similar in patients with autoimmune disease and in healthy individuals. This suggests that standard vaccination strategies remain appropriate for people with autoimmune conditions, and that genetics play a role in shaping individual vaccine responses.
The global impact of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic has underscored the critical importance of effective vaccines in controlling the spread of the virus and mitigating severe disease outcomes. Vaccines developed using mRNA and adenoviral vector technologies have shown high efficacy in preventing severe coronavirus disease 2019 (COVID-19).1 A key factor in this efficacy is the level of neutralizing antibodies that target the viral spike protein, particularly its receptor-binding domain (RBD).2,3 However, there is considerable variability in the levels of neutralizing antibodies among individuals after vaccination, which may be influenced by a range of factors, including genetic predisposition. Indeed, multiple genome-wide association studies (GWAS) have identified genetic association signals within the major histocompatibility complex (MHC) region, especially with the most significant association at HLA-DRB1 amino acid position 71.4,5
Variability in antibody responses could be even more pronounced in patients with systemic autoimmune diseases (SADs), who may exhibit altered antibody production due to preexisting immune dysregulation. Genetic variants affecting antibody levels in these patients could have distinct effect sizes compared to those in healthy individuals. Moreover, if variants associated with SAD risk influence antibody levels, it could suggest a need for tailored vaccine protocols to optimize efficacy for patients with SADs.
This study aimed to identify genetic variants associated with antibody levels after SARS-CoV-2 nucleic acid vaccination through a GWAS involving patients with SADs and healthy controls. Additionally, we investigated whether SAD diagnosis modifies the effect of these genetic variants by examining disease–single-nucleotide polymorphism (SNP) interactions.
METHODS
Subjects and vaccination. The analysis-ready sample comprised 165 participants, including patients with systemic lupus erythematosus (n = 61), rheumatoid arthritis (RA; n = 29), immune-mediated kidney disease (n = 11), Sjögren disease (n = 19), idiopathic inflammatory myopathy (n = 9), and antineutrophil cytoplasmic antibody–associated vasculitis (n = 9), as well as 27 healthy controls. All patients with SAD fulfilled established criteria for their specific SAD. The IgG antibody levels against the RBD and nucleocapsid of SARS-CoV-2 were measured in heat-inactivated serum samples using a previously optimized ELISA-based approach,6 between 28 and 365 days after the last vaccination. All participants received COVID-19 mRNA vaccines, specifically produced by Pfizer-BioNTech (BNT1622b2) or Moderna (mRNA-1273), with ≥ 2 doses of the same vaccine. Regarding vaccination doses, 125 participants received 2 doses, 36 received 3 doses, and 4 received 4 doses before antibody measurements. Patients were recruited at the National Institutes of Health Clinical Center and provided informed written consent to participate.
Genetic data. DNA was purified from peripheral blood. Genotype data were generated at the Cancer Genomics Research Laboratory of Division of Cancer Epidemiology and Genetics–National Cancer Institute as part of the COVNET consortium (dceg.cancer.gov/research/how-we-study/genomic-studies/covnet) on an Infinium Global Screening Array (GSAMD-24v2-0; Illumina), which includes 712,189 markers. To ensure high-quality genetic data, several criteria were applied for variant exclusion. Variants with low minor allele counts (MACs; MAC < 4; n = 177,431) and low call rates (< 97%; n = 21,035) were excluded. Indels and G/C or T/A SNPs (4035 variants) were removed, as well as duplicates (1482 variants). After inferring ancestral origin of each genetic variant for each individual, variants violating Hardy-Weinberg equilibrium (HWE; n = 711 variants; PHWE < 0.001 in any ancestry or PHWE < 0.01 in ≥ 2 ancestries) were also excluded. Sample quality checks confirmed that all subjects demonstrated high genotyping call rates, no relatedness (up to the third degree), and consistency between self-reported and genetic sex.
Whole-genome imputation was performed using Beagle5,7 with reference phases derived from the 1000 Genomes Project (1KGP) phase III study. Postimputation quality control retained imputed variants with minor allele frequencies (MAFs) ≥ 0.005 and imputation quality scores (R2) ≥ 0.7, yielding a total of 13,566,613 variants across 165 samples. HLA imputation was conducted using Eagle v2.4 and Minimac4 with a high-resolution multiethnic HLA imputation panel to impute HLA classical alleles and amino acid residues.8
Local ancestry inference. We performed ancestry-specific segment analysis using RFMix v2 (github.com/slowkoni/rfmix) to assign local ancestry segments along fully phased imputed data for each study subject, based on ancestry-informed reference haplotypes from 5 target ancestries: European, African, Latino/admixed American, Central/South Asian, and East Asian.
For the ancestry reference haplotypes, we curated a set of 3593 individuals from the gnomAD Human Genome Diversity Project (HGDP)/1KGP callset (gnomAD v3.1.2; gnomad.broadinstitute.org/downloads), filtering to include only individuals with an ancestry inference probability of 1 for one of the target ancestries. Ancestry inference probabilities for the gnomAD HGDP/1KGP callset were previously estimated (gnomad.broadinstitute.org/news/2021-10-gnomad-v3-1-2-minor-release/) and are available through the gnomAD project. For HLA imputation data, we constructed a separate ancestry reference dataset using the same gnomAD HGDP/1KGP samples, applying HLA imputation methods as described earlier.
This approach enabled high-resolution, ancestry-specific calls across the genome and within the MHC region, allowing us to calculate ancestry-specific haplotype dosage at each variant and the ancestry-specific dosage of its effect allele in everyone.
Genetic association analysis. The genetic association analysis was conducted using a multivariate linear regression framework, as described in the Tractor method,9 to deconvolute ancestry-specific effect sizes for each variant. The association model applied for each variant is as follows:
In this model, T represents the normalized antibody titer for anti-RBD IgG; Hi is the dosage of haplotypes from the ith ancestry, Ai is the dosage of effect alleles from the ith ancestry, and Cj denotes the jth covariate, which includes the top 5 genetic principal components (PCs), sex, age, SAD diagnosis, log-transformed days between last vaccination and antibody measurement, vaccine manufacturer, number of vaccine doses, and normalized antinucleocapsid IgG titer. In this model, α0 is the intercept term, βi is the ancestry-specific effect estimate of haplotype dosages, γi is the ancestry-specific effect estimate for a given genetic variant, θj represents the effect estimate for the jth covariate, and ε is the error term. This model includes ancestry-specific haplotype dosages from n – 1 ancestries to avoid collinearity. Antibody titers were transformed using a rank-based inverse normal transformation, which was applied to both anti-RBD and antinucleocapsid titers. This transformation was essential for conducting a metaanalysis with previously published datasets.
The deconvoluted ancestry-specific effect sizes (γ) for each variant were subsequently aggregated using inverse variance–weighted metaanalysis. Effect sizes from ancestries with a haplotype dosage < 5%, an MAF < 5%, or an MAC < 3 were excluded from the metaanalysis. The aggregated effect sizes were then metaanalyzed with results from a previous GWAS conducted with 1076 individuals in a UK vaccine efficacy trial.5 In this external dataset, IgG RBD levels, normalized using a rank-based inverse normal transformation, were assessed 28 days after initial vaccination with the AZD1222 adenoviral vector vaccine encoding the spike protein.
For variants reaching genome-wide significance for antibody titers, we investigated potential interactions between genetic variants and disease diagnosis. We compared model fits between the regression model described above and an expanded model that additionally included interaction terms, as follows:
In this model, δi represents the ancestry-specific effect estimate for the interaction between a genetic variant and SAD diagnosis in the ith ancestry. Disease status (D) was encoded as 1 for patients with SAD and 0 for healthy controls.
RESULTS
We analyzed > 13 million autosomal variants from whole-genome imputation along with HLA variants, including classical alleles and amino acid residues, in a cohort of 165 individuals (138 patients with SAD and 27 healthy controls). Our study cohort included participants from diverse ancestries, and some individuals exhibited mixed genetic ancestry. To adjust for ancestral background in our association analysis, we conducted variant-level local ancestry inference for everyone using RFMix v2, with the gnomAD HGDP/1KGP callset as the reference dataset (Figure 1A). Ancestry proportions inferred from this process were associated with the top 3 genetic PCs, which explained substantial portions of the genetic variance in the study cohort (Figures 1B,C).
(A) Ancestral composition of study participants (n = 165) across 5 major population groups: European (EUR), African (AFR), Latino/admixed American (AMR), Central/South Asian (CSA), and East Asian (EAS). Vertical bars represent genetic ancestry proportions estimated using local ancestry inference. (B) Relative variance explained by the top 20 PCs in the study cohort. (C) Scatter plots of PC1, PC2, and PC3, with data points color-coded by each individual’s ancestry-specific proportions, highlighting the clustering of participants by genetic ancestry. (D) Quantile-quantile plot of P values from the metaanalysis combining this study and the data from Mentzer et al.5 The red diagonal line represents perfect agreement between observed and expected distributions, with the 95% CI shaded in gray. (E) Manhattan plot of −log10Pmeta (y-axis) for each variant according to chromosomal position (x-axis). The genome-wide significance threshold is denoted by a horizontal gray line. (F) Forest plot for rs706785, depicting the effect sizes (dots) of the G allele with 95% CIs across deconvoluted ancestries, the prior GWAS, and the metaanalysis. (G) Regional association plot for the RACGAP1 locus, displaying −log10Pmeta by chromosomal position. Variants are color-coded based on their LD with the most significant variant (rs706785), based on the 1KGP EUR reference panel. 1KGP: 1000 Genomes Project; GWAS: genome-wide association study; LD: linkage disequilibrium; PC: principal component.
Genome-wide association analysis was performed to estimate ancestry-specific effect sizes for each variant. We then metaanalyzed these deconvoluted effect sizes to obtain global effect estimates, adjusting for various covariates described in our methods section, including genetic PCs (PC1-PC5), sex, SAD or control diagnosis, vaccine dose, vaccine manufacturer, and time since the last vaccination (in log-transformed days), as well as antinucleocapsid IgG antibody levels.
Among these covariates, time since vaccination and antinucleocapsid IgG levels showed significant associations with antibody titers (Table 1). Notably, the log-transformed days since the last vaccination were inversely associated with antibody levels, suggesting a decline in titers over time. In contrast, antinucleocapsid IgG antibody levels, used as an additional control for SARS-CoV-2 exposure, were significantly associated with higher levels of anti-RBD IgG antibodies. Although the comparison did not reach significance after accounting for multiple testing, individuals who received the mRNA-1273 vaccine tended to have higher antibody titers compared to those who received the BNT1622b2 vaccine, which is consistent with expectations due to the higher mRNA amount in the mRNA-1273 vaccine.
Multivariate linear regression for the normalized level of anti–SARS-CoV-2 RBD IgG antibody.
We conducted GWAS using variant-level local ancestry information to estimate ancestry-specific effect sizes of genetic variants in multivariate linear regression and to estimate combined effect sizes through inverse variance–weighted metaanalysis. The genomic inflation factor (λ) from our genome-wide association results was 1.05, indicating minimal systemic bias and no evidence of significant underlying population substructure. We found no variants surpassing the genome-wide significance threshold (P = 5 × 10−8), likely due to the limited sample size of the study.
To improve statistical power and result robustness, we conducted a metaanalysis by integrating our data with results from a recent GWAS5 that employed the same normalized antibody measurement method, enabling comparability of genetic effect sizes across studies. The genomic inflation factor (λ) was 1.06 (Figure 1D). In this metaanalysis, we replicated the previously reported association of HLA-DRB1 amino acid position 71 with antibody titers (Pmeta = 1.94 × 10−11; Table 2; Figure 1E), achieving a slightly stronger significance than the original study. The presence of arginine or serine at this position was associated with lower antibody titers (βmeta = −0.23), with consistent findings in both our study (β = −0.25, P = 2.63 × 10−2) and the previous GWAS (β = −0.23, P = 2.38 × 10−10). Notably, position 71 is a critical site that shapes the epitope-binding surface of HLA-DR molecules and determines the shared epitope associated with RA.10 Arginine at this position is not exclusively present in shared epitope alleles nor is it associated with RA risk.11 Serine at position 71 is very rare in multiple populations.
Association of HLA-DRB1 and RACGAP1 variants with postvaccination anti–SARS-CoV-2 IgG antibody levels.
Additionally, in a metaanalysis, we identified a novel association locus near the RACGAP1 gene (Figures 1E,F,G). The most significant association was detected at rs706785, located 3.5 kb upstream of RACGAP1 (βmeta = −0.30; Pmeta = 3.85 × 10−8; Table 2; Figure 1G). Consistent effect sizes were observed across both our study (β = −0.59, P = 3.90 × 10−4) and the prior GWAS (β = −0.27, P = 4.67 × 10−6). To further investigate the locus, we conducted a fine-mapping analysis and identified 6 variants within the 90% credible set for this locus,12 based on approximate Bayesian factors with a prior variance of 0.04 in allelic effects. Among these, 3 variants located in the upstream and intronic regions of RACGAP1 accounted for 86% of the posterior probability (Table 2). Notably, these 3 variants have been reported as expression quantitative trait loci (eQTLs) for multiple genes, including RACGAP1, CERS5, LIMA1, and COX14, across various tissues (Supplementary Table S1, available with the online version of this article). This observation is also supported by chromatin interaction data, which suggest potential regulatory interactions within this locus (Supplementary Figure S1). For instance, the G allele of rs706785, associated with lower antibody titers, is also linked to decreased RACGAP1 expression levels.
Although genetic variants of RACGAP1 and an alanine or serine residue at HLA-DRB1 amino acid position 71 have not been associated with SADs, we investigated potential interactions of these variants with disease diagnosis to examine whether the observed GWAS effect sizes are distinct in patients with SAD compared to healthy individuals. Logistic regression analysis, incorporating interaction terms between disease diagnosis and ancestry-specific allelic dosages, did not reveal any significant interaction effects associated with antibody levels (P > 0.05).
DISCUSSION
In this study, we investigated the genetic underpinnings of antibody responses to nucleic acid vaccines against SARS-CoV-2 by performing a GWAS and a subsequent metaanalysis with a recent study that used the same antibody measurement and transformation method.5 We confirmed the previously reported association between HLA-DRB1 amino acid position 71 and anti-RBD IgG titers. Further, we identified a novel signal near the RACGAP1 gene, which includes eQTLs for RACGAP1 and nearby genes, underscoring the potential role of this locus in vaccine-induced immunity.
By leveraging genetic data from participants of diverse ancestries, our analysis effectively accounted for ancestral differences and minimized bias through robust methodologies, including ancestry-specific effect size estimation and Bayesian fine-mapping. These approaches enhanced our ability to identify genetic variants associated with antibody responses, despite the relatively small sample size. The integration of our findings with data from a prior GWAS5 significantly increased statistical power, enabling cross-study validation and the discovery of a novel locus associated with antibody levels. This combined analysis not only reinforced the evidence for previously implicated genetic variants but also identified variants around RACGAP1 as a potential contributor to the variability in vaccine-induced immune responses. Our results should be pursued in subsequent studies designed to identify and characterize further genetic variants that contribute to antibody responses to SARS-CoV-2 in individuals diagnosed with specific SADs, as well as interactions with immunosuppressive medications. Importantly, the inclusion of covariates such as vaccine manufacturer, number of doses, and time since vaccination allowed us to adjust for nongenetic factors that influence antibody titers, providing a clearer view of the genetic effects on SARS-CoV-2 vaccine responses. Although the inclusion of multiple covariates may raise concerns about overfitting in a small cohort, the genomic inflation factor (λ = 1.05-1.06) supported that the model was relatively well calibrated.
We also investigated potential interactions between genetic variants and disease diagnosis. Our results did not reveal significant interactions, suggesting that the genetic variants associated with antibody responses operate similarly in both patients with SADs and healthy individuals.
We identified several candidate genes regulated by credible-set variants around RACGAP1, although direct evidence linking them to vaccine efficacy or antibody levels is currently lacking. Among these, RACGAP1, a member of the Rho GTPase–activating protein family, is primarily known for critical roles in cytokinesis, cell migration, and invasion through its regulation of cytoskeletal dynamics and signaling pathways.13,14 Although RACGAP1 has been extensively studied in cancer biology, where its overexpression is associated with poor prognosis across various malignancies, its role in immune responses remains poorly understood. Previous findings correlating RACGAP1 expression in various cancer cells with immune cell infiltration suggest a potential role in modulating immune cell dynamics and responsiveness.15 In the context of vaccine-induced immunity, which requires immune cell infiltration to trigger an effective immune response at the site of injection or in proximal lymph nodes, the observed association of RACGAP1 expression–reducing alleles with lower antibody titers suggests that these alleles may weaken immune cell activation or recruitment, potentially dampening vaccine efficacy. Further functional studies and independent GWAS are needed to validate our findings and to clarify potential mechanisms and their implications for immune responses.
Other candidate genes, associated with the credible-set variants based on eQTL annotations (Supplementary Table S1), include LIMA1 and FMNL3, which are functionally related to RACGAP1 and regulate actin cytoskeleton dynamics,16,17 suggesting possible relevance to immune cell motility and cytoskeletal remodeling during antigen recognition. In contrast, other genes regulated by credible-set variants (CERS5, SMARCD1, AQP6, ASIC1, FAM186A, and RP4-605O3.4) have limited functional characterization or are poorly understood in terms of their immunological roles. Although their roles in vaccine-induced immune responses remain unclear, we cannot exclude the possibility that some of these genes contribute to the observed effects.
In conclusion, this study provides new insights into the genetic basis of vaccine-induced antibody responses, marking the discovery of the first non-MHC genetic association with SARS-CoV-2 vaccine–induced immunity. Notably, the absence of detectable interactions between the identified genetic variants and disease diagnosis suggests that, at least for the loci investigated, the heritable genetic effects on vaccine-induced immune responses in patients with SADs may be comparable to those observed in the general population. Although the genetic variants identified in this and previous GWAS explain only a modest fraction of the variability in antibody responses, future studies integrating multiomics data (eg, transcriptomics, proteomics, and epigenomics), larger and more diverse cohorts, and longitudinal antibody measurements could further elucidate the molecular determinants of vaccine efficacy. Ultimately, such efforts may help to identify individuals at risk of rapid antibody waning and inform personalized booster schedules, immune monitoring strategies, and precision vaccination approaches as SARS-CoV-2 moves into an endemic phase.
ACKNOWLEDGMENT
The authors sincerely thank the patients who generously participated in this study. We are also deeply grateful to Dr. Alexander J. Mentzer from the Wellcome Centre for Human Genetics, Nuffield Department of Medicine, University of Oxford, Oxford, UK, for providing the GWAS summary statistics from his recent study.5 The content of this publication does not necessarily reflect the views or policies of the US Department of Health and Human Services (HHS), nor does mention of trade names, commercial products, or organizations imply endorsement by the US government. The contributions of the National Institutes of Health (NIH) author(s) were made as part of their official duties as NIH federal employees, are in compliance with agency policy requirements, and are considered works of the US government. However, the findings and conclusions presented in this paper are those of the author(s) and do not necessarily reflect the views of the NIH or HHS.
Footnotes
CONTRIBUTIONS
All authors contributed to at least one of the following manuscript preparation roles: conceptualization and/or methodology, software, investigation, formal analysis, data curation, visualization, and validation and drafting or reviewing/editing the final draft. As corresponding author, MJK confirms that all authors have provided the final approval of the version to be published and takes responsibility for the affirmations regarding article submission (eg, not under consideration by another journal), the integrity of the data presented, and the statements regarding compliance with institutional review board/Declaration of Helsinki requirements.
FUNDING
This work was supported by the National Institutes of Health (NIH) through the Intramural Targeted Anti-COVID Award from National Institute of Allergy and Infectious Diseases; the Intramural Research Programs at National Institute of Arthritis and Musculoskeletal and Skin Diseases (no. ZIA AR-041199), National Institute of Diabetes and Digestive and Kidney Diseases, and the National Institute of Environmental Health Sciences; by federal funds from the National Cancer Institute, NIH (under contract no. 75N91019D00024); and by the Sabbatical Research Program of Kyung Hee University.
COMPETING INTERESTS
The authors declare no conflicts of interest associated with this article.
ETHICS AND PATIENT CONSENT
This study was approved by the NIH Institutional Review Board (IRB; no. 000207 COVID-SAD). This study adhered to all relevant ethical regulations, including compliance with the Declaration of Helsinki.
DATA AVAILABILITY
The genetic datasets generated and/or analyzed during the current study are not publicly available to protect privacy but are available from the corresponding author upon reasonable request.
- Accepted for publication November 4, 2025.
- Copyright © 2026 by the Journal of Rheumatology
This is an Open Access article, which permits use, distribution, and reproduction, without modification, provided the original article is correctly cited and is not used for commercial purposes.
REFERENCES
SUPPLEMENTARY DATA
Supplementary material accompanies the online version of this article.
