key: cord-0805977-9xhbofva
authors: Moon, Chang; Schilder, Brian M.; Raj, Towfique; Huang, Kuan-lin
title: Phenome-wide and eQTL Associations of COVID-19 Genetic Risk Loci
date: 2021-05-18
journal: iScience
DOI: 10.1016/j.isci.2021.102550
sha: 2dd3346332057102c70a591eb0c61aafdb809e9a
doc_id: 805977
cord_uid: 9xhbofva

While several genes and clinical traits have been associated with higher risk of severe COVID-19, how host genetic variants may interact with these parameters and contribute to severe disease is still unclear. Herein, we performed phenome-wide association study (PheWAS), tissue and immune-cell specific expression/splicing quantitative trait loci (eQTL/sQTL), and colocalization analyses for genetic risk loci suggestively associated with severe COVID-19 with respiratory failure. Thirteen phenotypes/traits were associated with the severe COVID-19-associated loci at the genome-wide significance threshold, including monocyte counts, fat metabolism traits, and fibrotic idiopathic interstitial pneumonia. In addition, we identified tissue and immune subtype-specific eQTL associations affecting 48 genes, including several that may directly impact host immune responses, colocalized with the severe COVID-19 GWAS associations, and showed altered expression in single-cell transcriptomes. Collectively, our work demonstrates that host genetic variations associated with multiple genes and traits show genetic pleiotropy with severe COVID-19 and may inform disease etiology.

Upon emerging from Wuhan, China in December 2019, the outbreak of Coronavirus disease 2019 , caused by Severe Acute Respiratory Syndrome-Coronavirus-2 (SARS-CoV-2), has become a global pandemic, infecting, as of August 2020, more than 24 million people worldwide While most infected individuals either are asymptomatic or present with mild flulike symptoms (i.e., fever, cough, sore throat, malaise, etc.), few experience severe manifestations of COVID-19 that span multiple organ systems Guan et al., 2020; Huang et al., 2020) .

Although several clinical observational studies have identified several demographic and clinical factors associated with increased disease severity and mortality such as sex, age and other comorbidities such as cancer, chronic kidney disease and obesity, which and how host genetic factors may contribute to the pathogenesis and pathophysiology of severe COVID-19 are poorly understood (Williamson et al., 2020) While a recent effort by the Severe COVID-19 Genome Wide Association Study (GWAS) Group identified two significant loci (3p21.31 and 9q34.2) and multiple suggestive host genetic loci associated with increased risk of suffering from severe COVID-19 with respiratory failure, the functional effects of these loci and how they may influence disease severity remains unclear (Ellinghaus et al., 2020) .

Phenome-Wide Association Study (PheWAS) is a genotype-to-phenotype approach to identify associations between a particular genetic variant and a wide range of human traits (Hebbring, 2014) .

Since the first published PheWAS in 2010 (Denny et al., 2010) , PheWAS has been used to interrogate genetic pleiotropy, facilitate drug discovery and repositioning, and discover novel genetic etiologies of complex and polygenic disorders (Diogo et al., 2018; Hebbring, 2014; Semmes et al., 2020; Zhao et al., 2018) , and systematic compilations of GWAS results, such as that recently undertaken by GWAS Atlas, have facilitated PheWAS (Watanabe et al., 2019) . Considering the multi-organ involvement and manifold individual variations of severe COVID-19, it's likely that host susceptibility to severe COVID-19 is determined by multiple genetic factors. By connecting variants associated with both severe COVID-19 and other diseases or phenotypes, PheWAS analysis can use genetic pleiotropy to validate known comorbidities and unveil novel biological insights. Furthermore, PheWAS may offer new insights into the genetic etiology of severe COVID-19 by leveraging known biology behind seemingly disparate phenotypes.

In addition to PheWAS that can link categorical phenotypes, tissue-specific quantitative trait loci (QTL) analysis can help to prioritize gene-disease connections for mechanistic investigations. With the advent of immense functional genomic databases such as GTEx (Consortium, 2020) , the J o u r n a l P r e -p r o o f transcriptional effects of genetic variants can be investigated at tissue-specific resolution. The findings gathered from tissue-specific quantitative trait locus (QTL) analysis to severe can help to 1) identify candidate genes that may contribute to the pathophysiology of severe 2) focus the range of mechanistic hypotheses regarding the effect of specific gene expressions to organs most affected by severe COVID-19, 3) cross-reference with known biology of linked phenotypes uncovered by the PheWAS analysis to validate and refine hypotheses regarding the pathophysiological mechanisms of severe COVID-19.

In this present study, we conducted both PheWAS and QTL analyses to identify potential mechanisms through which host genomic factors contribute to the pathogenesis and pathophysiology of severe COVID-19. We extracted the lead Single Nucleotide Polymorphism (SNPs) with the smallest p-value from each locus showing suggestive associations with severe COVID-19 cases with respiratory failure using the recent GWAS (Ellinghaus et al., 2020) and

