The COVID-19 GWAS data used in this study are from the COVID-19 Host Genetics Initiative (HGI), specifically from the round seven release, using European population data. The study includes data on three different COVID-19 phenotypes: very severe respiratory confirmed COVID-19 vs. population (COVID-19-a2) with a total sample size of 1,072,442 and 13,769 cases; hospitalized COVID-19 vs. population (COVID-19-b2) with a total sample size of 2,062,805 and 32,519 cases; and general COVID-19 vs. population (COVID-19-c2) with a total sample size of 2,475,240 and 122,616 cases. The release date for this data is 8 April 2022. The study encompasses contributions from numerous research partners, including but not limited to Estonian Biobank (EstBB), UK Biobank (UKBB), FinnGen, Generation Scotland, deCODE, and the Million Veterans Program (MVP). Each dataset has been filtered and standardized to ensure consistency and reliability in the analysis. Specific details and descriptions of the datasets are provided in Supplementary Table S1.

For the study on hematological traits, data was sourced from the trans-ethnic and ancestry-specific blood-cell genetics study published in Cell. This comprehensive study investigated blood-cell genetics across 746,667 individuals from five global populations (Chen et al., 2020). For our research, we focused on the European subset, including data on lymphocyte count (524923 Europeans), basophil count (474,001 Europeans), eosinophil count (474237 Europeans), neutrophil count (519,288 Europeans), and monocyte count (521,594 Europeans). Specific details of the datasets are provided in Supplementary Table S2.

The cross-ancestry genetic correlation (Rg) between each pair of COVID-19 phenotypes and blood-cell traits was evaluated using LD Score Regression (LDSC) and GWAS summary statistics. LDSC is an established method that identifies genetic correlations between complex traits and diseases, providing etiological insights and helping to prioritize likely causal relationships. The primary challenges in estimating genetic correlation from GWAS data with traditional methods include the unavailability of individual genotype data and widespread sample overlap among meta-analyses. LDSC overcomes these challenges by using a technique called cross-trait LD Score regression, which requires only GWAS summary statistics and is not biased by sample overlap (Bulik-Sullivan et al., 2015). The formula used in LDSC is as follows: E β j γ j = N 1 N 2 r g M l + N s r N 1 N 2 where β j and γ j   represent the effect sizes of S N P j on the two traits being tested, N 1 and N 2 are the sample sizes for the two traits, N s is the number of overlapping samples between the two traits, r is the phenotypic correlation in the overlapping samples, and, l j is the LD score. In this analysis, precomputed linkage disequilibrium scores for HapMap3 SNPs, derived from individuals of European ancestry in the 1000 Genomes Project (Bulik-Sullivan et al., 2015), were used. SNP markers with an imputation INFO score lower than 0.9 were excluded from the analysis (Bulik-Sullivan et al., 2015). This method allows for robust estimation of genetic correlations despite the complexities of sample overlap and lack of individual-level data.

Research has demonstrated that certain functional categories within the genome disproportionately contribute to the heritability of complex diseases (Trynka et al., 2013). Stratified LD Score Regression (s-LDSC) is a method developed to analyze these contributions by partitioning heritability from GWAS summary statistics while accounting for linkage disequilibrium (LD). This approach is computationally efficient and capable of handling large sample sizes, making it suitable for genome-wide studies involving extensive datasets. It leverages genome-wide information without relying solely on SNPs reaching genome-wide significance, allowing for a more comprehensive analysis. s-LDSC enables the identification of cell type-specific elements and other functional genomic regions that significantly contribute to the polygenic architecture of complex traits and diseases, providing insights into their biological mechanisms and helping prioritize genomic regions for further functional characterization (Finucane et al., 2015).

This study utilized annotation data constructed from six chromatin marks (DHS, H3K27ac, H3K36me3, H3K4me1, H3K4me3, and H3K9ac) across 88 cell types or tissues from the Roadmap project to partition the SNP heritability of various traits. For each chromatin mark, the cell type-specific annotations were further divided into nine categories: adipose, central nervous system, digestive system, cardiovascular, musculoskeletal and connective tissue, immune and blood, liver, pancreas, and others (Finucane et al., 2018). The annotation-specific enrichment values for each trait were converted into a color scale and visualized through hierarchical clustering.

