ACE2 Netlas: In silico Functional Characterization and Drug-Gene Interactions of ACE2 Gene Network to Understand Its Potential Involvement in COVID-19 Susceptibility

Angiotensin-converting enzyme-2 (ACE2) receptor has been identified as the key adhesion molecule for the transmission of the SARS-CoV-2. However, there is no evidence that human genetic variation in ACE2 is singularly responsible for COVID-19 susceptibility. Therefore, we performed an integrative multi-level characterization of genes that interact with ACE2 (ACE2-gene network) for their statistically enriched biological properties in the context of COVID-19. The phenome-wide association of 51 genes including ACE2 with 4,756 traits categorized into 26 phenotype categories, showed enrichment of immunological, respiratory, environmental, skeletal, dermatological, and metabolic domains (p < 4e-4). Transcriptomic regulation of ACE2-gene network was enriched for tissue-specificity in kidney, small intestine, and colon (p < 4.7e-4). Leveraging the drug-gene interaction database we identified 47 drugs, including dexamethasone and spironolactone, among others. Considering genetic variants within ± 10 kb of ACE2-network genes we identified miRNAs whose binding sites may be altered as a consequence of genetic variation. The identified miRNAs revealed statistical over-representation of inflammation, aging, diabetes, and heart conditions. The genetic variant associations in RORA, SLC12A6, and SLC6A19 genes were observed in genome-wide association study (GWAS) of COVID-19 susceptibility. We also report the GWAS-identified variant in 3p21.31 locus, serves as trans-QTL for RORA and RORC genes. Overall, functional characterization of ACE2-gene network highlights several potential mechanisms in COVID-19 susceptibility. The data can also be accessed at https://gpwhiz.github.io/ACE2Netlas/.

INTRODUCTION SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) is the causative agent responsible for recent global spread of COVID-19 (coronavirus disease 2019) (Wu et al., 2020;Zhou et al., 2020). Millions of people have been infected with the virus, which caused global lockdowns and heavily restricted interpersonal contact. These measures were taken to reduce viral spread through respiratory droplet exchange between persons.
SARS-CoV-2 is capable of entering the host cells via ACE2 (angiotensin converting enzyme 2) (Walls et al., 2020;Mercurio et al., 2021). ACE2 is found on many different cell types, which normally helps regulate blood pressure and inflammation through cleavage of angiotensin II (ANG II) (Hamming et al., 2007). The virus occupies cell-surface of ACE2 leading to accumulation of angiotensin (ANGII), inflammation, and cell death (Walls et al., 2020). The interactions between the spike protein and ACE2 trigger pre-/post-fusion conformational changes at the spike protein, following spike cleavage in two main domains. The spike N-terminal domain forms a stable protein-protein complex with ACE2, whereas the C-terminal domain, namely the spike post-fusion protein, favors virus/hostcells membrane fusion (Mercurio et al., 2021). In the lungs, SARS-CoV-2 mediated ANGII accumulation leads to alveolar cell death and a reduction in oxygen uptake .
Although ACE2 is the cellular entry point, there is little evidence that genetic variation in ACE2 is singularly responsible for COVID-19 susceptibility (Ellinghaus et al., 2020;Pairo-Castineira et al., 2020;Shelton et al., 2020). Due to the functional role of ACE2 in SARS-CoV-2 infection, we hypothesize that genes interacting with ACE2 activity are enriched for molecular pathways relevant for COVID-19 susceptibility. Accordingly, we employed a top-down approach to analyze tissue-specific transcriptomic regulation, drug-gene interactions, and variant prioritization using genetic variants within the ACE2 gene-gene connectome and protein-protein interaction networks. With this approach we identified several biological processes and functional effects of ACE2-gene network relevant for the vast symptoms observed following SARS-CoV-2 infection.

RESULTS
A study overview is presented in Supplementary Figure 1.

The ACE2 Gene Connectome
A total of 60 genes were identified from six network databases that interact with ACE2 (Supplementary Table 1).

