Main

The design of the GenOMICC study and the rationale for focusing on critical illness has been previously described1,2. In brief, patients with confirmed COVID-19 requiring continuous cardiorespiratory monitoring or organ support (a generalizable definition for critical illness) were recruited in 2020–2022. We first performed ancestry-specific GWAS analyses according to the methods that we described previously1,2. Using the results of these GWAS analyses, previously reported results obtained using GenOMICC participants with whole-genome sequencing data2 and data from GenOMICC Brazil, we performed trans-ancestry and -platform meta-analyses within the GenOMICC study for a critically ill COVID-19 phenotype and a hospitalized COVID-19 phenotype (Extended Data Fig. 1). The results of these GenOMICC-only meta-analyses are presented for both critically ill and hospitalized phenotypes (Table 1 and Extended Data Fig. 2). To put these results into the context of existing knowledge, we performed comprehensive meta-analyses, drawing on further GWAS results, including data shared by the SCOURGE consortium and published data from the COVID-19 Human Genetics Initiative (HGIv6, 2021)4. The characteristics of the contributing studies are summarized in Supplementary Tables 13 and 14 for the critically ill and hospitalized phenotypes, with further details on each study provided in the Supplementary Information. We used a mathematical subtraction approach, as done in our previous work2, to remove signals of previous GenOMICC releases from HGIv6, yielding an independent dataset.

Table 1 Genome-wide significant associations with critical COVID-19, listing independent lead variants
Full size table

As no replication cohorts exist for these meta-analyses, we used the heterogeneity across studies to assess the reliability of individual findings (Supplementary Table 15). Owing to the unusually extreme phenotype in the GenOMICC study, some heterogeneity is expected for the strongest associations when compared with studies with more permissive inclusion criteria. Importantly, significant heterogeneity was not detected for any of the findings that we report here (Supplementary Table 15). Comparing effect estimates between studies using a regression approach that takes into account estimation errors (Methods), we detected systematic differences in effect sizes between studies (Extended Data Fig. 3). For example, effects for the HGI critical illness phenotype (which was designed to parallel the GenOMICC inclusion criteria) are smaller than those obtained using prospective recruitment in GenOMICC by a factor of 0.68. As the effect sizes in GenOMICC are consistently larger than other studies, and GenOMICC contributes a disproportionately large signal to meta-analyses of both critical and hospitalized phenotypes (Extended Data Fig. 4), between-study heterogeneity is likely to reflect the careful case ascertainment and extreme phenotype in GenOMICC compared with other studies.

We found 49 common genetic associations with critical COVID-19 meeting our criteria for genome-wide significance in the absence of heterogeneity (Extended Data Fig. 2 and Table 1). Findings from previous reports were consistently replicated (Extended Data Table 2). Conditional analysis revealed two additional lead variants (Table 1) and statistical fine-mapping provided credible sets of putative causal variants for a majority of lead variants (Supplementary Figs. 27–44 and Supplementary Table 5). Gene-level analyses found 196 significantly associated genes at a Bonferroni-corrected threshold (Supplementary Table 10). There were no genome-wide significant differences in the effects between sexes in a sex-stratified meta-analysis using a subset of cohorts (Supplementary Fig. 1).

Therapeutic implications

Our analysis is limited to common variants that are detectable on genotyping arrays and imputation panels. Although most lead variants are not directly causal, in some cases, they highlight molecular mechanisms that alter clinical outcomes in COVID-19, and may have direct therapeutic relevance. To investigate the disease mechanisms, we first quantified the effect of inferred gene expression on critical illness in three relevant tissue/cell types. Many of the genes that we have found to be implicated in critical COVID-19 (refs. 1,2) are highly expressed in the monocyte–macrophage system, which has poor coverage in existing expression quantitative trait loci (eQTL) datasets. For this reason, we constructed a new TWAS model in primary monocytes obtained from 176 individuals (Methods). We found significant associations after Bonferroni correction between critical COVID-19 and predicted gene expression in lung (33), blood (21), monocyte (37) and all-tissue (107) meta-analysis (Supplementary Table 2 and Supplementary Table 11). We extended these findings using generalized summary-level data Mendelian randomization (GSMR) for RNA expression (Fig. 2, Extended Data Table 1, Supplementary Figs. 11–18 and Supplementary Table 4).

