GWAS and meta-analysis identifies 49 genetic variants underlying critical COVID-19

Abstract

Critical illness in COVID-19 is an extreme and clinically homogeneous disease phenotype that we have previously shown¹ to be highly efficient for discovery of genetic associations². Despite the advanced stage of illness at presentation, we have shown that host genetics in patients who are critically ill with COVID-19 can identify immunomodulatory therapies with strong beneficial effects in this group³. Here we analyse 24,202 cases of COVID-19 with critical illness comprising a combination of microarray genotype and whole-genome sequencing data from cases of critical illness in the international GenOMICC (11,440 cases) study, combined with other studies recruiting hospitalized patients with a strong focus on severe and critical disease: ISARIC4C (676 cases) and the SCOURGE consortium (5,934 cases). To put these results in the context of existing work, we conduct a meta-analysis of the new GenOMICC genome-wide association study (GWAS) results with previously published data. We find 49 genome-wide significant associations, of which 16 have not been reported previously. To investigate the therapeutic implications of these findings, we infer the structural consequences of protein-coding variants, and combine our GWAS results with gene expression data using a monocyte transcriptome-wide association study (TWAS) model, as well as gene and protein expression using Mendelian randomization. We identify potentially druggable targets in multiple systems, including inflammatory signalling (JAK1), monocyte-macrophage activation and endothelial permeability (PDE4A), immunometabolism (SLC2A5 and AK5), and host factors required for viral entry and replication (TMPRSS2 and RAB2A).

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

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; P_xy = 1.3 × 10⁻⁹). The colour shows linkage disequilibrium (LD) with the missense variant rs117169628. b, Three cartoon views of an AlphaFold 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; P_xy = 9.7 × 10⁻⁶). Colour shows linkage disequilibrium with the missense variant rs721917. d, Three cartoon views of an AlphaFold 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.

Fig. 2. GSMR effect sizes.

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.

Extended Data Fig. 1. Pipeline of meta-analysis and post-GWAS analyses.

Red border indicates that the data is only available for the hospitalized phenotype, while a black border indicates that the analysis was performed for the critical illness phenotype.

Extended Data Fig. 2. Miami plots.

Meta-analysis results are shown for a) critical and b) hospitalized phenotypes. In each plot results obtained using all cohorts are shown at the top and using GenOMICC data only at the bottom. Independent lead variants in the analyses of all cohorts are annotated with associated genes. Genome-wide significant associations that have not been previously reported are indicated in bold.

Extended Data Fig. 3. Comparison of effect size estimates.

GenOMICC is compared with the critical and hospitalized phenotype definitions in the SCOURGE, 23andMe, and HGI analyses. The black line indicates the best linear fit, given by the equation in each plot, obtained using Orthogonal Distance Regression to account for estimation errors in both sets of effects in the comparison.

Extended Data Fig. 4. Study weightings for (a) critical and (b) hospitalized COVID-19.

Mean +/− standard deviation of weights assigned to each data source in meta-analyses for all significant SNPs.

Extended Data Fig. 5. Cartoon showing postulated roles for genes and mediators implicated in the pathogenesis of critical COVID-19 by GenOMICC GWAS, TWAS and Mendelian randomization.

Postulated roles for genetic variants are shown in a highly simplified format to illustrate potential roles in pathogenesis, with the shaded background indicating the hypothetical impact of the host immune response over time. Host immune processes are divided into those that are thought to play a role in controlling viral replication early in disease (orange section, showing “adaptive” response), and those implicated in driving hypoxaemic respiratory failure later in disease (green section, showing “maladaptive” response). Bold type gene names indicate a higher level of confidence in both the gene identification and the biological role.

Extended Data Fig. 6. Functional genomics analyses for TYK2.

(a) Effect size plot for effect of multiple variants on TYK2 expression (eQTLgen, x-axis) against increasing susceptibility to critical COVID-19 (${\beta }_{\mathit{xy}}$ = 0.53; $P_{\mathit{xy}}=1.2\times {10}^{-23}$). Colour shows linkage disequilibrium (LD) with the missense variant rs34536443. (b) Crystal structure of TYK2 kinase domain (Protein Data Bank ID 4GVJ) in two views that differ by a 45° rotation around a horizontal axis. The side chain of P1104 is shown as connected spheres with a nitrogen atom coloured in blue. Carbon, oxygen, nitrogen and phosphorus atoms of ATP are shown as magenta, red, blue and orange connected spheres, respectively. The N-terminal region of the kinase domain is not shown in the second view for clarity. The right-most panel shows a close view of P1104 and neighbouring residues with their side chains shown as sticks. Numbering of residues corresponds to UniProtKB entry P29597. P1104 is in the catalytic kinase domain and proximal to the ATP-binding site; TYK2 P1104A is catalytically impaired.

Extended Data Fig. 7. Steroid treatment and vaccination status.

Data are shown for a subset of GenOMICC cases who were also recruited to the ISARIC4C study in the UK.