To complement the genome-wide genetic correlation estimated by LDSC, which aggregates information across all variants in the genome, we employed ρ-HESS to quantify the local genetic correlation between pairs of traits (Finucane et al., 2018). ρ-HESS is designed to measure the correlation due to genetic variation at small regions in the genome, providing a more granular understanding of the genetic architecture. This technique only requires GWAS summary data and does not assume any specific distribution of causal variant effect sizes, while accounting for linkage disequilibrium (LD) and overlapping GWAS samples. The analysis using ρ-HESS involves several steps: first, it computes the eigenvalues of LD matrices and the squared projections of GWAS effect size vectors onto the eigenvectors of LD matrices for each trait. Next, it estimates the local SNP-heritability of each trait using these projections. Finally, ρ-HESS uses the local SNP-heritability estimates to calculate local genetic covariance estimates and their standard errors. In our study, all significant trait pairs identified in the global genetic correlation analysis (LDSC, P < 0.01) underwent local genetic correlation analysis using ρ-HESS. This allowed us to investigate which specific local genomic regions contributed to the global genetic correlation. A Bonferroni-corrected p-value of less than 0.05/n was considered statistically significant (n = Hypothesis testing quantity).

Multi-Trait Analysis of GWAS (MTAG) is a method designed for the joint analysis of summary statistics from genome-wide association studies (GWAS) of different traits, possibly including overlapping samples. MTAG enhances the power to detect loci by analyzing GWAS summary statistics from related traits together, accounting for sample overlap and incomplete genetic correlation. It begins by filtering variants to remove non-common SNPs, duplicated SNPs, or SNPs with strand ambiguity (Turley et al., 2018). MTAG then estimates pairwise genetic correlations between traits using LD Score Regression (LDSC) and utilizes these estimates to calibrate the variance-covariance matrix of the random effect component (Bulik-Sullivan et al., 2015). Following this, MTAG performs a random-effect meta-analysis to produce SNP-level summary statistics. In our application, we prioritized nominally significant pleiotropic SNPs that achieved genome-wide significance (P < 1.67 × 10⁻⁸) in the multi-trait analysis and suggestive significance (P < 1 × 10⁻³) in the original single-trait GWAS. This threshold accounts for the number of independent tests conducted based on the traditional GWAS threshold of (P < 5 × 10⁻⁸), ensuring that the significance level is adjusted to minimize the risk of false positives. This method not only boosts the discovery of associated loci but also improves the interpretability of genetic associations through more informative bioinformatics analyses, thereby providing deeper insights into the genetic architecture of complex traits. In this study, our new loci refer to genetic loci that have not been reported in the relevant subtypes.

Colocalization (Giambartolomei et al., 2014) and gene-based analysis are critical for understanding the molecular basis of associations identified in genetic studies. This methodology integrates multiple association datasets, such as those from GWAS and gene expression studies, to assess whether two association signals are consistent with a shared causal variant. In our study, we applied this approach to investigate the relationship between COVID-19 and blood-cell traits. We used the Coloc tool to perform Bayesian colocalization analysis, aiming to identify loci where association signals for COVID-19 and blood-cell traits colocalize. This analysis uses default priors and considers a colocalization significant when the posterior probability of a shared causal variant (PP.H4) is greater than 0.8.

To identify genes significantly associated with complex traits, we utilized S-PrediXcan with multiple eQTL datasets from GTEx v8 (a total of 49 tissues). The expression weights used were derived from GTEx v8 multi-tissue expression data, and were estimated and provided by Junghyun Jung from the Mancuso lab. In the S-PrediXcan analysis (Barbeira et al., 2018), we used gene expression and genetic variation (SNP) data from a small reference set of individuals to infer the expression of cis-acting genetic components in a larger set of phenotyped individuals based on their SNP genotype data. The predicted expression is presented in the form of a linear model, with weights determined by the correlation between SNPs and gene expression in the training data while accounting for linkage disequilibrium among SNPs. We associated gene expression with traits by conducting a transcriptome-wide association study (TWAS) to identify significant expression-trait associations. Based on Bonferroni correction, the genome-wide significance threshold for TWAS was set at P-TWAS <1 × 10⁻⁶ focusing on loci established to be associated with COVID-19 risk that also showed relevance in the TWAS analysis.