In parallel, we assessed the effect of genetically determined variation in circulating protein levels on the critical illness phenotype using GSMR5. We identified 15 unique proteins linked to critical illness, as summarized in Extended Data Table 1 (Supplementary Table 3). Of the significant results, we found causal evidence implicating five new proteins in comparison to our previous GSMR analysis2: QSOX2, CREB3L4, myeloperoxidase (MPO), ADAMTS13 and mannose-binding lectin-2 (MBL2) (Supplementary Fig. 10). These include well-studied biomarkers and potential drug targets in sepsis—the innate immune pattern recognition receptor MBL2 and the neutrophil effector enzyme MPO. ADAMTS13 modulates von Willebrand-factor-mediated platelet thrombus formation and may have a role in the hypercoagulable state in critical COVID-19 (Extended Data Fig. 5).

Three genes containing non-synonymous protein-coding changes associated with severe disease were also found to have significant effects from differential gene expression: SLC22A31 (ref. 2) (Fig. 1), SFTPD4 (Fig. 1) and TKY2 (ref. 1) (Extended Data Fig. 6). Further biological and clinical research will be required to dissect the genetic evidence at these loci. In the example of TYK2, there is now a therapeutic test of the genetic predictions. Our previous report of association between higher expression and critical illness1 led directly to the inclusion of a new drug, baricitinib, in a large clinical trial; the result demonstrated a clear therapeutic benefit3. This therapeutic signal is consistent across multiple trials, providing the first proof-of-concept for drug target identification using genetics in critical illness and infectious disease.

Fig. 1: Functional genomics analyses for SLC22A31 and SFTPD.
figure 1

a, Effect-size plot for the effect of multiple variants on SLC22A31 expression (eQTLgen, x axis) against increasing susceptibility to critical COVID-19 (βxy = 0.11; Pxy = 1.3 × 10−9). The colour shows linkage disequilibrium (LD) with the missense variant rs117169628. b, Three cartoon views of an AlphaFold22 model of putative solute carrier family 22 member 31 (SLC22A31; UniProtKB: A6NKX4). The side chains of Pro474 and interacting amino acids are shown as connected spheres. A putative channel for small-molecule transport across the cell membrane is indicated by a dashed circle. Pro474 is predicted to be located in the transmembrane helix and point towards a putative transport pathway of a small molecule. The risk variant, P474L (Ala at rs117169628) would be expected to introduce more flexibility to the transmembrane helix and might therefore affect the transport properties of SLC22A31. Pro474 is predicted to be in a tightly packed environment, and may therefore affect the folding of SLC22A31. c, Effect-size plot for effect of multiple variants on SFTPD expression (eQTLgen, x axis) against increasing susceptibility to critical COVID-19 (βxy = 0.16; Pxy = 9.7 × 10−6). Colour shows linkage disequilibrium with the missense variant rs721917. d, Three cartoon views of an AlphaFold22 model of pulmonary surfactant-associated protein D (SFTPD; UniProtKB: P35247). The side chain of the variant Met31 is shown as connected spheres. Met31 is predicted to be located in the secondary-structure-lacking region of SFTPD. In the diagram on the right, oxygen and nitrogen atoms are coloured red and blue respectively, and the sulfur atom is coloured yellow.

Full size image