conducted PheWAS using these lead SNPs by leveraging the GWAS Atlas database and identified a list of linked phenotypes. We systematically performed colocalization analyses between each locus identified by the PheWAS and 112 QTL datasets from GTEx (Consortium, 2020) and the eQTL Catalogue (Kerimov et al., 2020) , spanning 69 unique tissue types. This provided a robust list of candidate genes with their tissue-specific expression alterations. Together, the results presented here have generated multiple focused hypotheses regarding the pathophysiology of severe COVID-19 for future validation and investigation.

J o u r n a l P r e -p r o o f

To identify phenotypes sharing genetic etiologies with severe COVID-19, we conducted PheWAS of 22 pruned genetic risk loci significantly or suggestively associated (P < 1e-5) with severe COVID-19 based on GWAS results published by The Severe Covid-19 GWAS group (Accession numbers: GCST90000256) (Ellinghaus et al., 2020) SNPs, including rs2440652, rs3934992, rs12610495, rs11385942, and rs657152, showed an association passing the Bonferroni threshold (P < 0.05/22), and were further considered as "COVID-19 hg-replicated SNPs" in the following texts. To gauge the independence of the signal for each lead SNP, linkage disequilibrium (LD) analysis was performed using the LDmatrix tool (Machiela and Chanock, 2015) . LD analysis demonstrated that the 22 lead SNPs were independent from each other (R 2 < 0.01; Supplementary Table 2 ) and each replicated SNP was found within a LD block (R 2 > 0.8) consisting of its neighboring variants.