Gene-Drug Interaction and Over-Represented Biological Functions
To identify known drugs that interact with the ACE2-gene set, we investigated the drug-gene interaction database (DGIdb) (Griffith et al., 2013). Out of 61 genes, 29 had information about their drug-gene interaction in DGIdb, resulting in 238 unique drug-gene observations (Supplementary Table 4). Some of the notable drugs observed via this approach were spironolactone, dexamethasone, metformin, and hydrocortisone. To understand the role of these drugs in affecting biological processes, we performed drug-set enrichment analysis. DSEA (Napolitano et al., 2016) found gene-ontology mapping for 47 drugs and tested against REACTOME gene ontology database. Although the results did not survive Bonferroni correction, nominally significant enrichments were observed for platelet sensitization by low-density lipoprotein cholesterol (p = 0.003), IL-7 signaling (p = 0.004), glycerophospholipid biosynthesis (p = 0.005), and viral messenger RNA synthesis (p = 0.011) (Figure 2 and Supplementary Table 5).

Neanderthal LA Introgression Within ACE2 Network SNPs
Due to the Neanderthal introgression observed in 3p21 locus as risk to COVID-19 (Zeberg and Pääbo, 2020), we compared mean probability of Neanderthal LA between the ACE2-network SNP set (mean = 0.032) and 1,000 randomly
We further created a unique comprehensive network by prioritizing genes from the ACE2-network using transcriptomic profile specific to SARS-CoV-2, over-represented domains for traits associated with ACE2-network genes, and the gene-drug targets that were overrepresented for biological functions (Figure 5).