To assess the immediate therapeutic use of our results for repurposing of existing compounds, we considered the drug therapies under consideration by the UK COVID-19 Therapeutic Advisory Panel (UK-CTAP), a national independent review group supported by an expert due-diligence panel6. Consistent evidence from gene-level GWAS (Supplementary Table 6 and Supplementary Table 10) and post-GWAS analyses was identified for several licensed compounds (Supplementary Table 12). For example, we found an association in another gene encoding a protein that is inhibited by baricitinib and other JAK inhibitors—the intracellular signalling kinase, JAK1, which is stimulated by numerous cytokines including type I interferons and IL-6. Mendelian randomization analysis of RNA expression revealed a significant positive association between the expression of the gene encoding a canonical inflammatory cytokine, tumour necrosis factor (TNF), and severe disease (Fig. 2). This suggests that inhibition of TNF signalling may be an effective therapy in severe COVID-19.

Fig. 2: GSMR effect sizes.
figure 2

a,b, The predicted effect of change in protein concentration (a) and gene expression (b) on the risk of critical COVID-19 is shown for proteins and genes significantly linked to critical COVID-19 by GSMR (false-discovery rate (FDR) < 0.01). The bars show 95% confidence intervals.

Full size image

Our additional expression data in monocytes reveal a marked tissue-specific effect on expression of PDE4A. This phosphodiesterase regulates the production of multiple inflammatory cytokines by myeloid cells. In contrast to the negative correlations seen in the lungs and blood, we show that a genetic tendency for higher expression of PDE4A in monocytes is associated with critical COVID-19 (Supplementary Table 11). Inhibition of PDE4A by several existing drugs is under investigation in multiple inflammatory diseases7, reduces pulmonary endothelial permeability8 and appears to be safe in small clinical trials in patients with COVID-19.

The postulated biological role of genes associated with critical COVID-19 in GWAS, TWAS and GSMR results is shown in Extended Data Fig. 5, which highlights the preponderance of genes with expression or functions in the mononuclear phagocyte system. This includes SLC2A5, encoding the GLUT5 fructose transporter, which is strongly inducible in primary macrophages in response to inflammatory stimulation9, and XCR1, a dendritic cell receptor with a critical role in cytotoxic T cell-mediated antiviral immunity10. NPNT, a significant meta-TWAS association in the genome-wide significant region on chromosome 4 (chr4:105673359; Supplementary Table 11), encodes a pulmonary basement membrane protein that may have a protective role in acute lung injury11.

Host–pathogen interaction

Our results also demonstrate the capacity of host genetics to reveal core mechanisms of disease. Multiple genes implicated in viral entry are associated with severe disease. In addition to ACE2, we detect a genome-wide significant association in TMPRSS2, a key host protease that facilitates viral entry that we have previously studied as a candidate gene12. This effect may be viral-lineage specific13. A strong GWAS association is seen in RAB2A (Table 1), with TWAS evidence suggesting that more expression of this gene is associated with worse disease (Supplementary Table 11). RAB2A is highly ranked in our previous meta-analysis by information content14 study of host genes implicated in SARS-CoV-2 interaction using in vitro and clinical data15, and is consistent with CRISPR screen data showing that RAB2A is required for viral replication16.

Although our focus on critical illness enhances discovery power (Extended Data Fig. 4), it has the disadvantage of combining genetic signals for multiple stages in disease progression, including viral exposure, infection and replication, and development of inflammatory lung disease. From these data alone we cannot identify when in disease progression the causal effect is mediated, although clinical evidence helps to make some predictions17 (Extended Data Fig. 5). As most cases included were recruited before vaccinations and treatments became available (Extended Data Fig. 7), at present, our study does not have sufficient statistical power to dissect the genetic effects of treatments or vaccination. These effects may include the masking of true associations, or the detection of genetic effects mediated by vaccine or drug response, rather than COVID-19 susceptibility. However, the absence of divergent genetic effects between studies (Supplementary Figs. 2–5) or consistent changes in effect allele frequency among cases over time (Supplementary Figs. 45–48) suggests that treatment and vaccination have not substantially affected the association between the specific variants that we report and the risk of critical illness.