Upon conducting PheWAS of the 22 lead SNPs using the GWAS ATLAS (Watanabe et al., 2019) that has aggregated summary statistics of over 4,756 GWAS studies, we identified 5 SNPs, including the 2 reported COVID-19 GWAS SNPs (rs647800, rs11385942) (Ellinghaus et al., 2020) and two COVID-19 hg-replicated SNPs (rs3934992, rs12610495) , that showed associated traits reaching the genome-wide significance threshold (P < 5e-8) (a full list of PheWAS results were listed in Supplementary Table 3) . A total of 13 traits were significant ( Figure 1A )-8 of which were associated with rs647800, an intronic variant of the ABO gene (Yeung et al., 2017) . Significant associated traits were predominantly immunological or metabolic ( Figure 1B) . The immunological traits were mostly relevant to population counts/proportions of monocytes and erythrocytes. Top metabolic traits were biochemical and physiological metrics concerning coagulation ("Activated partial thromboplastin time") and fat metabolism ("total cholesterol", "legs-leg fat ratio (male)). Of all the significantly associated traits, fibrotic idiopathic interstitial pneumonia was the only disease linked through the COVID-19 GWAS lead SNPs. Taken together, PheWAS analysis of 22 lead SNPs identified several physiological metrics as well as a disease associated with severe COVID-19.

To evaluate the functional relevance of the lead GWAS loci, we used the 22 lead SNPs to query GTEx and compiled e/sQTLs associated with the lead SNPs. After filtering associations based on QTL FDR < 0.05, 11 lead GWAS SNPs were identified as significant e/sQTLs (Supplementary Table 4 ,5) with the vast majority of eQTLs (n = 227/258) and sQTLs (n = 140/143) being associated with one SNP (rs3934992; Figure 2A ). No one particular tissue type dominated the eQTL and sQTL associations; the top tissue types with associated eQTLs include the muscularis layer of esophagus (n = 10), tibial nerve (n = 10), sun-exposed skin at lower leg (n = 10) and lung (n=6). The eQTL analysis encompassing 22 lead SNPs identified 29 novel genes in addition to 4 genes (i.e., SLC6A20, FYCO1, CXCR6, CCR1) previously reported by The Severe COVID-19 GWAS Group to be associated with rs11385942 at locus 3p21.31 (Ellinghaus et al., 2020) . Two sets of three genes (CMB9-55F22.1, PIDD1, RPLP2), all of which are associated with the predominant lead COVID-19 hg-replicated SNP (rs3934992), constituted top 10% of statistically significant eGenes and sGenes, respectively.

In addition to these top genes, we also identified several genes with known immune and metabolic functions-paralleling our findings in the PheWAS analysis ( Figure 2B ). Corresponding to the top metabolic PheWAS traits related to fat metabolism, one distinct COVID-19 hg-replicated SNP was a genome-wide significant e/sQTL of PNPLA2, a critical enzyme involved in the first step of triglyceride hydrolysis (Taxiarchis et al., 2019) . In terms of immunity, expression changes in CD151, PIDD1 and DPP9-genes involved in immune cell activation and NF-kB-mediated immune response (Okondo et al., 2017; Rocha-Perugini et al., 2014; Tinel et al., 2007) -were associated with two distinct COVID-19 hg-replicated SNPs (rs3934992, rs12610495). Overall, the QTL analyses of the lead SNPs identified genes of which expression changes may be functionally associated with severe COVID-19 at a tissue-specific resolution.

Based on the identification of genes with known immune function in the eQTL analysis, as well as the reported involvement of the dysregulated immune response in severe COVID-19 (Lucas et al., 2020) , we further queried the 22 lead SNPs in the eQTL Catalogue to determine whether they affected gene expression in specific immune cell-subtypes. We identified 134 immune cell subtypespecific eQTL associations (FDR <0.05) that are distributed across 7 GWAS SNPs and 31 eGenes.

Several genes that were the most statistically significant in the GTEx analysis (RPLP2, PIDD1) reemerged in the immune subtype-specific eQTL analysis and were primarily seen in monocytes and J o u r n a l P r e -p r o o f macrophages (Mo/M) and several T cell subsets. Moreover, expression changes of a few immune genes (DPP9, CD151, PIDD1) that were previously highlighted were predominantly seen in both CD4 + and CD8 + T cells as well as Mo/M. Notably, 14 of these genes were not nominated by the GTEx eQTL analysis but were associated in various T, B, NK and myeloid cell subsets (Figure 3) , including several genes involved in leukocyte recruitment and migration such as CCR1, CCR2, CCR3, CCR5, CXCR6. These results highlight the resolution gained by conducting cell type-specific analyses into immune cell-expressed genes that may be linked to the genetic factors underlying severe COVID-19.

Due to linkage disequilibrium, overlapping significant SNPs between GWAS and QTL do not necessarily mean they share the same underlying causal signal. Nor does a lack of overlapping lead variants necessarily mean that they do not share the same signal being tagged by different variants.

To robustly test for GWAS-eQTL signal sharing, we performed a series of colocalization analyses (Wakefield, 2009) and PIDD1 are perhaps controlled by a common regulatory element that is disrupted by the severe COVID-19 risk signal in this region.

While locus chr1-161253626-C-T was also colocalized with eQTLs in monocytes (naïve), the signal modulated expression of CD48 and Signaling Lymphocytic Activation Molecule Family Member 1 (SLAMF1), both of which are immunoglobulin-like receptors of the CD2 subfamily ( Figure 4B ). Interestingly, a closer inspection revealed that the colocalizing GWAS signal was not the same one tagged by the previously identified proxy SNP showing the most significant association with severe COVID-19 (rs114173811), but instead a signal ~0.5Mb upstream tagged by J o u r n a l P r e -p r o o f the SNP rs1980606. LD analysis confirmed that these signals are independent (no SNPs in the upstream signal had r 2 >0 with the lead SNP in the downstream signal). This highlights the complexity of genetic association signals as well as the additional insights that colocalization analyses can reveal.

To determine if expression of any of these eQTL-nominated immune genes as well as genes involved in implicated pathways (i.e. inflammasome activation) was altered in severe COVID-19 patients compared to moderate patients or healthy controls, we performed a differential gene expression (DGE) analysis on a recently-published scRNA-seq dataset of bronchoalveolar lavage fluid of COVID-19 patients of varying disease severity and healthy controls (Liao et al., 2020b) .

Because the functions of these immune genes are either most well-characterized in Mo/Ms (DPP9, NLRP1, GSDMD) (Okondo et al., 2017) or known to be most relevant for Mo/Ms (CCR1, CCR2, CCR5) (Dyer et al., 2019) , we specifically focused our DGE analysis on Mo/Ms subsets within the dataset and evaluated the differential expression of six eQTL-nominated immune genes. Average expression of genes involved in inflammasome activation (DPP9, NLRP1, GSDMD) was significantly altered (adj. P < 0.001) between Mo/Ms of severe and moderate COVID-19 patients.

Moreover, average expression of several genes involved in leukocyte recruitment and migration (CCR1, CCR2, CCR5) was significantly higher (Bonferroni adj. P < 0.001) in Mo/Ms of severe COVID-19 patients when compared to both moderate COVID-19 patients as well as healthy controls (Supplementary Figure 3) . These results demonstrate in at least a subset of COVID-19 patients the differential expression of eQTL-nominated immune genes between Mo/Ms of severe COVID-19 patients and moderate patients or healthy controls.

In this study, we characterized the phenome-and transcriptome-wide associations of genetic loci implicated in severe COVID-19 with respiratory failure. We leveraged GWAS ATLAS database (Watanabe et al., 2019) containing summary statistics from 4,756 GWAS studies to perform a phenome-wide associated study of the lead SNPs and identified 13 significant traits linked to genetic risk loci of severe COVID-19. Using both eQTL databases, we performed eQTL analysis of J o u r n a l P r e -p r o o f the lead SNPs to characterize their effects on gene expression in tissue-and immune cell subsetspecific resolution. 18 lead SNPs were identified as eQTLs/sQTLs of 48 genes. Taken together, these results reveal insights into potential mechanisms that could drive the formation of focused hypotheses about the pathophysiology of severe COVID-19.

Immunological and metabolic traits were amongst significant phenotypes uncovered by the PheWAS analysis. Severe cases of COVID-19 are often characterized by dysregulated and hyperinflammatory immune system (Lucas et al., 2020; Merad and Martin, 2020; Vabret et al., 2020) and metabolic dysfunction (Ayres, 2020) -suggesting that these SNPs may either play a role in multisystem manifestations of severe COVID-19 or contribute to the development of comorbidities known to exacerbate the severity of disease. For instance, that top immunological traits were hematologic parameters such as monocyte and erythrocyte populations is concordant with the hematological findings in severe COVID-19 cases such as monocyte expansion (Merad and Martin, 2020; Wen et al., 2020) and decrease in hemoglobin levels (Chowdhury and Anwar, 2020; Li et al., 2020) .

Recent reports have suggested that cytokine profiles of severe COVID-19 patients show similarities to cytokine release syndromes such as macrophage activation syndromes (Giamarellos-Bourboulis et al., 2020; Mehta et al., 2020; Moore and June, 2020) and that monocytes/macrophages extensively accumulate in the lungs (Sánchez-Cerrillo et al., 2020) . Further investigations are required to disseminate how these lead SNPs (rs647800 and rs11385492) may similarly or differentially contribute to monocyte count and monocyte-driven inflammation seen in COVID-19.

Coupled with hyperinflammation driven partly by these mononuclear cells, another major prong of severe COVID-19 pathology is systemic coagulopathy. Patients suffering from severe COVID-19 often have lower levels of platelet counts and increased levels of D-dimers (Becker, 2020; Chan and Weitz, 2020; Iba et al., 2020; Liao et al., 2020a) . Numerous cases of widespread pulmonary thrombosis and disseminated intravascular coagulation (DIC) or sepsis-induced coagulopathy (SIC) have also been described in severe COVID-19 patients (Becker, 2020; Tang et al., 2020) . Possibly reflective of the predominant role of coagulopathies in the pathophysiology of severe COVID-19, the most statistically significant trait identified in our PheWAS analysis was activated partial thrombin time (p value = 2.29e-43; associated SNP: rs647800) which has been reported to be prolonged in patients with severe COVID-19. Notably, a recent review by Merad et al (Merad and Martin, 2020) has suggested a link between hyperinflammatory monocytes and coagulopathies-as even in the absence of any vascular injury, coagulation pathway can still be activated through the recruitment of tissue factor (TF)-expressing inflammatory monocytes. Taken together, it would be interesting to J o u r n a l P r e -p r o o f interrogate the genetic pleiotropy underlying rs647800 and how it affects monocyte count, activated partial thrombin time, and severe COVID-19.

In addition to coagulation, another metabolic traits highlighted by the PheWAS results were regarding cholesterol levels ("total cholesterol"; p-value = 3.29e-10) and physiological measures associated with obesity ("waist-hip ratio" and "legs-leg fat ratio"; p-values = 1.34e-10 and 1.68e-8, respectively). Both obesity and high cholesterol levels have been associated with increased risk of severe COVID-19 Zhu et al., 2020) , thereby suggesting that the PheWAS analysis can capture some of the comorbidities found by epidemiological studies.

