Detailed stratified GWAS analysis for severe COVID-19 in four European populations

Abstract

Given the highly variable clinical phenotype of Coronavirus disease 2019 (COVID-19), a deeper analysis of the host genetic contribution to severe COVID-19 is important to improve our understanding of underlying disease mechanisms. Here, we describe an extended genome-wide association meta-analysis of a well-characterized cohort of 3255 COVID-19 patients with respiratory failure and 12 488 population controls from Italy, Spain, Norway and Germany/Austria, including stratified analyses based on age, sex and disease severity, as well as targeted analyses of chromosome Y haplotypes, the human leukocyte antigen region and the SARS-CoV-2 peptidome. By inversion imputation, we traced a reported association at 17q21.31 to a ~0.9-Mb inversion polymorphism that creates two highly differentiated haplotypes and characterized the potential effects of the inversion in detail. Our data, together with the 5th release of summary statistics from the COVID-19 Host Genetics Initiative including non-Caucasian individuals, also identified a new locus at 19q13.33, including NAPSA, a gene which is expressed primarily in alveolar cells responsible for gas exchange in the lung.

Figure 1

Forest plot of candidates from the in-depth stratified analysis. The plots show the variants chr3:45823240, chr3:45848457:C:T, chr17:46142465:T:A which were significantly associated with age (interaction P-value [P_FDR(I) < 0.05]) when comparing age groups ≤60 and >60. We additionally show variants chr19:4717660:A:G and chr21:33242905:T:C with insignificant, though strong trends for association with age. ORs and their respective 95% confidence intervals (CIs) are visualized for each of our four cohorts separately. The size of the dots indicates the size of the respective cohort (N). The OR value is displayed in addition as a numerical value. Only cohorts in which N > 50 in both cases and controls are shown. The headers are built as follows: gene variant id (chr:pos(hg38):allele)—rs-id—effect allele ^*/**. ^* indicates a variant that was observed as associated from data in this study; ^**indicates a variant that was observed as associated in the release 5 of the COVID-19 HGI.

Figure 2

Association of the 17q21.31 locus with severe COVID-19 with respiratory failure. (A) Regional association plot showing the variant most strongly associated with severe COVID-19 (rs1819040, purple diamond) a ~ 0.9 Mb inversion polymorphism at 17q21.21 (56) (white line with blue rectangles representing the variable segmental duplication (SD) blocks at the breakpoints), and the large credible set obtained by statistical fine-mapping including 2178 SNPs in high LD (median [LD] = 0.97) with the inversion (Supplementary Material, Table S7). Pairwise LD values (r²) with lead variant rs1819040 were calculated from merged Italian, Spanish, German/Austrian and Norwegian GWAS main anlysis datasets. The dotted line indicates the genome-wide significance threshold (P = 5 × 10⁻⁸). Below, organization of the 17q21.31 inversion genomic region, with the extended haplotypes associated with each orientation (H1 and H2) shown as red and blue arrows, respectively, and breakpoint SDs as dark rectangles. Protein-coding genes for which the inversion is a lead eQTL in at least one GTEx tissue are shown as pointed rectangles indicating the direction of transcription. (B) Forest plot and extended meta-analysis of our first discovery analysis and the COVID-19 HGI release 5 analysis B2 dataset (Material and Methods) of the association between severe COVID-19 and the 17q21.31 inversion based on the presence relative to the absence of the inversion haplotype H2. We visualized the ORs and their 95% confidence invervals (CIs) across all analysed cohorts of the main analysis and the COVID-19 HGI release 5 analysis B2 data. In this analysis, the overlap between the main analysis cohort and the COVID-19 HGI data was excluded from the main analysis cohort. The OR value is displayed in addition as a numerical value. The size of the dots indicates the size of the respective cohort (N). (C) Phenome-wide association study (PheWAS) results for the 17q21.31 inversion allele H2 showing only potentially COVID-19-related phenotypes from the GWAS Catalog (P = 10⁻⁷) grouped by disease categories using different colors. The effect direction of known SNP-trait associations from the corresponding GWAS is shown using triangles pointing upward (increase) and downward (decrease), whereas dots represent unknown effect direction. Phenotypes shown were selected according to previously reported COVID-19 links with lung damage, blood cell alterations and exacerbated immune response, as well as some potential co-morbidities. The whole list of phenotypic associations is included in Supplementary Material, Table S11.

Figure 3

Expression analysis of the most plausible candidate genes associated with the 17q.21.31 and 19q.13.33 loci in organ tissues and COVID-19 relevant cell types. (A) GTEx tissue-specific expression QTL (eQTL, upper panel) and splicing QTL (sQTL, middle panel) effects of the 17q21.31 and 19q.13.33 loci on selected candidate genes as well as expression of these genes in GTEx (20) tissues (lower panel). The direction of the normalized eQTL and sQTL effect size (NES) of the lead SNP rs1405655 and the inversion tagger rs62055540 in perfect LD with the inversion is represented by color intensities, and statistical significance by dot size. Black rectangles indicate genes for which the expression colocalizes (regional probability > 0.9) with GWAS loci in a given human tissue from the GTEx dataset. Heatmap displays gene-wise centered median by tissue expression values (represented by color intensities), showing in which tissues candidate genes are mostly enriched. (B) Expression levels of candidate genes in scRNA-seq datasets from healthy upper airways (nasal, bronchi) and lung (parenchyma) cells (47) and adult human brain cells from recently deceased, non-diseased donors (48). Figure displays log-normalized mean expression (represented by color) and fraction of cells expressing those genes (indicated by dot size). Processed and cell-type-annotated gene expression levels from studies were retrieved from COVID-19 Cell Atlas (49). (C) The figure shows differential expression of candidate genes in lung cells of COVID-19 patients compared with healthy controls. Log₂ fold change (log2FC) values are presented as color gradient. Nominal P-values in −log₁₀ scale are shown proportionally to dot size. Black-bordered circles indicate significantly differentially expressed genes after FDR correction. Results were obtained from pseudo-bulk differential expression analysis by Delorey et al. (21). More detailed figures are shown in Supplementary Material, Figures S11, S12 and S15.