As we performed a meta-analysis of multiple studies that may have slightly different definitions of the phenotype, effect sizes differ between studies (Supplementary Figs. 2–5). This, together with ancestry-specific effects1, may explain the heterogeneity in strong GWAS signals, such as the LZTFL1 signal in Table 1. Different studies also have sets of variants that are not completely overlapping, so P values between variants in high linkage disequilibrium are more different than expected. Although most of the studies contain individuals from multiple ancestries, a large majority of the individuals are of European ancestry. In future research, there is a scientific and moral imperative to include the full diversity of human populations.

Together, these results deepen our understanding of the pathogenesis of critical COVID-19 and highlight new biological mechanisms of disease, several of which have immediate potential for therapeutic targeting.

Methods

Hospitalization meta-analysis

The hospitalized phenotype includes patients who were hospitalized with a laboratory-confirmed SARS-Cov2 infection. In this analysis we included GenOMICC, GenOMICC Brazil, GenOMICC Saudi Arabia, ISARIC4C, HGIv6 B2 phenotype with subtraction of GenOMICC data, SCOURGE hospitalized versus population and mild cases, and 23andMe broad respiratory phenotype. A summary description of each analysis is given above, a table with the included studies can be found in Supplementary Table 14 and an extended description can be found in Supplementary Table 1.

Critical illness meta-analysis

The critically ill COVID-19 group included patients who were hospitalized owing to symptoms associated with laboratory-confirmed SARS-CoV-2 infection and who required respiratory support or whose cause of death was associated with COVID-19. In the critical illness analysis, we included GenOMICC, patients with critical illness from ISARIC4C, HGIv6 phenotype A2 with subtraction of GenOMICC data, SCOURGE severity grades 3 and 4 versus population controls, and 23andMe respiratory support phenotype. A summary description of each analysis can be found above, a table with the included studies can be found in Supplementary Table 13 and an extended description can be found in Supplementary Table 1.

Meta-analyses

All meta-analyses across studies were performed using a fixed-effect inverse-variance weighting method and control for population stratification in the METAL software23. Allele frequency was calculated as the average frequency across studies with the METAL option AVERAGEFREQ. P values for heterogeneity in effect sizes between studies were calculated using a Cochran’s Q-test implemented in METAL. For variants in the same position with different REF and ALT alleles across studies, the GenoMICC variant in the European population was selected and the rest were removed. Finally, variants with switched ALT and REF alleles between HGIv6 and GenOMICC were also removed on the basis of differences in allele frequency of the alternative allele. Variants were annotated to the closest genes using dbsnp v.b151 GRCh38p7 and bionrRt R package (v.2.46.3)24. As each single-nucleotide polymorphism (SNP) of the meta-analysis can be present in different subsets of cohorts, there may be large differences in P values in SNPs with a high level of linkage disequilibrium, which may have an effect on downstream analyses. For this reason, variants that were not present in one of the three biggest studies—GenOMICC European ancestry, HGIv6 or SCOURGE—were filtered out from post-GWAS analysis.

Conditional analysis

We performed a step-wise conditional analysis to find independent signals. As European-specific data are not available in some cohorts but European ancestry is largely predominant (87.2% of cases with critical illness), we performed the conditional analysis using a European reference panel and the meta-analysis results of the whole cohort. To perform the conditional analysis, we used the GCTA (v.1.9.3) --cojo-slct function25. The parameters for the function were P = 5 × 10−8, a distance of 10,000 kb and a co-linear threshold of 0.9 (ref. 26), and the reference population for the conditional analysis was individuals of European ancestry with whole-genome sequence available in the GenOMICC study and whole genomes from the 100,000 Genomics England project2.

Credible set fine-mapping

