Main

As critical illness in patients with COVID-19 is caused, in part, by inflammatory injuries that affect the lungs and lung blood vessels1, there are at least two distinct biological components to the mortality risk: susceptibility to viral infection and propensity to develop harmful pulmonary inflammation. Susceptibility to life-threatening infections4 and immune-mediated diseases are both strongly heritable. In particular, susceptibility to respiratory viruses5 such as influenza6 is heritable and known to be associated with specific genetic variants7. In the case of COVID-19, one genetic locus—on chromosome 3p21.31—has been repeatedly associated with hospitalization8,9. As with other virus-associated diseases10, there are several examples of loss-of-function variants affecting essential immune processes that lead to severe disease in young people, such as TLR711 and some genes implicated in type 1 interferon signalling, including the receptor subunit IFNAR212. Genome-wide studies have the potential to reveal previously undescribed molecular mechanisms of critical illness in patients with COVID-19, which may provide therapeutic targets to modulate the host immune response to promote survival3.

Strong evidence indicates that critical illness caused by COVID-19 is qualitatively different from mild or moderate disease, even among hospitalized patients. There are multiple distinct disease phenotypes with differing patterns of symptoms13 and marked differential responses to immunosuppressive therapy2. In patients without respiratory failure, there is a trend indicating that treatment with corticosteroids is harmful, whereas among patients with critical respiratory failure, there is a substantial benefit2. On this basis, we consider patients with critical COVID-19 respiratory failure to have a distinct pathophysiology.

In the UK, the group of patients admitted to critical care is relatively homogeneous, with profound hypoxaemic respiratory failure being the archetypal presentation14. The active disease process in these patients is markedly responsive to corticosteroid therapy15 and is characterized by pulmonary inflammation including diffuse alveolar damage, influx of monocytes and macrophagess, mononuclear cell pulmonary artery vasculitis and microthrombus formation1,16.

Host-directed therapies have long been an aim for the treatment of the severe disease caused by respiratory viruses17. Identification of genetic loci associated with susceptibility to COVID-19 may lead to specific targets for the development or repurposing of drugs3.