The PheWAS also identified rs12610495 to be associated with fibrotic idiopathic interstitial pneumonias (FIIP; p-value: 9.57e-9), a disease directly connected to respiratory failure. FIIP shares a similar risk profile with severe COVID-19 disease and thus is a potential comorbidity of severe COVID-19 (George et al., 2020) . Notably, in addition to sharing a similar risk profile with FIIP, some patients with severe COVID-19 suffer from lung function abnormalities that strongly suggest the presence of pulmonary fibrosis (Raghu and Wilson, 2020) , implicating similar etiologies underlying COVID-19 and FIIP.

In line with the immunological and metabolic traits that emerged from the PheWAS analysis, our QTL analyses identified several genes with related functions that contribute to those traits. For instance, we found that rs12610495, an intron variant of DPP9, is within a susceptibility locus for FIIP (Fingerlin et al., 2013) and modulates both the expression and splicing of DPP9. DPP9 can inhibit NLRP1 and repress downstream inflammasome activation (Zhong et al., 2018) . Conversely, pharmacological inhibition of DPP9 has been shown to induce NLRP1 inflammasome-mediated pyroptosis in monocytes and macrophages (Gai et al., 2019; Okondo et al., 2017; Vasconcelos et al., 2019; Zhong et al., 2018) . NLRP1-dependent pyroptosis in the lung alveolar macrophages has been shown to cause progressive lung damage that is similar to one seen in acute respiratory distress syndrome (ARDS) (Kovarova et al., 2012) . Furthermore, NLRP1-dependent pyroptosis is sufficient to cause progressive lung injury, independent of NLRP3 activation or IL-1b secretion, which suggests a relevant role of DPP9 in mediating lung injury seen in severe COVID-19 (Kovarova et al., 2012) . Lastly, it was recently shown that human NLRP1 senses double-stranded RNA (dsRNA) and