We performed fine-mapping using the SuSiE model27 to construct credible sets for the independent signals identified using conditional analysis. As for conditional analysis, we used a European reference panel and the meta-analysis results of the whole critical illness cohort. We performed analyses in 1 Mb windows centred on the lead variants identified through conditional analysis. In cases in which windows for multiple variants overlapped, they were joined into a single window. For each window, we fitted the SuSiE summary statistics model setting the expected number of independent signals to the number of identified though conditional analysis. Models for three windows did not converge in 500 iterations and have been excluded. As a reference, we used the publically available linkage disequilibrium information for non-Finish Europeans from the GNOMAD 2.1.1 release. Full data for all variants included in credible sets are included in Supplementary Table 5.

Gene-level analysis

We performed an analysis summarizing the genetic associations at the gene level using the mBAT-combo method28. We used the COVID ‘all critical cohorts’ meta-analysis (GenOMICC, HGIv6 phenotype A2, SCOURGE and 23andMe) summary statistics. As this is a trans-ethnic meta-analysis, we used a mixed ancestry linkage disequilibrium reference panel, consisting of 3,202 1000 Genomes phase 3 samples. We considered a list of protein-coding genes with unique ensemble gene ID based on the release from GENCODE (v.40) for hg38, which can be found on the mBAT-combo website (https://yanglab.westlake.edu.cn/software/gcta/#mBAT-combo). A gene region was taken to span 50 kb upstream to 50 kb downstream of the gene’s untranslated regions.

Sex-stratified meta-analysis

To test for differences in genetic effects, we performed sex-stratified GWAS of the COVID-19 critical illness phenotype in the European ancestry GenOMICC WGS and genotyped cohorts and SCOURGE. We then performed a meta-analysis for each sex following the same methods as for the main analysis. We tested for differences in effects between the meta-analyses of the two sexes following previously described methods29.

Mendelian randomization

GSMR5 was performed. We used the COVID ‘all critical cohorts’ meta-analysis (GenOMICC, HGIv6 phenotype A2, SCOURGE and 23andMe) as the outcome, protein expression quantitative-trait loci (pQTLs) from ref. 30 and RNA expression quantitative-trait loci (eQTLs) from eQTLgen31 (2019-12-23 data release) as exposures, and 10,000 individuals of European ancestry randomly sampled from the UK Biobank as the linkage disequilibrium reference cohort (50,000 for linkage disequilibrium to missense variant plots). GSMR was performed for all exposures for which we were able to identify two or more suitable SNPs. SNPs were chosen to meet the following criteria: (1) SNP to exposure association P < 5 × 10−8; (2) linkage disequilibrium clumping lead SNPs only (±1 Mb, r2 < 0.05); (3) SNP not removed by HEIDI-outlier filtering (for the removal of SNPs with evidence of horizontal pleiotropy) at the default threshold value of 0.01. eQTLGen effect sizes and standard errors were estimated as described in supplementary note 2 of ref. 32. We considered as significant those exposure–outcome pairs with FDR < 0.05.

TWAS analysis

To perform TWAS analysis in GTExv8 tissues33, we used the MetaXcan framework and the GTExv8 eQTL and sQTL MASHR-M models available for download online (http://predictdb.org/) and the ‘all critical cohorts’ meta-analysis. We first calculated individual TWAS for whole blood and lungs using the S-PrediXcan function34,35. We next performed a metaTWAS including data from all tissues to increase the statistical power using s-MultiXcan36. We applied Bonferroni correction to the results to choose significant genes and introns for each analysis.

Monocyte gene expression

To detect eQTLs, untreated primary monocytes were prepared from 174 healthy individuals of Northern European (British) ancestry recruited through the Oxford Biobank. Poly(A) RNA was paired-end 100 bp sequenced in the Oxford Genome Centre using the Illumina HiSeq-4000 machines (median = 47,735,438 reads per sample). Reads were aligned to CRGh38/hg38 using HISAT2 with the default parameters. High mapping quality reads were selected on the basis of MAPQ score using bamtools. Duplicate reads were marked and removed using picard (v.1.105). Samtools was used to pass through the mapped reads and calculate statistics. Read count information was generated using HTSeq and normalized using DESeq2. Sample contamination and swaps were detected by comparing the imputed SNP-array genotypes with genotypes called from RNA-seq using verifyBamID. Genotyping was performed with Illumina HumanOmniExpress with coverage of 733,202 separate markers. Genotypes were pre-phased with SHAPEIT2, and missing genotypes were imputed with PBWT. Poly(A) RNA was paired-end sequenced at the Oxford Genome Centre using the Illumina HiSeq-4000 machines. vcftools (v0.1.12b) was applied on genetic variation data in the form of variant call format (VCF) files to filter out indels and SNPs with a minor allele frequency of less than 0.04.

TWAS analysis for monocyte data was performed using genotyping and monocyte RNA-sequencing data from 174 individuals. Using a region of 500 kb around each gene, we calculated gene expression models using the Fusion R package37. For each gene, three models were calculated adding as covariates the two first principal components calculated from the genotype: blup, elastic networks and lasso. The model with a better r2 between predicted and measured expression in a fivefold cross-validation was chosen. Then SNP genetic heritability was calculated for the 500 kb region for each gene and those genes with a nominal significant SNP heritability estimate (P ≤ 0.01) were chosen for the TWAS analysis. Summary statistics for the ‘all critical cohorts’ meta-analysis and the best model for each gene were then used to perform the TWAS.

Colocalization

Significant genes in the TWAS and metaTWAS were selected for a colocalization analysis using the coloc R package. The lead SNPs and a region of 200 Mb around the gene were used to colocalize with significant genes in the TWAS with eQTL summary statistics data on the region from GTExv8 lung, GTExv8 whole blood, eQTLgen or monocyte eqtl. As in our previous analysis2, we first performed a sensitivity analysis of the posterior probability of colocalization (PPH4) on the prior probability of colocalization (P12), going from P12 = 10−8 to P12 = 10−4, with the default threshold being P12 = 10−5. eQTL signal and GWAS signals were deemed to colocalize if these two criteria were met: (1) at P12 = 5 × 10−5 the probability of colocalization PPH4 > 0.5; and (2) at P12 = 10−5 the probability of independent signal (PPH3) was not the main hypothesis (PPH3 < 0.5). These criteria were chosen to allow eQTLs with weaker P values, owing to lack of power in GTEx v.8, to be colocalized with the signal when the main hypothesis using small priors was that there was not any signal in the eQTL data.

Effect comparison

We compared the estimates of effect sizes between the individual GWASs used in the meta-analysis, for all variants that were genome-wide significant in at least one of the individual GWASs. To this end, we regressed the effects obtained using critical illness and hospitalization in the SCOURGE and 23andMe cohorts, as well as the HGI meta-analyses on the effect estimates obtained using the GenOMICC cohort. To account for estimation errors present in both the dependent and independent variables of the regression we used orthogonal distance regression38.

Weight of studies

To calculate the weight of GenOMICC, we downloaded the leave-one-out data of HGIv7. As the meta-analysis is performed using a variance-weighted method, we can recover the variance for each SNP as \(v=\frac{1}{{{\rm{s.e.}}}^{2}}\), for the meta-analysis of all of the cohorts and for each one of the leave-one-out analysis. The total weight is \({w}_{{\rm{tot}}}=\frac{1}{v}\) and the weight leaving out a specific study is \({w}_{{\rm{loo}}}=\frac{1}{{v}_{{\rm{loo}}}}\). The weight of a cohort is then \({w}_{{\rm{tot}}}-{w}_{{\rm{loo}}}\). We calculated the weight for each the significant SNPs in our analysis for each study and normalized it using the total weight. Finally, we calculated the mean and s.d. from the significant SNPs for each cohort.

Forest plots

To compare effects between cohorts, we first performed a trans-ancestry meta-analysis for GenOMICC and 23andMe using METAL23. Then, we used the metagen and forest functions of the meta R package to produce forest plots for critical illness and hospitalization separately.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.