We investigated the genetic correlation between COVID-19 and blood cell counts (Figure 1). Interestingly, while the correlation is weak among general patients (Rg = −0.03, p = 0.36), it becomes significantly stronger in severe patients (Rg = −0.07, p = 0.04) and hospitalized patients (Rg = −0.11, p = 2 × 10⁻⁴), suggesting a possible relationship with the severity of the disease. Importantly, our findings indicate a significant negative correlation between lymphocyte counts and COVID-19. Beyond lymphocytes, no significant correlations were found between COVID-19 and white blood cell counts (Supplementary Table S3). Most viral infections are associated with an increase in lymphocyte counts. Interestingly, several epidemiological studies have shown that COVID-19 infection leads to a decrease in lymphocyte levels (Chen et al., 2023; Niu et al., 2022; Xiong et al., 2020). Notably, this finding aligns with the results we obtained from our genetic analyses.

For pulmonary COVID-19 and its subtypes, we further divided the heritability of SNP according to six chromatin markers and nine cell types. We found that each subtype of COVID-19 patients showed a similar pattern in the six chromatin markers of lymphocyte count (Supplementary Figure S1). It is worth noting that the respiratory system tissues or cell types, such as lungs, significantly enrich the chromatin markers H3K27Ac, H3K4me3 and H3K4me1 of COVID-19 or its subtypes. We also observed significant genetic enrichment of COVID-19 and its subtypes in musculoskeletal/connective tissue, digestive system, nervous system and other tissue or cell types. Interestingly, we also observed significant genetic enrichment of T cells, NK cells, and monocytes related to immune factors in the blood system, consistent with previous research findings (Ma et al., 2021; Ma et al., 2022).

Building on these significant findings, we conducted a comprehensive whole genome scan to identify specific genomic regions that may influence the shared heritability of various genetically related traits (Supplementary Table S4). During our statistical analysis, we adjusted for multiple testing to ensure the accuracy and reliability of our results. Notably, we discovered a significant correlation between COVID-19 and lymphocyte counts, particularly in the specific genomic region 19p13.2. This finding is particularly noteworthy given that previous studies have indicated multiple GWAS signals in the 19p13.2 region associated with COVID-19 (Figure 2). This research not only reinforces the importance of the 19p13.2 locus in the genetic architecture of COVID-19 but also provides new insights and a deeper understanding of the genetic mechanisms underlying this disease (Ferreira et al., 2022).

The evidence obtained from genetic correlation analysis has prompted us to further investigate the potential shared pleiotropic loci between COVID-19 and lymphocyte count. We employed Multi-Trait Analysis of GWAS (MTAG) to perform a multi-trait analysis of COVID-19 and its subtypes in relation to lymphocyte count, with special attention to SNPs achieving genome-wide significance (P < 1.66 × 10⁻⁸). To ensure the reliability of our results, we conducted comprehensive quality control (QC) analyses. The inflation lambda values for the original GWAS datasets were 1.055, 1.067, and 1.078 for severe COVID-19 (COVID-19-a2), hospitalized COVID-19 (COVID-19-b2), and general COVID-19 (COVID-19-c2), respectively, and 1.140 for lymphocyte count. The slightly higher lambda value for lymphocyte count may be attributed to its highly polygenic nature and the large sample size of the dataset, which can increase the observed inflation coefficient. In comparison, the MTAG results showed lower lambda values of 1.048 (COVID-19-a2), 1.041 (COVID-19-b2), and 1.023 (COVID-19-c2), indicating no significant inflation. Manhattan and QQ plots for these analyses are provided in Supplementary Figures S2, S3, confirming the robustness and consistency of our findings. These QC results further validate the reliability of our MTAG analysis and its ability to identify pleiotropic loci shared between COVID-19 and lymphocyte count.

