{'fraction_non_alphanumeric': 0.04033765792306385, 'fraction_numerical': 0.03671180103110305, 'mean_word_length': 5.07311199036642, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 2, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phenome-wide and eQTL Associations of COVID-19 Genetic Risk Loci


 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.



INTRODUCTION
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

PheWAS of genetic risk loci associated with severe COVID-19
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.

J o u r n a l P r e -p r o o f
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.

Immune-specific eQTL analysis
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.

GWAS/eQTL Colocalization Analysis
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.

Differential gene expression in bronchoalveolar lavage fluid using single-cell analyses
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.

DISCUSSION
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.,

2021) can be sensed by NLRP1 and set off an inflammatory cascade mediated by inflammasome
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.

LIMITATIONS OF THE STUDY
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.

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

Lead Contact
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).

Materials Availability
This study did not generate new unique reagents.

Data and Code Availability
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

METHOD DETAILS
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.

J o u r n a l P r e -p r o o f
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.

GWAS/eQTL colocalization analyses.
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.

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