DISCUSSION
ACE2 is expressed in several tissues and plays a key role in host-entry of SARS-CoV-2 (Hoffmann et al., 2020). However, the genomic profile of ACE2 is limited in explaining the vast symptomology observed for COVID-19. Understanding ACE2 associated molecular networks presents several functional insights between genetic targets based on gene expression, topology, and protein and signaling relationships (Huang et al., 2018). Due to the well-characterized role of ACE2 in SARS-CoV-2 infection, we generated novel information regarding the molecular and phenotypic characteristics of ACE gene network in the context of their potential involvement in COVID-19 susceptibility. Our PheWAS-based analysis showed Frontiers in Genetics | www.frontiersin.org 6 August 2021 | Volume 12 | Article 698033 FIGURE 5 | Network of ACE2 interacting genes. The overview of ACE2-network genes (blue) as drug targets (green), phenotype domains (red), and gene expression category (purple) wherein the genes are significant from the CZBioHub (see text for details). The pink edges represent upregulation and blue edges represent downregulation, and yellow edges connect genes to overrepresented phenotype domains.
that genetic variation within ACE2 gene network is associated with immunity (white blood cell, neutrophil count, lymphocyte count), respiratory (FVC, asthma, DVT), and metabolic traits (BMI, cholesterol, body measurements). This is in line with known epidemiology of COVID-19 and its comorbidities (Ejaz et al., 2020;Gardinassi et al., 2020). The expression of ACE2-network genes was enriched for regulatory mechanisms related to small intestine, colon, kidney, and liver. It is hypothesized that furin, a serine protease present in lungs but also highly expressed in small intestine, and is involved in the cleavage of S-spike for attachment of the ACE2 receptor (Mönkemüller et al., 2020). Patients with kidney disease have higher risk for COVID-19 severe symptoms (Ajaimy and Melamed, 2020). Additionally, the inflammation and cytokine storm from COVID-19 is observed to damage kidney tissues . Lastly, modest increase in liver enzymes has been associated with COVID-19, and returning to baseline during the recovery phase (Pawlotsky, 2020).
Understanding the genes that interact with ACE2 receptor has potential to understand drug-targets and molecular processes that might play a role in susceptibility and treatment response of COVID-19. The drug-gene interaction analysis within ACE2 network identified dexamethasone, reported to lower mortality in COVID-19 cases requiring mechanical ventilation (RECOVERY Collaborative Group et al., 2020). Drugs-spironolactone and hydrocortisone target the androgen system. Androgen signaling modulates ACE2 expression and elevated androgen levels have been associated with severe symptomology of COVID-19 (Samuel et al., 2020). Spironolactone is a diuretic and alleviates respiratory symptoms by reducing fluid from the lungs (Cadegiani et al., 2020). The use of spironolactone is currently being tested for acute respiratory distress syndrome in COVID-19 patients (Dumanlı et al., 2020). Hydrocortisone is currently under clinical trials for treating COVID-19 related hypoxia symptoms (Petersen et al., 2020). Among the other compounds identified, metformin, a known drug for treating diabetes, can also affect respiratory outcomes (Yen et al., 2020). A recent study reported protective effects of metformin in women with diabetes and obesity who were admitted with COVID-19 diagnosis (Bramante et al., 2020). Lastly, melatonin has been hypothesized to improve general immunity and lower oxidative stress generated from SARS-CoV-2 infection (Shneider et al., 2020).
The miRNA target sites altered by ACE2-network SNPs identified miR-302b and miR-181d as over-represented miRNA clusters. The downregulated expression of miR-302b has been observed to reduce survival rates in chronic obstructive pulmonary disease (COPD) patients (Keller et al., 2019). A metaanalysis showed that COPD diagnosis increased susceptibility to COVID-19 (Lippi and Henry, 2020). The miRNA-181 cluster has been associated with regulation of TNF-alpha (Zhu et al., 2017), T-cell aging (Ye et al., 2018) and emphysema (Osei et al., 2015). miRNA-17 and 106 belong to same miRNA family, miRNA-17 is upregulated in bronchoalveolar stem cells to lower SARS-CoV replication (Mallick et al., 2009). An in silico study of miRNA targets for SARS-CoV-2 genomic sequence found miRNA-17 as one of the targets with experimental evidence of its upregulation in H7N9 Influenza virus infection (Khan M.A.A.K. et al., 2020). The top over-represented diseases in miRNA-ACE2-network-SNPs were diabetes, hepatitis C viral infection, heart failure and Alzheimer's disease. A greater number of diabetic individuals with COVID-19 have been reported to require hospitalization than non-diabetic individuals (Apicella et al., 2020). Furthermore, SARS-CoV-2 infection contributes in the development of ketosis in diabetic individuals resulting in longer length of hospitalization stay . Triglyceride and glucose index was associated with severity of COVID-19 (Ren et al., 2020). While there are limited studies about hepatitis C in COVID-19 patients (Richardson et al., 2020), heart failure was reported by multiple studies as being associated with COVID-19 severity (Hanley et al., 2020;Yancy and Fonarow, 2020). Alzheimer's disease is another condition associated with COVID-19 susceptibility (Chang et al., 2020), including APOE4 carrier status with increased risk of severe COVID-19 (Kuo et al., 2020).
In contrast to specific enrichment of Neanderthal LA in a COVID-19 risk locus on chromosome 3 (Zeberg and Pääbo, 2020), there is no evidence of increased Neanderthal LA in the ACE2 network investigated here. This suggests that, although some loci conferring risk for COVID-19 severity, such as the one identified on chromosome 3, may have originated from Neanderthal admixture events, this mechanism did not shape the genetic architecture of the ACE2 network responsible for entry of SARS-CoV-2 into host cellular machinery.
Lastly, among ACE2-network-SNPs, potential COVID-19 risk alleles were observed in RORA gene with respect to multiple COVID-19 phenotypes. RORA protein product is involved in immune response, cancer and metabolism (Cook et al., 2015). RORA plays a role in the activation of T helper cells during lung inflammation by regulating tumor necrosis factor and interleukins (Nejati Moharrami et al., 2018;Haim-Vilmovsky et al., 2019), and was upregulated in cardiomyocytes infected with SARS-CoV-2 (Hachim et al., 2020). The hypothesis-free approach of genome-wide association of hospitalized COVID-19 vs. the population highlighted SLC6A20 with genome-wide significance on chromosome 3 locus. The SLC12 (SLC12A6) class is responsible for transport of inorganic ions such as sodium and chloride while the SLC6 class (SLC6A19, identified via network approach and SLC6A20, identified via genome-wide approach) are responsible for transport of amino acids such as glutamate and glycine which are important for neurotransmitter activity (Lin et al., 2015). SLC6A19 (among other SLC-class genes) serves similar function to SLC6A20, both are expressed in the intestinal tissue and contingent upon ACE2 expression (Vuille-Dit-Bille et al., 2020). Multiple studies report more than 10% of the COVID-19 confirmed patients exhibit gastrointestinal symptoms Lian et al., 2020).
Although we provided a wide range of information highlighting the molecular and phenotypic characteristics of ACE2 gene network and their putative implications with COVID-19 risk, the findings reported have to be considered exploratory. We used appropriate computational methods and statistical approaches to generate reliable evidence useful to open new directions in COVID-19 research. We also highlighted when the results reported did not survive stringent multiple testing correction. This limitation is particularly relevant with respect to the ACE2 network genetic associations. Due to the limited statistical power of the genome-wide data available for the Freeze 3 data from the COVID-19 Host Genetics Initiative, none of the risk alleles identified as functionally relevant survive genome-wide testing correction. Future work from the HGI will potentially lead to more risk loci being identified. Further analyses will be needed to validate our current findings.