In severe cases, we identified four loci, of which two are not reported previously. Notably, a signal at rs73009538 was discovered in the CARM1 gene, which is implicated in regulating gene transcription, growth, and RNA splicing, among other biological processes (Lai et al., 2021). We also identified a locus (rs55779981) within the RAVER1 gene region. Subsequently, among the four loci detected between hospitalized COVID-19 patients and lymphocyte count, two are not reported previously, specifically rs115545251 (1q22.1) and rs3181049 (19p13.3). The locus rs115545251 near the GFI1 gene. It is worth noting that another locus associated with hospitalized patients (rs3181049) and a locus identified in severe patients (rs55779981) are both located in the RAVER1 gene region. Interestingly, previous animal studies have suggested that other loci within the RAVER1 region, such as rs74956615 (19q13.2), increase susceptibility to coronavirus through complex mechanisms (Fink-Baldauf et al., 2022). In other words, these three loci, including two newly identified and one known locus, are densely located within the RAVER1 region yet remain independent of each other. For the general patient population, our LDSC analysis did not confirm a significant genetic correlation between COVID-19 and lymphocyte count. Despite the lack of significant genetic correlation, we conducted an MTAG analysis between hospitalized COVID-19 and lymphocyte counts, which revealed a previously unreported nominally significant loci (rs11065822) in the general population.

In this study, we employed the Coloc method to assess the Bayesian colocalization between COVID-19 and lymphocyte counts for each pleiotropic locus identified in the Multi-Trait Analysis of GWAS (MTAG). We identified three regions with PP.H4 values greater than 0.9, suggesting a high level of colocalization between COVID-19 and lymphocyte counts at these loci. Of the three colocalized loci, 19p13.2 was also identified in previous local genetic correlation analyses (Table1), reinforcing the consistency of our findings. This alignment not only bolsters the credibility of our genome-wide association results but also sheds new light on the genetic factors contributing to COVID-19 susceptibility, highlighting the potential for shared genetic pathways in immune response. To further explore the 19p13.2 (rs73009538) region in detail, we utilized the Locuscompare tool to visualize the association signals within this region. This analysis further confirmed the colocalization of genetic signals for COVID-19 susceptibility and lymphocyte counts at 19q13.2 (Figure 3) (Pruim et al., 2010).

We conducted a Transcriptome-Wide Association Study (TWAS) to investigate gene expression profiles associated with COVID-19 susceptibility and related phenotypes. TWAS integrates genotype and expression data, allowing for the identification of genes significantly linked to specific traits, thereby offering insights into gene regulation and disease association. Additionally, we have included Manhattan plots and QQ plots for the TWAS analyses, which can be accessed on Figshare (DOI: 10.6084/m9.figshare.28005845). In most tissues, the TWAS inflation factor (λ > 1.2) indicates moderate inflation, likely driven by the large sample size of the COVID-19 dataset. Similar levels of inflation have been observed in Alzheimer’s disease studies (Wingo et al., 2021), reflecting the increased statistical power to detect polygenic associations in large-scale datasets. These findings suggest that the observed inflation is primarily attributable to true signals rather than systematic bias.

The findings were largely consistent with earlier studies (Pathak et al., 2020) (Supplementary Tables S5–S7). In severe case A2, the eQTL rs10059611 was associated with ACSL6 in brain tissue, LncRNA in leukocytes and the stomach, microRNA in lung, and RAD50 in the uterus, pancreas, and skeletal muscle. These associations suggest that rs10059611 might have pleiotropic effects across multiple tissues, influencing diverse biological pathways relevant to severe COVID-19 outcomes. In hospitalized case B2, the eQTL rs7664615 was associated with PI4K2B in the pituitary, a gene involved in early T cell activation and immune regulation. In the thyroid, it was linked to SLC34A2, which plays a role in phosphate transport and cellular metabolism, and in the cerebellum, it was connected to ZCCHC4, a zinc finger protein involved in RNA processing and regulation. Additionally, rs72669986, located in the ARHGEF38 gene, was associated with GSTCD in the breast, esophagus, skin, and tibial nerve; SLC34A2 in the thyroid; and INTS12 in leukocytes, fibroblasts, prostate, and thyroid.

Furthermore, the TWAS analysis identified loci that provide additional insights into discoveries from GWAS studies. In severe case A2, the locus rs2236645 on chromosome 21 was associated with ATP5PO in the lung, brain, and skin, as well as MRPS6 in the esophagus and left atrium. Both genes are involved in mitochondrial function, a pathway that may play a critical role during early COVID-19 infection, as the virus relies on enhanced mitochondrial metabolism to support replication. In hospitalized case B2, rs3785632 was linked to SNHG26, a long non-coding RNA (lncRNA), and rs717624 was associated with NCOR1 in the testis, which regulates transcriptional activity. Additionally, rs3910266 was connected to SNX19 in the brain, a gene potentially involved in intracellular trafficking and cellular response mechanisms.