The GenOMICC (Genetics Of Mortality In Critical Care, https://genomicc.org/) study has been recruiting patients with critical illness syndromes, including influenza, sepsis and emerging infections, for 5 years. To better understand the host mechanisms that lead to life-threatening disease in patients with COVID-19, we performed a genome-wide association study (GWAS) comparing critically ill patients with COVID-19 with controls from population genetic studies in the UK.

Characteristics of participants

Critically ill cases were recruited through the GenOMICC study in 208 UK intensive care units and hospitalized cases were recruited through the International Severe Acute Respiratory Infection Consortium (ISARIC) Coronavirus Clinical Characterization Consortium (4C) study.

DNA was extracted from whole blood and array-based genome-wide genotypes of good quality were obtained for 2,734 unique individuals (Methods). Genetic ancestry was inferred using principal component analyses and individuals from the 1000 Genomes Project were included as population references (Methods). After quality control and matching to ancestry groups, 2,244 individuals were included in the GWAS for analysis. Clinical and demographic features of these cases are shown in Extended Data Table 1. Additional clinical details for a subset of 1,069 cases for whom additional data were available are presented in Supplementary Figs. 713 and Supplementary Table 2. Cases were representative of critically ill patients with COVID-19 in the UK population14. Imputation in this multi-ancestry cohort was performed using the TOPMed reference panel.

Ancestry-matched control individuals were selected from the large population-based cohort of UK Biobank (five controls were included for each case). Controls with a known positive COVID-19 test were excluded. The inevitable presence of individuals in the control group, who may exhibit the critical illness phenotype if exposed to SARS-CoV-2 is expected to bias any associations towards the null hypothesis. GWAS was carried out separately for each ancestry group using logistic regression in PLINK and accounting for age, sex, deprivation decile based on the individual’s residential postal code and principal components of ancestry. In addition to several standard filters to minimize spurious associations (Methods), whole-genome sequencing analysis of a subset of 1,613 cases was used to filter out variants that are likely to have been badly called or imputed; 83,937 out of the 4,469,187 imputed variants that passed other quality-control filters after GWAS were thus removed. There was a high level of residual inflation in the South Asian and East Asian ancestry groups, rendering results in these subgroups unreliable (Extended Data Fig. 1, Supplementary Fig. 4). The largest ancestry group contained 1,676 individuals of European descent; this group was used for the primary analyses presented below.

GWAS results

In the primary analysis (cases of European descent from GenOMICC versus controls from the UK Biobank), after linkage-disequilibrium-based clumping, 15 independent association signals were genome-wide significant at P < 5 × 10−8 (Supplementary Fig. 1). Eight of these signals were successfully validated in GWAS analyses using two independent population genetic studies (100,000 Genomes Project and Generation Scotland: Scottish Family Health Study, hereafter Generation Scotland) as controls (Table 1) and were therefore taken forward for replication. A sex-specific GWAS among this group found no sex-specific associations (Supplementary Table 1). A trans-ethnic meta-analysis did not reveal any additional associations (Supplementary Fig. 3).

Table 1 Lead variants from independent genome-wide significant regions

Replication

As no study of critical illness in patients with COVID-19 of sufficient size is available, replication was carried out in a meta-analysis of data from 2,415 hospitalized patients with COVID-19 and 477,741 population controls from the COVID-19 Host Genetics Initiative (HGI, individuals of mixed ancestry from which cases and controls of the UK Biobank were excluded) and 1,128 cases and 679,531 controls of European ancestry from the ‘broad respiratory phenotype’ study of 23andMe Inc, which includes cases that were reported as having been placed on a ventilator, administered oxygen or having had pneumonia compared with control individuals who did not report COVID-19-positive tests. In addition to the locus on chromosome 3 that has previously been reported (rs73064425, odds ratio = 2.14, discovery P = 4.77 × 10−30), we found robust replication for previously undescribed associations in four loci from GenOMICC: a locus on chromosome 12 in the OAS gene cluster (rs10735079, odds ratio = 1.3, discovery P = 1.65 × 10−8), near TYK2 on chromosome 19 (rs74956615, odds ratio = 1.6, discovery P = 2.3 × 10−8), in DPP9 on chromosome 19 (rs2109069, odds ratio = 1.36, discovery P = 3.98 × 10−12) and a locus on chromosome 21, which contains the gene IFNAR2 (rs2236757, odds ratio = 1.28, discovery P = 4.99 × 10−8) (Fig. 1, Extended Data Table 2).

Fig. 1: Discovery GWAS and meta-analysis results.
figure 1

Miami plot showing the P values for the GenOMICC GWAS in individuals of European ancestry after validation (top) and meta-analysis including patients from the COVID-19 HGI and 23andMe (bottom). Uncorrected P values from GWAS analysis are shown. Top, the red horizontal line indicates genome-wide significance for common variants at −log10(5 × 10−8). Bottom, the red horizontal line indicates a more-stringent genome-wide significance threshold for meta-analysis variants at −log10(1 × 10−8). Insets, quantile–quantile (QQ) plots of observed versus expected P values that show genomic inflation (λ) for each analysis: GenOMICC European cohort, λ = 1.099; GenOMICC HGI–23andMe meta-analysis, λ = 1.017. Full QQ plots are shown in Extended Data Fig. 1.

Three variants, all in a region of chromosome 6 for which population stratification is difficult to control (the major histocompatibility complex) did not replicate (Extended Data Table 2, Supplementary Fig. 2). Further studies will be required to determine whether these associations are real.

To increase power for exploratory analyses, inverse-variance meta-analysis was performed between critically ill patients of European descent from GenOMICC (n = 1,676 cases, n = 8,380 controls), hospitalized patients with COVID-19 versus population controls (B2, version 2) without UK Biobank participants from the HGI (n = 2,415 cases, n = 477,741 controls) and the participants with the broad respiratory phenotype from 23andMe (n = 1,128 cases, n = 679,531 controls). This revealed one additional (unreplicated) locus in CCHCR1 at genome-wide significance (using a more-stringent threshold of P < 10−8 because of the absence of replication opportunities for the meta-analysis) (Table 2).

Table 2 Meta-analysis of overlapping SNPs between GenOMICC, HGI and 23andMe studies

Mendelian randomization

Mendelian randomization provides evidence for a causal relationship between an exposure variable and an outcome, given a set of well-characterized assumptions18. We used two-sample summary-data Mendelian randomization to assess the evidence in support of causal effects of RNA expression (Genotype-Tissue expression (GTEx) v.7, whole blood) of various genes on the odds of patients becoming critically ill due to COVID-19.

We specified an a priori list of target genes that relate to the mechanism of action of many host-targeted drugs that have been proposed for the treatment of patients with COVID-19 (Supplementary Table 3). Seven of these targets had a suitable locally acting expression quantitative trait locus (eQTL) in GTEx (v.7). Of these, IFNAR2 remained significant after Bonferroni correction for multiple testing for seven tests (β = −1.49, s.e. = 0.52, P = 0.0043) (Supplementary Table 4). There was equivocal evidence of heterogeneity (HEIDI19 P = 0.015), indicating that the effect of this variant on critical illness in COVID-19 may be mediated through another mechanism, which may lead to an under- or overestimation of the effect of IFNAR2 expression on the risk of critical illness in COVID-19.

We next performed transcriptome-wide Mendelian randomization to quantify support for unselected genes as potential therapeutic targets. Instruments (genetic variants that affect gene expression) were available for 4,614 unique Ensembl gene IDs. No genes were statistically significant after correcting for multiple comparisons in this analysis (4,614 tests). After conservative filtering for heterogeneity (HEIDI P > 0.05), the smallest Mendelian randomization P = 0.00049 was associated with a variant on chromosome 19: 10,466,123 that affects expression of TYK2. Nine other genes with nominally significant Mendelian randomization P values (P < 0.0051) were also taken forward for further analysis (Supplementary Table 5).

To replicate these findings, we tested for external evidence using a separate eQTL data set (eQTLgen)20 and GWAS (HGI B2, excluding UK Biobank participants). Mendelian randomization signals with consistent directions of effect were significant for IFNAR2 (P = 7.5 × 10−4) and TYK2 (P = 5.5 × 10−5) (Supplementary Table 6).

Transcriptome-wide association study

We performed a transcriptome-wide association study (TWAS)21,22 to link GWAS results to tissue-specific gene expression data by inferring gene expression from known genetic variants that are associated with transcript abundance (eQTL). For this analysis, we used GTEx v.8 data for two disease-relevant tissues chosen a priori: whole blood and lung samples (Fig. 2). We selected genes with P < 0.05 in these tissues and performed a combined meta-TWAS analysis23, incorporating eQTL data from other tissues in GTEx, to optimize power to detect differences in predicted expression in lung or blood.

Fig. 2: Summary of TWAS results.
figure 2

a, Gene-level Manhattan plot showing raw P-value results from meta-TWAS analysis across tissues (see Methods). The red horizontal line shows gene-level genome-wide significance at −log10(5 × 10−6). b, Z-scores showing the direction of effect for the genotype-inferred expression of transcripts that encode protein-coding genes in lung tissue (GTEx v.8). Red circles indicate genes with genome-wide significance at P < 5 × 10−6.

We discovered five genes with genome-wide significant differences in predicted expression compared to control individuals (Supplementary Table 7). This included four genes with differential predicted expression in lung tissue (three on chromosome 3, CCR2, CCR3 and CXCR6; and one on chromosome 5, MTA2B) (Supplementary Tables 810).

We used meta-analysis by information content (MAIC)24 to put these results in the context of existing biological knowledge about host–virus interactions associated with COVID-19. We combined the top 2,000 genes in metaTWAS with previous systematically compiled experimental evidence implicating human genes in SARS-CoV-2 replication and host response. MAIC derives a data-driven weighting for each gene from a range of experimental data sources in the form of gene lists, and outperforms other approaches to providing a composite of multiple lists24. We found that the GenOMICC TWAS results had greater overlap with results from transcriptomic, proteomic and CRISPR studies of host genes implicated in COVID-19 than any other data source(Extended Data Fig. 2).

Genetic correlations

We used a high-definition likelihood method25 to provide an initial estimate of heritability based on single-nucleotide polymorphisms (SNPs) (that is, the proportion of phenotypic variance captured by additive effects at common SNPs) for severe COVID-19, which was found to be 0.065 (s.e. = 0.019). We were not able to detect a significant signal for heritability in two additional analyses: first, using controls from the 100,000 Genomes Project (in which matching to the cases from GenOMICC is less close, which may limit heritability estimation) and second, in a smaller GWAS in which some cases from GenOMICC were compared with controls from the UK Biobank, using matched data for body-mass index and age where possible. This second analysis was less powerful because of the lack of close matches for many cases (n = 1,260 cases, n = 6300 controls) (Supplementary Fig. 14). Including rare variants in future analyses, which have larger numbers of cases of COVID-19, will provide a more comprehensive estimate of heritability. We also tested for genetic correlations with other traits—that is, the degree to which the underlying genetic components are shared with severe COVID-19. Using the high-definition likelihood method, we identified significant negative genetic correlations with educational attainment and intelligence. Significant positive genetic correlations were detected for a number of adiposity phenotypes including body-mass index and leg fat (Supplementary Fig. 19).

Consistent with GWAS results from other infectious and inflammatory diseases, there was a significant enrichment of strongly associated variants in promoters and enhancers26, particularly those identified by the EXaC study as under strong evolutionary selection27 (Supplementary Fig. 18). The strongest tissue-type enrichment was in spleen (which may reflect enrichment in immune cells), followed by pancreas (Supplementary Fig. 20).

Discussion

We have discovered and replicated significant genetic associations that are associated with life-threatening COVID-19 (Fig. 1). Our focus on critical illness increases the probability that some of these associations relate to the later, immune-mediated phase of COVID-19 associated with respiratory failure that requires invasive mechanical ventilation2. Notably, the GWAS approach is unbiased and genome-wide, enabling the discovery of previously undescribed pathophysiological mechanisms. Because genetic variation can be used to draw a causal inference, genetic evidence in support of a therapeutic target substantially improves the probability of successful drug development28. In particular, Mendelian randomization occupies a unique position in the hierarchy of clinical evidence29.

Patients admitted to intensive care units in the UK during the first wave of COVID-19 were, on average, younger and less burdened by comorbid illnesses than the hospitalized population14. The population studied here was defined by their propensity to critical respiratory failure due to COVID-19. GenOMICC recruited in 208 intensive care units (covering >95% of the capacity of intensive care units in the UK), ensuring that a broad spread across the genetic ancestry of patients from the UK was included (Extended Data Fig. 3).

For external replication, the nearest comparison is the analysis of hospitalized patients versus population controls in the COVID-19 HGI, and the broad respiratory phenotype dataset from 23andMe, which have been generously shared with the international community. We have also immediately made full summary statistics from GenOMICC data openly available at https://genomicc.org/data/.

Despite the differences in case definitions, new associations from our study of critical illness replicate robustly in combined data from studies of hospitalized cases of COVID-19 (Extended Data Table 2). Separately, the Mendelian randomization results that suggest a causal role for IFNAR2 and TYK2 are also statistically significant in confirmatory analyses. Our findings reveal that critical illness in COVID-19 is related to at least two biological mechanisms: innate antiviral defences, which are known to be important early in disease (IFNAR2 and OAS genes), and host-driven inflammatory lung injury, which is a key mechanism of late, life-threatening COVID-19 (DPP9, TYK2 and CCR2)2.

Interferons are canonical mediators of antiviral signalling in the host and stimulate release of many essential components of the early host response to viral infection30. Consistent with a beneficial role for type I interferons, increased expression of the interferon receptor subunit IFNAR2 reduced the odds of severe COVID-19 with Mendelian randomization discovery P = 0.0043 (seven tests); replication P = 7.5 × 10−4 (one test). Within the assumptions of Mendelian randomization, this represents evidence for a protective role of IFNAR2 in COVID-19. Rare loss-of-function mutations in IFNAR2 are associated with severe COVID-1912 and many other viral diseases31,32. This suggests that administration of interferon may reduce the probability of critical illness in COVID-19, but our evidence cannot distinguish when during disease progression of COVID-19 such a treatment may be effective. Exogenous interferon treatment did not reduce mortality in hospitalized patients in a large-scale clinical trial33, suggesting that this genetic effect may be mediated during the early phase of disease when the viral load is high.

The variant rs10735079 (chromosome 12, P = 1.65 × 10−8) lies in the interferon-inducible oligoadenylate synthetase (OAS) gene cluster (OAS1, OAS2 and OAS3) (Fig. 1). Our TWAS detected significant associations with predicted expression of OAS3 (Fig. 2). OAS1 variants were implicated in susceptibility to SARS-CoV in candidate gene association studies in Vietnam34 and China35. These genes encode enzymes that produce a host antiviral mediator (2′,5′-oligoadenylate (2-5A)) that activates an effector enzyme, RNase L. RNase L degrades double-stranded RNA36, a replication intermediate of coronaviruses37. The betacoronaviruses OC43 and MHV make viral phosphodiesterases that cleave 2-5A38, but SARS-CoV-2 is not known to have this ability. The OAS genes therefore also provide a potential therapeutic target: endogenous phosphodiesterase 12 (PDE-12) activity degrades 2-5A. Therapeutic PDE-12 inhibitors are available and augment OAS-mediated antiviral activity39.

The association on chromosome 19p13.3 (rs2109069, P = 3.98 × 10−12) is an intronic variant in the gene that encodes dipeptidyl peptidase 9 (DPP9). Variants in this locus are associated with idiopathic pulmonary fibrosis40. DPP9 encodes a serine protease that has diverse intracellular functions, including the cleavage of the key antiviral signalling mediator CXCL1041, and key roles in antigen presentation42 and inflammosome activation43.

As opportunities for therapeutic intervention, particularly experimental therapy, are more abundant in later, more-severe cases of disease, it is important that our results also reveal genes that may act to drive inflammatory organ injury. TYK2 is one of four gene targets for JAK inhibitors such as baricitinib44, one of the nine candidate drugs that we used in the creation of our a priori target list (Supplementary Table 3). The association between TYK2 expression and critical illness was also confirmed in an external dataset (Table 2).

We replicate the association at chromosome 3p21.31 that has previously been described8. The extremely small P value at this locus (P = 4.77 × 10−30) may reflect the large size of our study, and our focus on extreme severity, as we see a larger effect size in the GenOMICC dataset than in the replication studies (Extended Data Fig. 4). A number of genes in this locus could plausibly explain an association. Our systematic review and meta-analysis of experimental data on infections with betacoronaviruses from other sources provides moderate biological support for FYCO1, although this additional information comes mostly from in vitro model systems45. Our TWAS results show that variants in this region confer genome-wide significant differences in predicted expression of CXCR6, CCR2 and CCR3 (Fig. 2a); it is likely that at least one of these genes is an important mediator of critical illness in patients with COVID-19.

Association with critical illness for genotype-inferred CCR2 (CC-chemokine receptor 2) expression is particularly strong in lung tissue (Fig. 2b). CCR2 promotes chemotaxis of monocytes and macrophages towards sites of inflammation, and there is increased expression of the canonical ligand for CCR2 (monocyte chemoattractant protein (MCP-1)), in bronchoalveolar lavage fluid from the lungs of patients with COVID-19 during mechanical ventilation46. Concentrations of circulating MCP-1 are associated with more-severe disease47. Anti-CCR2 monoclonal antibody therapy for the treatment of rheumatoid arthritis is safe48.

The ABO locus was also previously associated with COVID-198, but did not reach genome-wide significance in the GenOMICC cohort of critically ill patients with COVID-19. Notably, there is a signal close to genome-wide significance at this locus in the combined meta-analysis (Fig. 1), suggesting that this variant may be associated with susceptibility to COVID-19, but not critical illness (Extended Data Fig. 4).

Analysis of shared heritability highlights a positive correlation with adiposity. This does not confirm a causal relationship, as a number of biases may affect this analysis. Two effects are likely to be present: first, increased body-mass index and lower socio-economic status are strong risk factors for severe COVID-1914, and second, UK Biobank participants are disproportionately drawn from social groups in which obesity is underrepresented compared with the general population49.

Because of the urgency of completing and reporting this study, we have included controls from population genetic studies with systematic differences in population structure, demographics and comorbid diseases, who were genotyped using different technologies compared with the cases that we report. Residual confounding is reflected in the genomic inflation (λ0.5) value of 1.099 for the primary analysis (Extended Data Fig. 1). We mitigated the consequent risk of false-positive associations driven by genotyping errors by genotyping the majority of our participants using two different methods (microarray and whole-genome sequencing analyses), and by verifying significant associations using two additional control groups (the 100,000 Genomes Project and Generation Scotland). The success of these mitigations is demonstrated by robust replication of our sentinel SNPs in external studies. Our meta-analysis, combining GenOMICC with multiple additional sources of genome-wide associations, has a λ0.5 = 1.017 (Extended Data Fig. 1).

There is an urgent need to extend these findings through further studies. Our MAIC results show that highly ranked genes in the GenOMICC dataset are more likely to be implicated in COVID-19 in other studies (Extended Data Fig. 2). We continue to recruit to the GenOMICC study, in the expectation that additional associations exist and can be detected with larger numbers of cases of COVID-19. Future studies using whole-genome sequencing will search the rarer end of the allele frequency spectrum for variants that increase susceptibility to COVID-19. Effect sizes are likely to be greater in the GenOMICC study because the cohort is strongly enriched for immediately life-threatening disease in patients who are either receiving invasive mechanical ventilation, or considered by the treating physicians to be at high risk of requiring mechanical support.

We have discovered new and highly plausible genetic associations with critical illness in COVID-19. Some of these associations lead directly to potential therapeutic approaches to augment interferon signalling, antagonize monocyte activation and infiltration into the lungs, or specifically target harmful inflammatory pathways. Although this adds substantially to the biological rationale that underpins specific therapeutic approaches, each treatment must be tested in large-scale clinical trials before entering clinical practice.

Methods

Data reporting

No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

Recruitment of cases

In total, 2,636 patients who were recruited to the GenOMICC study (https://genomicc.org/) had confirmed COVID-19 according to local clinical testing and were deemed, in the view of the treating clinician, to require continuous cardiorespiratory monitoring. In UK practice this kind of monitoring is undertaken in high-dependency or intensive care units. An additional 135 patients were recruited through ISARIC4C (https://isaric4c.net/)—these individuals had confirmed COVID-19 according to local clinical testing and were deemed to require hospital admission. Both studies were approved by the appropriate research ethics committees (Scotland, 15/SS/0110; England, Wales and Northern Ireland, 19/WM/0247). Current and previous versions of the study protocol are available at https://genomicc.org/protocol/. All participants gave informed consent.

Genotyping

DNA was extracted from whole blood using the Nucleon Kit (Cytiva) with the BACC3 protocol. DNA samples were resuspended in 1 ml TE buffer pH 7.5 (10 mM Tris-Cl pH 7.5, 1 mM EDTA pH 8.0). The DNA yield was measured using Qubit and normalized to 50 ng μl−1 before genotyping.

Genotyping was performed using the Illumina Global Screening Array v.3.0 + multi-disease bead chips (GSAMD-24v3-0-EA) and Infinium chemistry. In summary, this consists of three steps: (1) whole-genome amplification; (2) fragmentation followed by hybridization; and (3) single-base extension and staining. For each of the samples, 4 μl of DNA normalized to 50 ng μl−1 was used. Each sample was interrogated on the arrays against 730,059 SNPs. The arrays were imaged on an Illumina iScan platform and genotypes were called automatically using GenomeStudio Analysis software v.2.0.3, GSAMD-24v3-0-EA_20034606_A1.bpm manifest and cluster file provided by the manufacturer.

In 1,667 cases, genotypes and imputed variants were confirmed with Illumina NovaSeq 6000 whole-genome sequencing (WGS). Samples were aligned to the human reference genome hg38 and variant called to GVCF stage on the DRAGEN pipeline (software v.01.011.269.3.2.22, hardware v.01.011.269) at Genomics England. Variants were genotyped with the GATK GenotypeGVCFs tool v.4.1.8.150, filtered to minimum depth 8× (95% sensitivity for heterozygous variant detection51,) merged and annotated with allele frequency with bcftools v.1.10.2.

Quality control

Genotype calls were carefully examined within GenomeStudio using manufacturer’s and published52 recommendations, after excluding samples with low initial call rate (<90%) and reclustering the data thereafter. In brief, X and Y marker calls were all visually inspected and curated if necessary, as were those calls of autosomal markers with minor allele frequency (MAF) > 1% that had a low Gentrain score, cluster separation and an excess or deficit of heterozygous calls. Genotype-based sex determination was performed in GenomeStudio and samples were excluded if they did not match the expectation based on clinical records. Five individuals with XXY genotypes were also detected and excluded for downstream GWAS analyses. Genotypes were exported, in Genome Reference Consortium human build 37 (GRCHb37) and Illumina ‘source’ strand orientation, using the GenotypeStudio PLINK input report plugin v.2-1-4. A series of filtering steps was then applied using PLINK 1.9, leaving 2,790 individuals and 479,095 variants for further analyses (after exclusion of samples with a call rate of <95%, selection of variants with call rate of >99% and MAF > 1% and the final selection of samples using a call rate of >97%).

Kinship

Kinship and ancestry inferences were calculated following the UK Biobank49 and 1 Million Veteran Program53. First, King 2.154 was used to find duplicate individuals who have been recruited by two different routes. The analysis flagged 56 duplicate pairs, from which one was removed according to genotyping quality (GenomeStudio p50GC score and/or individual call rate). This lef a set of 2,734 unique individuals.

Regions of high linkage disquilibrium (LD) defined in the UK Biobank49 were excluded from the analysis, as well as SNPs with MAF < 1% or missingness >1%. King 2.1 was used to construct a relationship matrix up to the third degree using the King command --kinship --degree 3 and then the function largest_independent_vertex_set() from the igraph tool http://igraph.sf.net was used to create a first set of unrelated individuals. Principal component analysis (PCA) was conducted with gcta 1.955 in the set of unrelated individuals with pruned SNPs using a window of 1,000 markers, a step size of 80 markers and an r2 threshold of 0.1. SNPs with large weights in PC1, PC2 or PC3 were removed, keeping at least two-thirds of the number of pruned SNPs to keep as an input for the next round of King 2.1. The second round of King 2.1 was run using the SNPs with low weights in PC1, PC2 and PC3 to avoid overestimating kinship in individuals of non-European ancestry. After this round 2,718 individuals were considered unrelated up to the third degree.

Genetic ancestry

Unrelated individuals from the 1,000 Genome Project dataset were calculated using the same procedure as described above, and both datasets were merged using the common SNPs. The merged genotyped data was pruned with PLINK using a window of 1,000 markers, a step size of 50 and a r2 of 0.05, leaving 92,017 markers that were used to calculate the 20 first principal components with gcta 1.9. Ancestry for GenOMICC individuals was inferred using ADMIXTURE56 populations defined in the 1000 Genomes Project. When one individual had a probability >80% of pertaining to one ancestry, then the individual was assigned to this ancestry, otherwise the individual was assigned to ‘admixed’ ancestry as in the 1 Million Veteran cohort53. According to this criterion, there are 1,818 individuals of European ancestry, 190 individuals of African ancestry, 158 individuals of East Asian ancestry, 254 individuals of South Asian ancestry and 301 individuals with admixed ancestry (2 or more ancestries).

Imputation

Genotype files were converted to plus strand and SNPs with Hardy–Weinberg equilibrium (HWE) P < 10−6 were removed. Imputation was calculated using the TOPMed reference panel57 and results were given in the GRCh38 human reference genome and plus strand. The imputed dataset was filtered for monogenic and low imputation quality score (r2 < 0.4) using BCFtools 1.9. To perform GWAS analyses, files in VCF format were further filtered for r2 > 0.9 and converted to BGEN format using QCtools 1.358.

UK Biobank imputed variants with imputation score >0.9 and overlapping our set of variants (n = 5,981,137) were extracted and merged with GenOMICC data into a single BGEN file containing cases and controls using QCtools 1.3.

GWAS

Related individuals to third degree were removed. In addition, 13 individuals with American ancestry were removed as the sample size provided insufficient power to perform a reliable GWAS analyses for this group. The final dataset includes 2,244 individuals. Using PCA to infer genetic ancestry, there were 1,676 individuals of European ancestry, 149 individuals of East Asian ancestry, 237 individuals of South Asian ancestry and 182 individuals of African ancestry (Extended Data Table 1). If age or deprivation status were missing for some individuals, the value was set to the mean of their ancestry. GWAS analyses were performed separately for each ancestry group.

Tests for association between case–control status and allele dosage at individual SNPs were performed by fitting logistic regression models using PLINK59. Independent analyses were performed for each ethnic group. All models included sex, age, mean-centred age2, deprivation score decile of residential postcode, and the first 10 genomic principal components as covariates.

Genomic principal components were computed on the combined sample of all UK Biobank and GenOMICC participants. Specifically, 456,750 genetic variants were identified that were shared between the variants contained in the called genotypes in the GenOMICC dataset and imputed UK Biobank genotypes, which had an imputation info score >0.95 and MAF > 1%. After merging genotypes at these variants, variants were removed that had MAF < 2.5%, a missingness rate >1.5%, showed departure from HWE P < 10−50 or were within previously identified regions of high linkage disequilibrium within the UK Biobank. After LD pruning of the remaining variants to a maximum r2 of 0.01 based on a 1,000-variant window moving in 50 variants steps, using the PLINK indep-pairwise command and yielding 13,782 SNPs, the leading 20 genomic principal components were computed using FlashPCA260.

GWAS results for individuals of European ancestry were filtered for MAF > 0.01, HWE P > 10−50 and genotyping rate >0.99. An extra filter was added to avoid bias for using a different genotyping method and imputation panel between controls and cases. This could not be controlled for using regression because all cases and all controls were genotyped using different methods. MAFs for each ancestry were compared between UK Biobank European controls and gnomAD hg38 non-Finnish European individuals downloaded in August 202061. SNPs were removed from the GWAS results following two rules: (1) In SNPs with MAF > 10% in gnomAD, an absolute difference in MAF of 5% between gnomAD and UK Biobank controls; (2) in SNPs with MAF < 10% in gnomAD, a difference in MAF of >25% between UK Biobank controls and gnomAD. GWAS analyses of individuals of non-European ancestries were filtered for a MAF in UK Biobank controls corresponding to the same ancestry >5% and then for the SNPs that passed quality control in the European GWAS. To calculate differences between UK Biobank European individuals and gnomAD allele frequencies, gnomAD allele frequencies for individuals of non-Finnish European descent were used, as European UK Biobank controls are mainly non-Finnish.

Filtered GWAS analyses for each ancestry, containing a total of around 4.7 million SNPs, were combined in a trans-ethnic meta-analysis using METAL62 standard error mode and controlling for population stratification (genomic control on). The nearest genes were defined using the SNP2GENE function in FUMA v.1.3.663, using LD r2 > 0.6 and UK Biobank release 2 reference panel.

A sex-specific GWAS within European individuals was performed using 1,180 cases of unrelated men and 496 cases of unrelated women and 5 UK Biobank random control individuals matched by sex and ancestry for each case. Tests for association between case–control status and allele dosage at individual SNPs were performed by fitting a logistic regression model with PLINK. Age, mean age squared, deprivation decile of residential postcode and the first 10 principal components were added as covariates in the models.

Deprivation score

The UK Data Service provides measures of deprivation based on census data and generated per postcode. The latest version of the deprivation scores were published in 2017 and are based on the 2011 census. As only partial postcodes were available for most samples we were unable to use these indices directly. However, we generated an approximation to the scores by calculating an average weighted by population count across the top-level postcode areas.

The initial input file was part of the aggregated census data identified by in the 2011 Census aggregate85 and the postcode data were downloaded from http://s3-eu-west-1.amazonaws.com/statistics.digitalresources.jisc.ac.uk/dkan/files/Postcode_Counts_and_Deprivation_Ranks/postcodes.zip.

Population count and deprivation score for each published postcode were extracted and the weighted average score was calculated for each top-level postcode. We further categorized each top-level postcode score into decile and quintile bins for more coarse-grained analyses.

WGS

WGS gVCF files were obtained for the 1,667 individuals for whom we had WGS data. Variants overlapping the positions of the imputed variants were called using GATk and variants with depth <8× (the minimum depth for which 95% coverage can be expected) were filtered. Individual VCF files were joined in a multi-sample VCF file for comparison with imputed variants. We used 1,613 of the 1,667 individuals in the final GWAS. Samples were filtered and variants were annotated using bcftools 1.9. VCF files obtained from imputation were processed in an identical manner. Alternative allele frequencies were calculated with PLINK 2.064 for both WGS and imputed data.

Controls

UK Biobank

UK Biobank participants were considered as potential controls if they were not identified by the UK Biobank as outliers based on either genotyping missingness rate or heterogeneity, and their sex inferred from the genotypes matched their self-reported sex. For these individuals, information on sex (UKBID 31), age, ancestry and residential postcode deprivation score decile was computed. Specifically, age was computed as age on 1 April 2020 based on the participant’s birth month (UKBID 34) and year (UKBID 52). The first part of the residential postcode of participants was computed based on the participant’s home location (UKBID 22702 and 22704) and mapped to a deprivation score decile as described above for GenOMICC participants. Ancestry was inferred as described above for GenOMICC participants.

After excluding participants who had received PCR tests for COVID-19, based on information downloaded from the UK Biobank in August 2020, five individuals with matching inferred ancestry were sampled for each GenOMICC participant as controls. After sampling each control, individuals related up to the third degree were removed from the pool of potential further controls. An additional analysis with more stringent matching on individual characteristics was also performed (Supplementary Information, ‘Matched controls’).

The 100,000 Genomes Project

Following ethical approval (14/EE/1112 and 13/EE/032), consenting participants from the 100,000 Genomes Project with a broad range of rare diseases, cancers and infection were enrolled by 13 regional NHS Genomic Medicine Centres across England and in Northern Ireland, Scotland and Wales and whole blood was drawn for DNA extraction. After quality assurance, WGS at 125 or 150 base pairs was performed by Illumina Laboratory Services on either Hiseq 2500 or Hiseq X sequencers in the Genomics England Sequencing Centre, followed by detection of small variants (single-nucleotide variants and small insertions or deletions) using Starling.

Tests for association between cases and control status were performed by running mixed model association tests using SAIGE (v.0.39). In total, 1,675 individuals from the GenOMICC study and 45,875 unrelated participants of European ancestry were included. Genomic principal components were calculated for the combined dataset of GenOMICC participants and WGS data from the 100,000 Genomes Project. PCA was performed with GCTA software using approximately 30,000 SNPs selected with MAF > 0.005 and after LD pruning (r2 < 0.1 with a window size of 500 kb). Fitting of the null logistic mixed model was performed using the SNPs used for PCA and included age, sex, squared age, age × sex and the first 20 genomic principal components as covariates.

Tests for association using SAIGE were performed after filtering variants in the WGS dataset for genotype quality and MAF ≥ 0.05. GWAS-specific quality filtering was performed to include variants with minor allele count ≥20 for each phenotype, differential missingness between cases and controls (P < 1 × 10−5) and departure from HWE (P < 1 × 10−5).

Generation Scotland

The Generation Scotland: Scottish Family Health Study is a population-based cohort of 24,084 participants sampled from five regional centres across Scotland65. A large subset of participants were genotyped using either Illumina HumanOmniExpressExome-8v1_A or v1-2, and 20,032 passed quality control criteria that have previously been described66,67. Genotype imputation using the TOPMed reference panel was recently performed (freeze 5b) using Minimac4 v.1.0 on the University of Michigan server (https://imputationserver.sph.umich.edu)68. Imputation data from 7,689 unrelated (genomic sharing identical-by-descent estimated to be <5% using PLINK 1.9) participants were used as control genotypes in a GWAS using GenOMICC cases of European ancestry, for quality control purpose of associated variants. GWAS was performed in a logistic regression framework implemented in the glm function of PLINK 2 (https://www.cog-genomics.org/plink/2.0/), adjusting for age, sex and the first 10 principal components of European ancestry. These coordinates were obtained from projection to the principal component space of the European population samples of the 1000 Genomes Project using KING v.2.2.554 and a LD-pruned subset of target genotyped markers that passed quality control and intersecting with the reference populations.

Validation

Clumped hits in discovery GWAS analyses were validated using controls from Generation Scotland and the 100,000 Genomes Project. To consider a hit validated, the direction of effect should be the same in all three GWAS datasets and the P value in both Generation Scotland and the 100,000 Genomes Project had to be P < 0.05/nvalidations, where nvalidations is the number of significant independent loci in our analysis at the discovery threshold of P < 5 × 10−8.

Replication

Loci in GenOMICC individuals of European ancestry were defined using the clump function of PLINK 1.964 and clumping parameters r2 = 0.1, P = 5 × 10−8 and P2 = 0.01; distance to the nearest gene was calculated using ENSEMBL GRCh37 gene annotation.

No GWAS has been reported of critical illness or mortality in COVID-19. As a surrogate, to provide some replication for our findings, replication analyses were performed using HGI build 37, version 2 (July 2020) B2 (hospitalized patients with COVID-19 versus the population) GWAS. Summary statistics were used from the full analysis, including all cohorts and GWAS without UK Biobank, to avoid sample overlap. The replication P value was set to 6.25 × 10−4 (0.05/8, where 8 is the number of loci significant in the discovery analysis).

Genome-wide meta-analysis

Meta-analysis of the GenOMICC, HGI and 23andMe datasets was performed using fixed-effect inverse-variance meta-analysis in METAL62, with correction for genomic control on. The 23andMe study comprises cases and controls from an European genetic ancestry group. The HGI B2 analysis is a trans-ancestry meta-analysis, with the great majority of cases being multi-ethnic European (European and Finnish), with 238 cases of non-European ancestry (176 individuals of admixed American descent from the BRACOVID study and 62 individuals of South Asian ancestry from the GNH study).

Post-GWAS analyses

TWAS and Meta-TWAS

We performed transcriptome-wide association using the MetaXcan framework23 and the GTEx v.8 eQTL MASHR-M models (http://predictdb.org/). To increase SNP coverage to perform TWAS, GWAS summary statistics for European ancestry were first imputed using the fizi69 impute function (https://github.com/bogdanlab/fizi), the European population from the 1000 Genomes Project as LD reference and 30% as minimum proportion of SNPs for a region (--min-prop 0.3). Then, imputed GWAS results were harmonized, lifted over to hg38 and linked to the 1000 Genomes Project reference panel using GWAS tools (https://github.com/hakyimlab/summary-gwas-imputation/wiki/GWAS-Harmonization-And-Imputation).

Imputed and harmonized GWAS summary statistics were used to perform TWAS for whole-blood samples (Supplementary Fig. 16) and lung tissues (Fig. 2) in GTEx v.8 with the S-PrediXcan function. Resulting P values were corrected using the Bonferroni correction to find significant gene associations. To overcome the limitations of sample size in GTEx v.8 lung tissues and whole-blood samples, we performed a meta-TWAS prioritizing genes with small P values in these tissues and using GTEx v.8 gene expression in all tissues and S-Multixcan70.

Mendelian randomization

Mendelian randomization19 based on two-sample summary data was performed using the results of GenOMICC and GTEx v.771 (using SMR/HEIDI pre-prepared data from https://cnsgenomics.com/software/smr/#DataResource), with Generation Scotland65,72 forming a LD reference. GenOMICC results from individuals of European ancestry were used as the outcome; and GTEx (v.7) whole-blood expression results as the exposure. Additional data pertaining to GTEx v.7 were downloaded from GTEx: https://gtexportal.org/ (accessed 20 February, 5 April and 4 July 2020), and SMR/HEIDI from https://cnsgenomics.com/software/smr/ (accessed 3 July 2020). Analyses were conducted using Python v.3.7.3 and SMR/HEIDI v.1.03 (plots were made using SMR/HEIDI v.0.711). An LD reference was created using data from the population-based Generation Scotland cohort (used with permission; described previously67): from a random set of 5,000 individuals, using PLINK v.1.9 (www.cog-genomics.org/plink/1.9/), a set of individuals with a genomic relatedness cut-off of <0.01 was extracted; 2,778 individuals remained in the final set. All data used for the SMR/HEIDI analyses were limited to autosomal biallelic SNPs: 4,264,462 variants remained in the final merged dataset.

Significant (as per GTEx v.7; nominal P value below the nominal P-value threshold) local (distance to transcriptional start site <1 Mb) eQTLs from GTEx v.7 whole-blood samples for protein-coding genes (as per GENCODE v.19) with a MAF > 0.01 (GTEx v.7 and GenOMICC) were considered to be potential instrumental variables. Per variant, we first selected the Ensembl gene ID with which the eQTL was most strongly associated followed by selecting the variant with which each Ensembl gene ID was most strongly associated. Instruments were available for 4,614 unique Ensembl gene IDs.

Results were assessed based on a list of genes selected a priori as of interest (Supplementary Table 3), and together as a whole. Replication of Bonferroni-corrected significant results was attempted for the results of COVID-19 HGI (https://www.covid19hg.org/) with the UK Biobank excluded (2 July 2020 data release) using the eQTLgen expression dataset20. Hospitalized patients with COVID-19 versus population (ANA_B2_V2) was selected as the phenotype most similar to our own, and therefore the most appropriate to use as a replication cohort.

To further validate the analyses above, generalized summary-data Mendelian randomization (GSMR)73 was performed using exposure data from https://www.eqtlgen.org/index.html (accessed 26 October 2020)20 and the publicly available GenOMICC European data for TYK2 and IFNAR2 (Supplementary Fig. 15). GSMR was performed using GCTA v.1.92.1 beta6 Linux. Pleiotropic SNPs were filtered using HEIDI-outlier test (threshold = 0.01) and instrument SNPs were selected at a genome-wide significance level (PeQTL < 5 × 10−8) using LD clumping (LD r2 threshold = 0.05 and window size = 1 Mb). The imputed genotypes for 50,000 unrelated individuals (based on SNP-derived genomic relatedness <0.05 using SNPs from HapMap 3) from the UK Biobank were used as the LD reference for clumping. GSMR accounts for remaining LD not removed by LD clumping.

Genomic region plots

Genomic region plots were created using https://github.com/Geeketics/LocusZooms (Supplementary Figs. 5, 6).

Gene-level and pathway analyses

The gene-level burden of significance in the results of the European ancestry group was calculated using MAGMA v.1.0874 (Supplementary Fig. 17). SNPs were annotated to genes by mapping based on genomic location. SNPs were assigned to a gene if the SNPs location was within 5 kb up- or down-stream of the gene region (defined as the transcription start site to transcription stop site). The MAGMA SNP-wise mean method was applied, which utilizes the sum of squared SNP Z-statistics as the test statistic. The European reference panel of the 1000 Genomes Project was used to estimate LD between SNPs.

Auxiliary files were downloaded from https://ctg.cncr.nl/software/magma on 1 September 2020. Gene location files for protein-coding genes were obtained from NCBI (ftp.ncbi.nlm.nih.gov) on 29 April 2015 (gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz) and 25 May 2016 (genomes/Homo_sapiens/ARCHIVE/ANNOTATION_RELEASE.105/mapview/seq_gene.md.gz).

The reference data files used to estimate LD are derived from Phase 3 of the 1000 Genomes Project.

Competitive gene set enrichment analysis was conducted in MAGMA using a regression model that accounts for gene–gene correlations to reduce bias resulting from clustering of functionally similar genes on the genome74. Gene sets were queried from the databases KEGG 2019, Reactome 2016, GO Biological Process 2018, Biocarta 2016 and WikiPathways 2019. The Benjamini–Hochberg procedure was used to control false-discovery rate (<0.05).

MAIC

To put these results in the context of existing biological data about host genes in SARS-CoV-2 replication and response, we performed MAIC24 analysis, which integrates gene-level results from the GenOMICC metaTWAS and an existing systematic review of host factors implicated in SARS-CoV-2 viral replication and host response in COVID-1945.

We developed MAIC to evaluate and integrate gene-level data from diverse sources24. Multiple in vitro and in vivo studies have identified key host genes that either directly interact with SARS-CoV-2, or define the host response to SARS-CoV-2. We have previously conducted a systematic review of these studies45. To put the new associations from this GWAS into context, we performed a data-driven meta-analysis of gene-level results combined with pre-existing biological data using MAIC24.

In brief, MAIC aggregates both ranked and unranked lists and performs better than other methods, particularly when presented with heterogeneous source data. The input to MAIC is a list of named genes. MAIC assigns a score to each gene according to how many source datasets have reported that gene, and then creates a data-driven weighting for each data source (usually an individual experiment) based on the scores of the genes that are highly ranked on that list. This procedure is performed iteratively until the scores and weightings converge on stable values. To prevent a single type of experiment from unduly biasing the results, input gene lists are assigned to categories, and a rule applied that only one weighting from each category can contribute to the score for any given gene.

Tissue and functional genomic enrichment

We downloaded the mean gene expression data summarized from RNA sequencing by the GTEx Project (https://gtexportal.org/). The GTEx v.7 dataset contains gene expression data of 19,791 genes in 48 human tissues. Gene expression values were normalized to numbers of transcripts per million reads. To measure the expression specificity of each gene in each tissue, each gene expression specificity was defined as the proportion of its expression in each tissue among all the tissues—that is, a value ranging between 0 and 1. SNPs within the 10% most specifically expressed genes in each tissue were annotated for subsequent testing of heritability enrichment. For functional genomic enrichment analysis, we considered the inbuilt primary functional annotations v.2.2 provided in the ldsc software (https://alkesgroup.broadinstitute.org/LDSCORE/) to annotated the SNPs.

With the annotated SNPs, we used stratified LD score regression75 to test whether any human tissue or specific functional genomic feature is associated with severe COVID-19. Our GWAS summary statistics were harmonized by the munge_sumstats.py procedure in ldsc. LD scores of HapMap3 SNPs (MHC region excluded) for gene annotations in each tissue were computed using a 1-cM window. The enrichment score was defined as the proportion of heritability captured by the annotated SNPs divided by the proportion of annotated SNPs.

Genetic correlations

We applied both the LD score regression76 and high-definition likelihood25 methods to evaluate the genetic correlations between severe COVID-19 and 818 GWAS-analysed phenotypes stored on LD-Hub77. GWAS summary statistics were harmonized by the munge_sumstats.py procedure in the ldsc software. In the high-definition likelihood analysis, we estimated the SNP-based narrow-sense heritability for each phenotype, and for the 818 complex traits GWAS analyses, those with SNPs with less than 90% overlap with the high-definition likelihood reference panel were removed.

Genome build

Results are presented using Genome Reference Consortium human build 37. Imputed genotypes and WGS data were lifted over from Genome Reference Consortium Human Build 38 using Picard liftoverVCF mode from GATK 4.0, which is based on the UCSC liftover tool (chain file obtained from ftp://ftp.ensembl.org/pub/assembly_mapping/homo_sapiens/GRCh38_to_GRCh37.chain.gz)78.

Reporting summary

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