CONCLUSION
ACE2 is one of the few molecular targets recognized to play a key role in the COVID-19 pathogenesis. We conducted a comprehensive analysis leveraging multiple resources (e.g., drug-gene interactions, tissue-specific transcriptomic profile, and phenome-wide and genome-wide datasets) to expand our understanding of the genomic characteristics of the host ACE2 gene network. Overall, our findings incorporate multi-tiered epigenomic, transcriptomic, and genomics of the known ACE2network which highlight the potential mechanisms linking ACE2 systems biology to COVID-19 susceptibility and its possible comorbidities.

Gene Network Collection
Information regarding ACE2 gene network was mined from GeneMANIA (Franz et al., 2018), Stringdb (Szklarczyk et al., 2017), Agile Protein Interactomes Database (APID) (Prieto and De Las Rivas, 2006), GeneNetwork (Deelen et al., 2019), Biogrid (Oughtred et al., 2019), and FunctionalNet (Lee et al., 2011) for Homo sapiens organism, last searched on June 27, 2020. Since all resources use different algorithms, using ACE2 as query gene, immediate genes connections that were available in each databank were retrieved with their default settings. Removing overlapping genes across the six databases, resulted in 61 unique genes (60 genes plus ACE2) (Supplementary Figure 1 and Supplementary Table 1). Specifically, GeneMANIA uses automatically selected weighting method for 20 max resultant genes with 10 max resultant attributes. In Stringdb, the ACE2 gene query with medium confidence interaction score, and including all active interaction source except text mining. For APID, single query of ACE2 generated a network of 11 interacting genes. In GeneNetwork, we used ACE2 as query genes and selected co-regulated genes with evidence of p < 2e-12. In Biogrid, we used ACE2 in Homo Sapiens as query, and used the network information. In FunctionalNet, the ACE2 was searched using its ENTREZ id (59272) and selected the genelist that interact with ACE2. We focused on the genes that interact immediately with ACE2 because including more genes that interact through intermediate genes will result in large volume data and difficult to interpret. Several genes were identified by multiple sources listed in Supplementary Table 1, and a total of 61 genes including ACE2 were investigated for their characteristics. The genomic coordinates for the genes were annotated using biomart (Durinck et al., 2009), ensemble GRCh37/hg19. The analysis and visualization were performed in R 3.6.

Tissue-Specific Transcriptomic Regulation
The tissue specificity was tested for 60 ACE2-interacting genes in FUMA (Watanabe et al., 2017). The input genes were tested for pre-calculated tissue-specific differentially expressed genes from the GTEx v8 (Aguet et al., 2019). We also considered the t-statistic sign for up and down-regulated genes against protein coding genes as background. Enrichments were performed using hypergeometric tests and significant enrichments were defined according to Bonferroni corrected p-value < 0.05.

Gene Expression of ACE2-Interacting Genes in Upper Respiratory Tissue for SARS-CoV-2 and Other Viruses
To understand which genes from the ACE2-interacting genes are differentially expressed in SARS-CoV-2 and other viruses, we extracted these genes from transcriptomic study of acute respiratory illnesses for COVID-19 patients (N = 93), other viral (N = 41) or non-viral (N = 100) in the upper respiratory tract tissue (Mick et al., 2020). Their data was extracted from https://github.com/czbiohub/covid19-transcriptomicspathogenesis-diagnostics-results. We used ENSEMBL identifiers of 61 ACE2 interacting genes and were able to extract 35/61 genes. The study performed three gene expression comparisons, SARS-CoV-2 vs. no-virus, SARS-CoV-2 vs. other-virus, othervirus vs. no-virus, and genes with Benjamin-Hochberg adjusted p-value < 0.05 were assigned as significant to the respective comparison.