activates an inflammasome complex upon dsRNA sensing (Bauernfried et al., 2021 )-thereby raising a distinct possibility that dsRNA replication intermediates of SARS-CoV2 (V'kovski et al.,

activation. Considering the evidence together, it is possible that decreased DPP9 expression in the J o u r n a l P r e -p r o o f lung can predispose COVID-19 patients to develop lung injury mediated by increased NLRP1dependent pyroptosis. Concordant with this hypothesis, the lead SNP rs12610495 is associated not only with decreased expression of DPP9 specifically in the lung and in macrophages but also with increased risk of severe COVID-19 with respiratory failure. Furthermore, our DEG analysis of scRNA-seq of bronchoalveolar fluid of COVID-19 patients showed that expression of DPP9, NLRP1 as well as GSDMD, an effector protein for pyroptosis (Mathur et al., 2018) , was significantly altered in Mo/Ms of severe COVID-19 patients compared to their moderate counterparts. Overall, our results suggest the potential role of inflammasome activation and pyroptosis in exacerbating respiratory distress in severe COVID-19-a point that was recently echoed in a recent review of inflammasomes in COVID-19 (Yap et al., 2020) .

Beyond DPP9, our QTL analyses identified several immune genes that may play a role in potential mechanisms through which host immune response against SARS-CoV-2 can be affected.

These potential mechanisms include NF-kB-mediated immune response (i.e. PIDD1) (Bock et al., 2013; Tinel et al., 2007) , immune synaptic interaction involved in T cell activation (i.e. CD151) (Rocha-Perugini et al., 2014) and T cell exhaustion in viral infection (i.e. CD244) (Agresta et al., 2018; Pacheco et al., 2013) . In our immune subset-specific QTL analysis, we were able to resolve changes in expression of PIDD1 and CD151 in diverse T cell subsets as well as monocytes. Moreover, the QTL overlap and colocalization analyses showed that the genetic signal associated with increased severe COVID-19 risk also modulated immune genes (i.e. PDDC1, CD150, CD48) in monocytes.

For instance, one of the colocalized genes, CD48, is known to interact with CD244 to modulate natural killer (NK) cell and CD8 + T cell effector function (Agresta et al., 2018; McArdel et al., 2016) .

Taken together, both QTL overlap and colocalization analyses highlight strong associations between risk of severe COVID-19 and transcriptional changes that underpin monocyte-lymphocyte interactions.

Furthermore, we also identified several leukocyte adhesion and chemokine receptors (CCR3, CCR5) in addition to the chemokine receptors (CCR1, CCR2, CXCR6) previously reported by the Severe COVID-19 GWAS Group (Ellinghaus et al., 2020 )-all of which can mediate the increased leukocyte recruitment seen in the lungs of severe COVID-19 patients (Sánchez-Cerrillo et al., 2020) In line with our immune subset-specific QTL analysis and the aforementioned histological findings, we observed in our DEG analysis a significant increase in the expression of these chemokine receptors (CCR1, CCR2, CCR5) in Mo/Ms of severe COVID-19 patients compared to those of both moderate patients and healthy controls. While current scRNA-seq dataset did not include genotype J o u r n a l P r e -p r o o f information, future single-cell datasets with paired DNA/RNA data can help confirm the link between genetic risk allele and their eQTL effects at a single-cell resolution.

Concordant with several metabolic traits identified in the PheWAS such as obesity and cholesterol levels (Kanai et al., 2018; Pulit et al., 2019; Rask-Andersen et al., 2019) , the QTL analyses identified two genes, PNPLA2 and USF1, whose functions are directly connected to those traits.

Furthermore, we were able to resolve changes in expression of USF1 to the level of various immune cell compartments and found that USF1 expression was altered predominantly in Mo/Ms and T cell subtypes. USF1 is a transcription factor that modulates several genes involved in lipid and glucose metabolism such as apolipoproteins and lipases including PNPLA2, a triglyceride lipase that catalyzes the initial step in triglyceride hydrolysis (Aguilar-Salinas et al., 2014; Taghizadeh et al., 2019) .

Knock-out studies of USF1 in mice have been shown to enhance cholesterol efflux in macrophages and downregulate secretion of pro-inflammatory cytokines such as MCP-1 and IL-1 (Ruuth et al., 2018) and USF1 variants in humans have been associated with familial combined hyperlipidemia and coronary artery disease (Fan et al., 2014; Laurila et al., 2016; Ruuth et al., 2018) . Interestingly, a recent paper by Dias et al showed that monocytes from COVID-19 have increased lipid droplet accumulation and that pharmacological inhibition of lipid droplet formation led to a significant inhibition of replication of SARS-CoV-2 (Dias et al., 2020) . In addition to lipid metabolism, kynurenine metabolism has also been implicated by one of our nominated genes, KYAT3, an aminotransferase that transaminates kynurenine to kynurenic acid. Considering that kynurenine pathway plays an integral role in immunosuppression (Routy et al., 2016; Wirthgen et al., 2018; Zhai et al., 2020) , it would be interesting to investigate how expression changes in KYAT3 can affect the immune response against SARS-CoV2. Notably, alterations in both metabolic pathways implicated in the QTL analyses-fatty acid and kynurenine metabolism-have also been observed in COVID-19 patients (Thomas et al., 2020) , further supporting that these genes may play a role in immune dysfunction seen in severe COVID-19.

In summary, we conducted PheWAS and QTL analyses to investigate how host genetic factors may contribute to pathophysiology of severe COVID-19 and uncovered numerous gene-diseasetissue connections. Our results provide plausible consequences of the host genetic risk factors underlying severe disease, informing future investigations into immune, metabolic, and other aspects of COVID-19 pathophysiology.

We acknowledge the limitations of our approaches, such as the restriction of discovering traits studied by available GWAS or QTL datasets and their correlational nature, which may be addressed by techniques such as Mendelian randomization as COVID-19 GWAS expands in size (T.C.-19 H.G., 2020). Widespread linkage disequilibrium also poses challenges when trying to distinguish causal variants from their close correlates-an issue which fine-mapping techniques aim to resolve (Broekema et al., 2020; Hutchinson et al., 2020; Schaid et al., 2018) . Furthermore, the suggestive genetic risk loci of severe COVID-19 used for analyses herein as well as the identified PheWAS and QTL associations also require validation. 

The authors declare no competing interests. The eQTLs were compiled from the eQTL Catalogue (Kerimov et al., 2020) 

Information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Kuan-lin Huang (kuan-lin.huang@mssm.edu).

This study did not generate new unique reagents.

All relevant data are available from the Lead Contact upon request. All code use in the eQTL Catalogue and colocalization analyses can be found in the following GitHub repository: https://github.com/RajLabMSSM/COVID19_PheWAS

Additional Supplemental Items are available from Mendeley Data at http://dx.doi.org/0.17632/nyn84yyms7.1 and https://data.mendeley.com/datasets/nyn84yyms7/1

Extraction of severe COVID-19 associated lead SNPs. We downloaded the summary statistics of the genome-wide meta-analysis that was corrected for age and sex as reported by The Severe

Covid-19 GWAS group (Accession numbers: GCST90000256) (Ellinghaus et al., 2020) . Aside from J o u r n a l P r e -p r o o f the reported SNPs, namely rs11385942 at locus 3p21.31 and with rs657152 at locus 9q34.2, we extracted the SNPs with suggestive associations (P<1e-5), resulting in 319 SNPs from GCST90000256. We then pruned the associated SNPs based on the most significant P value per locus (±100Kb regions surrounding the lead SNP) recursively within the study, resulting in 22 nonadjacent loci, including 2 genome-wide significant SNPs (P < 5e-8) (Altshuler and Donnelly, 2005; Pe'er et al., 2008) and 20 suggestive SNPs associated with severe COVID-19 with respiratory failure (Supplementary Table 1 ). We then cross-referenced 22 lead SNPs against the meta-analysis of GWAS studies between hospitalized COVID-19 patients and the population controls (B2_ALL, dated 09/30) performed by the COVID-19 Host Genetics Initiative (COVID-19 hg) (T.C.-19 H.G., 2020) We considered these SNPs' association as "replicated" in the independent GWAS if they surpassed the Bonferroni-corrected p-value for the 22 queried SNPs (P < 0.05/22), validating associations for 3 additional SNPs (rs2440652, rs3934992, rs12610495) in addition to the two previously reported COVID-19 GWAS SNPs. In order to evaluate the level of independence of each SNP, linkage disequilibrium (LD) analysis was performed using the LDmatrix tool (Machiela and Chanock, 2015) (1) between any pair of the 22 lead SNPs that are located on the same chromosome ( (Supplementary Table 2) and (2) PheWAS analysis. Phenotype associations of the 22 lead SNPs were queried using the PheWAS database available at the GWAS ATLAS (Watanabe et al., 2019) and the PheWAS results of each SNP was compiled into the candidate list of PheWAS hits (Supplementary Table 3 ). Upon compilation, the list of associations was filtered using the genome-wide significance threshold (P < 5e-8) (Altshuler and Donnelly, 2005; Pe'er et al., 2008) , resulting in 17 significant phenotypes. Considering the varying statistical thresholds used in the PheWAS literature (Denny et al., 2013; Diogo et al., 2018; Hebbring, 2014) , we focused on a concise description of the associations passing the genome-wide significance threshold in the main texts, while noting that the full list of PheWAS results are available in Supplementary Table 3 enabling examination of other associations.

Tissue-and immune-specific eQTL/sQTL analysis. Similar to the PheWAS analysis, we used our list of 22 lead SNPS to query tissue-specific expression QTL (eQTL) and splicing QTL (sQTL) datasets from the GTEx database (Consortium, 2020) . In order to correct for multiple testing, e/sQTL associations were then filtered using the Benjamini-Hochberg(Benjamini and Hochberg, 1995) , resulting in 258 eQTLs and 143 sQTLs with FDR<0.05. In addition to the GTEx analysis, we used the R package catalogueR (https://github.com/RajLabMSSM/catalogueR) (Schilder et al., 2020) to query these same 22 SNPs in the eQTL Catalogue (Kerimov et al., 2020) , which contains 63 eQTL datasets across 20 studies (excluding GTEx). GWAS hits within sex chromosomes were excluded as the eQTL Catalogue only contains results for autosomes. eQTL Catalogue contains many immune subtype-specific and stimulation condition-specific datasets, permitting greater resolution of our biological interpretations. 190 significant eVariants associated with 9 lead GWAS SNPs (Supplementary Table 6 of which were immune subtype-or state-specific.

Using the summary statistics of the original GWAS study (GCST90000256) and the eQTL Catalogue datasets, colocalization analysis was performed on 9 lead SNPs whose eVariants identified in the eQTL Catalogue were statistically significant (FDR<0.05).

The locus for each lead SNP was defined as ± 1Mb from the lead SNP (2Mb windows). Since the original GWAS summary statistics do not provide both minor allele frequency (MAF) or the sample size for each SNP, MAF of each SNP was assumed to be the same as the one provided by the eQTL Catalogue. Each SNP was assumed to have the same effective sample size calculated from the total number of participants in the GWAS while accounting for the proportion of cases/controls using the formula:

A series of Bayesian colocalization tests of were run to evaluate the probability of shared genetic signal between each GWAS locus and each eQTL locus. Specifically, we used the run_coloc() function within cataloguer (Schilder et al., 2020) , which runs Approximate Bayes Factor colocalization (Wakefield, 2009) to first infer the probability that each SNP is causal for the phenotype within each respective dataset and then calculate the probability that the same underlying genetic signal is shared J o u r n a l P r e -p r o o f in both datasets. eQTL pseudogenes were identified via Ensembl and removed from downstream analysis. For each GWAS-QTL locus comparison, the signals were considered to be colocalized if they met the following criterion (where PP.H4 is the posterior probability that the signals are associated with their respective traits and shared, and PP.H3 is the posterior probability that the signals are associated with their respective traits but not shared): controls was downloaded and analyzed for differential expression of six immune genes nominated from the immune-specific eQTL analyses above (DPP9, NLRP1, GSDMD, CCR1, CCR2, CCR5) in monocyte/macrophage subsets using the analytical pipeline described previously (Liao et al., 2020b) Briefly, using the cluster annotations described in the original publication to isolate monocyte/macrophage subsets, the macrophage clusters were then integrated using IntegrateData() function with default parameters from the Seurat v3 package (Stuart et al., 2019) and differential gene expression analysis between disease severity groups was performed using MAST (Finak et al., 2015) in Seurat v3 (Stuart et al., 2019) to generate a list of differentially expressed genes between two groups (i.e. Severe vs. Moderate, Moderate vs. Healthy, Severe vs. Healthy). The list was then queried for specific immune genes of interest for their statistics. In accordance with Seurat's best practices, Bonferroni adjusted p-values were computed taking into account the full number of genes tested in the dataset (Stuart et al., 2019) . Violin plots comparing their expression levels as well as the UMAP plot of the integrated macrophage subset are shown in Supplementary Figure 3 .

All quantification and statistical analyses were performed as described in the Methods Details section of the STAR Methods. 

A C 0 0 5 5 9 4 . 3 A C 0 0 7 2 9 2 . 4 A P 0 0 6 6 2 1 . 5 A P 0 0 6 6 2 1 . 6 C 1 1 o r f 3 5 

The Emerging Role of CD244 Signaling in Immune Cells of the Tumor Microenvironment

Genetic and environmental determinants of the susceptibility of Amerindian derived populations for having hypertriglyceridemia

A haplotype map of the human genome

A metabolic handbook for the COVID-19 pandemic

Human NLRP1 is a sensor for double-stranded RNA

COVID-19 update: Covid-19-associated coagulopathy

Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing

Loss of PIDD limits NF-κB activation and cytokine production but not cell survival or transformation after DNA damage

A practical view of fine-mapping and gene prioritization in the post-genome-wide association era

COVID-19 coagulopathy, thrombosis, and bleeding

Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study

Management of Hemoglobin Disorders During the COVID-19 Pandemic

The GTEx Consortium atlas of genetic regulatory effects across human tissues

PheWAS: demonstrating the feasibility of a phenome-wide scan to discover gene-disease associations

Systematic comparison of phenome-wide association study of electronic medical record data and genome-wide association study data

Lipid droplets fuel SARS-CoV-2 replication and production of inflammatory mediators

Phenome-wide association studies across large population cohorts support drug target validation

An interactive web-based dashboard to track COVID-19 in real time

Chemokine Receptor Redundancy and Specificity Are Context Dependent

Genomewide Association Study of Severe Covid-19 with Respiratory Failure

Upstream Transcription Factor 1 (USF1) allelic variants regulate lipoprotein metabolism in women and USF1 expression in atherosclerotic plaque

MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data

Genome-wide association study identifies multiple susceptibility loci for pulmonary fibrosis

DPP8/9 inhibitors are universal activators of functional NLRP1 alleles

Pulmonary fibrosis and COVID-19: the potential role for antifibrotic therapy

Complex Immune Dysregulation in COVID-19 Patients with Severe Respiratory Failure

Clinical Characteristics of Coronavirus Disease 2019 in China

The challenges, advantages and future of phenome-wide association studies

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Fine-mapping genetic associations

Coagulopathy in COVID-19

The COVID-19 Host Genetics Initiative, a global initiative to elucidate the role of host genetic factors in susceptibility and severity of the SARS-CoV-2 virus pandemic

Genetic analysis of quantitative traits in the Japanese population links cell types to complex human diseases

eQTL Catalogue: a compendium of uniformly processed human gene expression and splicing QTLs

NLRP1-Dependent Pyroptosis Leads to Acute Lung Injury and Morbidity in Mice

USF1 deficiency activates brown adipose tissue and improves cardiometabolic health

Hematological features of persons with COVID-19

Haematological characteristics and risk factors in the classification and prognosis evaluation of COVID-19: a retrospective cohort study

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

Longitudinal analyses reveal immunological misfiring in severe COVID-19

LDlink: a web-based application for exploring populationspecific haplotype structure and linking correlated alleles of possible functional variants

Molecular mechanisms of inflammasome signaling

Roles of CD48 in regulating immunity and tolerance

COVID-19: consider cytokine storm syndromes and immunosuppression

Pathological inflammation in patients with COVID-19: a key role for monocytes and macrophages

Cytokine release syndrome in severe COVID-19

DPP8 and DPP9 inhibition induces pro-caspase-1-dependent monocyte and macrophage pyroptosis

Simultaneous TCR and CD244 Signals Induce Dynamic Downmodulation of CD244 on Human Antiviral T Cells

Estimation of the multiple testing burden for genomewide association studies of nearly all common variants

Meta-analysis of genome-wide association studies for body fat distribution in 694 649 individuals of European ancestry

COVID-19 interstitial pneumonia: monitoring the clinical course in survivors

Genome-wide association study of body fat distribution identifies adiposity loci and sex-specific genetic effects

Tetraspanins CD9 and CD151 at the immune synapse support T-cell integrin signaling

The Kynurenine Pathway is a Double-Edged Sword in Immune-Privileged Sites and in Cancer: Implications for Immunotherapy

USF1 deficiency alleviates inflammation, enhances cholesterol efflux and prevents cholesterol accumulation in macrophages

COVID-19 severity associates with pulmonary redistribution of CD1c+ DC and inflammatory transitional and nonclassical monocytes

From genome-wide associations to candidate causal variants by statistical fine-mapping

echolocatoR: an automated end-to-end statistical and functional genomic fine-mapping pipeline

Leveraging Genome and Phenome-Wide Association Studies to Investigate Genetic Risk of Acute Lymphoblastic Leukemia

Comprehensive Integration of Single-Cell Data

Familial combined hyperlipidemia: An overview of the underlying molecular mechanisms and therapeutic strategies

Anticoagulant treatment is associated with decreased mortality in severe coronavirus disease 2019 patients with coagulopathy

PNPLA2 influences secretion of triglyceride-rich lipoproteins by human hepatoma cells

COVID-19 infection alters kynurenine and fatty acid metabolism, correlating with IL-6 levels and renal status

Autoproteolysis of PIDD marks the bifurcation between pro-death caspase-2 and prosurvival NF-κB pathway

Immunology of COVID-19: current state of the science

DPP8/DPP9 inhibition elicits canonical Nlrp1b inflammasome hallmarks in murine macrophages

Coronavirus biology and replication: implications for SARS-CoV-2

Bayes factors for genome-wide association studies: comparison with P-values

The role of high cholesterol in age-related COVID19 lethality

A global overview of pleiotropy and genetic architecture in complex traits

Immune cell profiling of COVID-19 patients in the recovery stage by single-cell sequencing

Factors associated with COVID-19-related death using OpenSAFELY

Kynurenic Acid: The Janus-Faced Role of an Immunomodulatory Tryptophan Metabolite and Its Link to Pathological Conditions. Front Immunol 8, 1957

Inflammasomes and Pyroptosis as Therapeutic Targets for COVID-19

Comorbidity of Alcohol Use Disorder and Chronic Pain: Genetic Influences on Brain Reward and Stress Systems

Immunosuppressive IDO in Cancer: Mechanisms of Action

An integrative functional genomics framework for effective identification of novel regulatory variants in genome-phenome studies

Human DPP9 represses NLRP1 inflammasome and protects against autoinflammatory diseases via both peptidase activity and FIIND domain binding

Obesity & genetic predisposition with COVID-19