Phenome-Wide Analysis of ACE2 Gene Network
A phenome-wide association study (PheWAS) was performed for 51 of 61 genes that were present in GWASAtlas (Watanabe et al., 2019) using all traits available per gene. Statistical significance was determined by applying a Bonferroni multipletesting correction accounting for the number of GWAS traits (4,765 traits) available in the GWASAtlas (p < 1.05 × 10 −5 ). Each trait was grouped into a domain (Supplementary Table 5) which was tested for enrichment using one-sided Fisher's exact test for high proportion of significant traits vs. all others tested. A significant domain enrichment was defined considering a Bonferroni-corrected threshold accounting for the number of domains tested (p-value < 0.0019; 0.05/26).

Gene-Drug Interactions and Biological Functions
Information on drugs that interact with ACE2 network genes were extracted from The Drug-Gene Interaction database (DGIdb) (Griffith et al., 2013) followed by drug-set enrichment for over represented biological functions using DSEA (Drug-Set Enrichment Analysis) (Napolitano et al., 2016).

Characterization of SNPs
Single nucleotide polymorphism (SNPs) were extracted based on the genomic coordinates of the genes (± 10 kb) for GrCh37; dbSNP153 from the UCSC browser (Haeussler et al., 2019) using bigbed utilities (Karolchik et al., 2004), and referred to as "ACE2-network SNPs." ACE2-network SNPs were annotated for global allele frequency, Combined Annotation-Dependent Depletion (CADD) score (Rentzsch et al., 2019), deep learning based algorithm framework (DeepSEA) (Zhou and Troyanskaya, 2015), and target miRNAs using SNPnexus (Dayem Ullah et al., 2018). DeepSEA is a deep learning-based algorithmic framework for predicting the chromatin effects of sequence alterations with single nucleotide sensitivity (Zhou and Troyanskaya, 2015). The identified miRNAs were tested for over-represented miRNA clusters, functions, and diseases using TAM 2.0 .

Neanderthal Introgression
Motivated by evidence of a chromosome 3 COVID-19 risk locus enriched of Neanderthal local ancestry (LA) (Zeberg and Pääbo, 2020), we compared the distribution of probability of Neanderthal LA in our COVID-19 ACE2-network SNP set and 1,000 randomly sampled SNP sets comprised on SNPs across the genome with comparable genomic features. ACE2-network SNPs were mapped using previously defined Neanderthal LA data (Sankararaman et al., 2014;Durvasula and Sankararaman, 2019). A total of 6,822 LD-independent pairwise SNPs (r 2 = 0.1 and p = 0.1 in 250 kb window size) were used as standard input for SNPsnap (Pers et al., 2015). In SNPsnap, 1,249/6,822 independent ACE2 network SNPs could be matched based on the following genomic features relative to the input SNP list: minor allele frequency within 2%, gene density within 50%, nearest gene within 50%, and number of linkage disequilibrium groups within 50%. SNPsnap was instructed to exclude the ACE2network SNP list from the pool of eligible feature-matched SNPs. Non-parametric Wilcoxon rank sum tests were used to compare the Neanderthal LA of our ACE2 network SNP list to that of all 1,000 random SNP sets and multiple testing correction was applied to adjust for a false discovery rate of 5%. The SNPs of the ACE2 network were extracted and pruned for LD and p-value using plink 1.9. The multiple testing correction was applied using Bonferroni p-value < 0.05. These significant SNPs were annotated further for pathogenicity using Combined Annotation Dependent Depletion (CADD) score and their role as quantitative trait loci (QTL) for gene expression using GTEx, and methylation using QTLbase . The trans-eQTL relationship of GWAS-reported locus-3p21.31 were identified from eQTLgen (Võsa et al., 2018).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher. The data presented is available in Supplementary Files and also on https://gpwhiz.github.io/ ACE2Netlas/.

ETHICS STATEMENT
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. The patients/participants provided their written informed consent to participate in their original studies reported.

AUTHOR CONTRIBUTIONS
GP conceptualized the study design, analyzed, and drafted the manuscript. FW contributed to analysis and manuscript writing. AG, DK, FD, and RP contributed to result interpretation, manuscript drafting and revision. RP supervised the study and finalized the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
We would like to acknowledge support from the National Institutes of Health (R21 DC018098, R21 DA047527 and F32 MH122058).