SARS-CoV-2 meta-interactome suggests disease-specific, autoimmune pathophysiologies and therapeutic targets

Virus- Host interaction protein networks

We considered the interactomes from 7 human infecting members of the Coronaviridae family and from other viruses causing respiratory diseases (Influenza A H1N1, Influenza A H7N9, Influenza B, HRSV) (Extended data¹³).

The SARS-CoV-2 host-protein network was recently uncovered by an affinity-purification mass spectrometry experiment⁹. Other host-pathogen protein interaction maps were derived from p-HIPSTER, a structure-informed atlas of human-virus interactions¹⁴, which contains both interactions present in the RCSB Protein Data Bank and predicted pathogen-host PPIs.

SARS-CoV-2 network expansion

We performed the expansion using the latest upgrade of Human Integrated Protein-Protein Interaction rEference (HIPPIE)¹⁵ which integrates confidence scored and functionally annotated human protein-protein interactions (PPIs) from the public databases BioGRID, DIP, MIPS, HPRD, IntAct, MINT and BIND.

To guarantee the interaction confidence, we accounted only for physical interactions, starting from a very high PPI confidence score (0.95), and progressively expanded the network coverage until reaching a medium confidence score of 0.65 (second quartile of the HIPPIE score distribution).

Disease-associated genes and statistical enrichement analysis

We extracted, from the NHGRI-EBI Catalog of human genome-wide association studies¹⁶, the mapped genes of SNPs significantly associated (p-value ì≤ 5 × 10⁻⁸) with immune-mediated disorders [multiple sclerosis (MS), systemic lupus erythematosus (SLE), rheumatoid arthritis (RA), Crohn’s disease (CD), ulcerative colitis (UC), type 1 diabetes mellitus (T1D)] and with conditions being relevant to COVID-19 clinical manifestations and prognosis [obesity, aging, type 2 diabetes mellitus (T2D), asthma, chronic obstructive pulmonary disease (COPD), hypertension and coronary artery disease (CAD)]. Specifically, we searched for the corresponding trait in the database and collected the reported loci (Extended data¹³).

We overlapped the interactomes and the disease-associated genes and evaluated the statistical significance using geneOverlap package version 1.24.0¹⁷ in R, version 3.6.3. We took into consideration the Jaccard index, which represents the similarity between two sets; p-value significance cutoff was set at p<0.05 after correction with Benjamini-Hochberg method.

COVID-19 meta-interactome construction and pathway analysis

We selected the genes included in the statistically significant intersection of SARS-CoV-2 network and disease associated gene sets, constructing a COVID-19 meta-interactome. To better understand the immune mechanism shared by COVID-19 and immunological diseases, we also defined a COVID-19 immunological subnetwork collecting genes from immune-mediated disease only.

We then performed a functional enrichment analysis of each significant overlapping subset of genes among GO Biological Processes and Reactome Gene sets using Metascape¹⁸, a web-based utility. This tool uses the MCODE (Molecular Complex Detection) algorithm to perform module analysis and detect dense regions of protein interaction networks¹⁹. Network visualization and node connection analysis was performed with Cytoscape 3.8.0²⁰.

The Circle plot representing COVID-19 meta-interactome genes was produced using the Circos table viewer tool²¹. The UpSet plot was produced using the UpSetR package for R²², version 1.4.0.

Gene expression analysis

We analyzed the COVID-19 meta-interactome genes using the GENE2FUNC of the FUMA tool²³, an integrative web-platform to perform post-GWAS annotation. We specifically performed a gene expression analysis based on GTEx v8 and Ensemble version v92, retrieving tissue specificity. Gene-enrichment statistical analysis was set as follows: the adjusted P value cut-off was 0.05, minimum overlapping genes with gene-sets was 2, method for multiple test correction was Benjamini-Hochberg (FDR).

We looked for the expression of COVID-19 meta-interactome genes in samples derived from lungs of infected humans. To do this, we selected from the NCBI Gene Expression Omnibus Database²⁴ the series GSE147507²⁵, including transcriptional profiling of post-mortem lung samples of two male COVID-19 patients compared with healthy lung biopsies, and transcriptional profiling of normal human bronchial epithelial (NHBE) cells after the treatment with human IFN-I-Β.

Briefly, at first we extracted the complete results using the GEO RNA-seq Experiments Interactive Navigator (GREIN) web tool²⁶; we then selected the COVID-19 meta-interactome genes and plot the final data with EnhancedVolcano R package, version 1.4.0.

To assess the potential of COVID-19 meta-interactome genes of being regulated by Type 1 interferon-β, we exploited Interferome 2.0.1, an online database of interferon-regulated genes (IRGs)²⁷. The search settings were as follows: Interferon type I, subtype IFN-β, any treatment concentration and time, in vitro and in vivo cell lines derived from Homo sapiens, any sample type, Fold change cutoffs 1.0.

COVID-19 meta-interactome drug target screening

To assess the therapeutic potential of COVID-19 meta interactome, we used the Scalable Precision Medicine Oriented Knowledge Engine (SPOKE)²⁸. SPOKE is a heterogeneous knowledge network that includes data from 29 publicly available biomedical databases. We searched the SPOKE Neighborhood Explorer for the COVID-19 meta-interactome genes with the highest Betweenness Centrality score in the whole network and in subnetworks, filtering results with the following settings: “compound” and “gene” nodes; “compound-binds-protein”, “compound-downregulates-gene”, “compound-upregulates-gene”. We filtered nodes representing compounds at development phase ≥3; Edges Attributes and Limits were as set from the database.

SARS-CoV-2 interactors enrichment

In this work, we downloaded the host-protein networks of SARS-CoV-2, six other HCovS (SARS-CoV, MERS-CoV, HCoV-229E, HCoV-NL63, HCoV-OC43, HCoV-HKU1) and four other viruses causing respiratory diseases in humans (Influenza A H1N1, Influenza A H7N9, Inluenza B, Human Respiratory Syncitial Virus). We extended the experimentally determined SARS-CoV-2-host interactions⁹ by including second neighbors, i.e. molecules directly interacting with the host proteins engaged by the virus extracted from HIPPIE. We progressively lowered the confidence score threshold to construct 7 more sets of increasing size (Extended data¹³).

We then looked for an enrichment of human proteins interacting with SARS-CoV-2 within sets of genes mapped in genome-wide studies of immune-mediated conditions and comorbidities impacting on COVID-19 prognosis.

The analysis revealed a significant enrichement of SARS-CoV-2 interactome with genes associated with several autoimmune diseases (Figure 1). The strongest significance was reached by SARS-CoV-2/MS overlap (p=0.0002), consolidated at multiple expansion levels; significant associations with CD, SLE and RA were of weaker magnitude.

When looking at clinically relevant comorbidities, there was a significant overlap between SARS-CoV-2 interactors and T2D associated genes, that followed a positive trend in parallel to the network expansion, with a maximum of 65 shared interactors (Figure 1 and Extended data¹³). The other COVID-19 relevant comorbidities reaching statistical significance were CAD and Asthma.

Figure 1. Viral interactomes enrichment in GWAS loci.

The heatmap shows the Jaccard Index of the overlaps of gene sets associated with immunological diseases and COVID-19 relevant comorbidities with viral interactomes. SARS-CoV-2 section shows in the first line the host protein network derived from experimental evidence⁹, followed by results obtained overlapping the expanded networks generated with HIPPIE settings at different confidence levels (0,95 to 0,65). HRSV: Human Respiratory Syncitial Virus; MS: Multiple Sclerosis; RA: Rheumatoid Arthritis; SLE: Systemic Lupus Erythematosus; CROHN: Crohn’s Disease. U.COLITIS: Ulcerative Colitis; CELIAC D.: Celiac Disease; T1D: Type 1 Diabetes Mellitus; T2D: Type 2 Diabetes Mellitus; CKD: Chronic Kidney Disease; COPD: Chronic Obstructive Pulmonary Disease; HTN: Hypertension; CAD: Coronary Artery Disease.

Neither the other HCovs nor the other respiratory viruses (H1N1, H7N9, Influenza B, HRSV) reached significance levels as high as SARS-CoV-2.

In summary, SARS-CoV-2 host protein network is enriched in sets of genes associated with multiple autoimmune diseases and shares significant interactions with three conditions influencing COVID-19 prognosis: asthma, T2D and CAD.

COVID-19 meta-interactome construction and analysis

We next constructed a COVID-19-human meta-interactome, selecting genes deriving from the significant overlaps with immunological diseases and comorbidities. The resultant network was composed of 248 genes (Figure 2A and Extended data¹³), with a non-specific distribution on autosomes, and a single gene located on Chromosome X (MECP2) and none in chromosome Y (Figure 2B).

To assess the prominence of interactors within the network, we ranked the gene according to the Betweenness score, a network centrality metric which reflects the number of shortest paths that go through a given node, thus representing the importance of that node on the information flow through the network (Figure 2C). Top ranking nodes were: MYC, a widely acting oncogene with a key role in immune regulation; XPO1, whose encoded protein mediates the nuclear export of cellular proteins and of RNAs, including viral ones; SMAD3, an intracellular transducer of transforming growth factor-β signaling pathway (Figure 2D). We then performed a functional enrichment analysis of each significant overlapping subset of genes using Metascape¹⁸ a web-based utility which uses the MCODE (Molecular Complex Detection) algorithm, to perform module analysis and detect dense regions of protein interaction networks¹⁹. Densely interconnected genes subsets pertain to central immunological processes (NFĸB, immune response and TLR3 pathway, podosome assembly, MAPK cascade), cell cycle and transcriptional regulation (TCF/LEF and WNT pathway), cholesterol and metabolism, and ionic channels regulation (Figure 3). In line with this data, the COVID-19 meta-interactome pathway enrichment analysis (Figure 4) showed that the top 20 clusters belonged to immunological (lymphocyte differentiation, neutrophils pathway, immune response-regulating signaling pathway, regulation of cytokine production), and cell cycle processes; interestingly, we found an enrichment of the hemostasis process. This is in line with evidence showing a coagulation dysfunction in severe COVID-19²⁹

Figure 2. COVID-19 meta-interactome.

A) Network Representation: node size is proportional to Betweenness Centrality score. B) Circle plot showing the chromosome distribution of the 248 COVID-19 meta-interactome genes. C) Top 10 nodes ranked by Betweenness Centrality D) Detail of the network view showing the top nodes.

Figure 3. COVID-19 meta-interactome clusters.

The M-CODE algorithm was applied to detect densely interconnected node subsets. The network representation displays a different color for each cluster, associated with the enriched biologic processes in the relative boxes.

Figure 4. COVID-19 meta-interactome pathway analysis.

Heatmap of top 20 enriched terms across input gene lists, colored by p-values. MS: Multiple Sclerosis; RA: Rheumatoid Arthritis; SLE: Systemic Lupus Erythematosus; CROHN: Crohn’s Disease. U. COLITIS: Ulcerative Colitis; CELIAC D.: Celiac Disease; T1D: Type 1 Diabetes Mellitus; T2D: Type 2 Diabetes Mellitus; CAD: Coronary Artery Disease.

Notably, in COVID-19 meta-interactome composition, most of the genes deriving from autoimmune diseases’ sets appeared to be shared among the traits: the most relevant example is seen for MS and RA, whose SARS-CoV-2 interacting proteins derived from the same genes (specifically, 49 of the 61 GWAS-associated identical genes). These data suggested the presence of a “core” immune mechanism driving COVID-19 physiopathology shared with most of autoimmune conditions. Conversely, COVID-19 meta-interactome elements belonging to comorbidities’ subsets resulted more specific (T2D 52/66, CAD 12/20, asthma 29/61 unshared genes) (Figure 5A,B,C). This suggested the possibility that the pathobiological mechanisms underlying the worse prognosis in COVID-19 patients affected by these comorbidities are at least partly independent from the imbalanced immune response to SARS-CoV-2 infection.

Figure 5. COVID-19 meta-interactome composition.

A) The table shows the subsets of genes from each trait contributing to the network composition. B) The circle plot shows the overlap between gene lists at the gene level, with purple curves linking identical genes. C) The UpsetPlot displays a schematic representation of the intersections among gene sets (x axis), with the respective size (y axis). MS: Multiple Sclerosis; RA: Rheumatoid Arthritis; SLE: Systemic Lupus Erythematosus; CROHN: Crohn’s Disease. U. COLITIS: Ulcerative Colitis; CELIAC D: Celiac Disease; T1D: Type 1 Diabetes Mellitus; T2D: Type 2 Diabetes Mellitus; CAD: Coronary Artery Disease.

To assess this hypothesis, we next performed a separate analysis of the immune subnetwork, containing genes derived from immunological diseases’ subsets (n=150), and the networks of each comorbidity (Asthma, n=61; T2D, n=66, CA D, n=20). The immune subnetwork pathway analysis confirmed a global enrichment of immunological functions. Betweenness Centrality ranking highlighted again MYC as the most relevant node, followed by RHOA, encoding for a small GTPase involved in cytoskeleton assembly and nuclear mechanotransduction, and TRAF6, a multiacting member of the toll-like receptor (TLR) family, capable of activating MAPK, PI3K, and interferon regulatory factor (IRF) pathways. These nodes also showed a high Betweenness Centrality score in the whole COVID-19 meta-interactome, thus highlighting the importance of immune processes in COVID-19 pathobiology. CAD and T2D only shared a few genes with the immune subnetwork, while, as expected, the overlap with asthma was more extended (Figure 6). Indeed, pathway and process enrichment analysis of asthma network showed the presence of immunological functions in the top clusters, in contrast to CAD (“extracellular structure organization”,”response to wounding”, “Plasma lipoprotein assembly, remodeling, and clearance”) and T2D (“Diseases of signal transduction by growth factor receptors and second messengers”, “Signaling by Receptor Tyrosine Kinases”, “PID HES HEY PATHWAY”) (Figure 7). In conclusion, comorbidities subnetworks are largely independent from the immune subnetwork, suggesting that distinct biological processes may interactively contribute to worsening COVID-19 clinical outcomes in these traits.

Figure 6. Venn diagram of the overlap between the COVID-19 immune subnetwork and comorbidities’ subnetworks, with the respective number of genes.

T2D: Type 2 Diabetes; CAD: Coronary Artery Disease.

Figure 7. COVID-19 Immune subnetwork an COVID-19 associated comorbidities subnetworks.

Each square contains the network representation, the top nodes ranked by Betweenness Centrality (if >0,00) and the pathways with the highest enrichment. Node size is proportional to betweenness centrality; node color represents densely interconnected subsets.

COVID-19 meta-interactome gene expression analysis

To assess the tissue specificity of COVID-19 meta-interactome genes, we conducted a tissue expression pattern analysis based on GTEx v8 RNA-seq data exploiting the FUMA web tool. Results showed a significantly higher expression of COVID-19 meta-interactome genes in lungs, adipose tissue, blood vessels, blood, cervix and uterus. In contrast, the lowest expression was seen in liver, pancreas, kidney and heart (Figure 8A). Furthering the analysis at the gene level, we found a more prominent tissue specificity of subsets of genes: arteries showed a marked expression of ITGA8, ERG, IRS1; AGAP2, PKIA, GNG4, KCNJ3 resulted overexpressed in the brain, in contrast to the marked underexpression of MYC, AHR, CDK2, PLAU, TNFAIP3 in this tissue. Results also highlighted a core subset of genes with a marked expression throughout the whole body, as NFKBIA, STAT3 and FURIN, a protease known to contribute to SARS-CoV-2 S protein cleavage³⁰ (Figure 8B). Overall, these expression patterns appeared to relate with the multifaceted clinical manifestation of the disease.

Figure 8. COVID-19 meta-interactome gene expression across human tissues.

A) Differentially Expressed Genes (DEGs) in GTEx v8 30 tissue types. Significantly enriched DEG sets (Pcorr < 0.05) are highlighted in red: the first graph shows DEGs overexpression, the second DEGs underespression, the third sums both side of regulation. B) Gene Expression Heatmap across GTEx v8 30 tissue types. Genes and tissues are ordered by clusters.

We then investigated if the expression pattern of COVID-19 meta-interactome genes was affected by SARS-CoV-2 infection in humans. To this aim, we reanalyzed the publicly available GSE147507 dataset, containing the transcriptional profiling of lung biopsies of a COVID-19 patient. SARS-CoV-2 infection caused an overall downregulation of COVID-19 meta-interactome genes, with the exception of a small group of genes including TNF (Figure 9A).

Recent studies suggested that SARS-CoV-2 impairs type I Interferon production, hampering the host’s antiviral response while triggering a sustained chemokine storm^10,25,31,32. To explore if the administration of IFN-β could revert the transcriptional profile of COVID-19 meta-interactome genes, we extracted from the GSE14750 the transcriptional profiling results of normal human bronchial epithelial cells treated with human recombinant IFN-β. Indeed, the COVID-19 meta-interactome genes showed a slight upregulation trend, in contrast to the effect of SARS-CoV-2, (Figure 9B). We then investigated the component of COVID-19 meta-interactome in the Interferome database, an online collection of experimental data regarding Interferon Regulated Genes (IRG): we retrieved experimental evidence that all the 248 COVID-19 meta-interactome genes are regulated by type I interferon.

Figure 9.

COVID-19 meta-interactome gene expression in SARS-CoV-2 infected human lungs (A) and normal human bronchial epithelial cells treated with IFN-Β (B). Volcano plots displaying the gene expression patterns. The y axis shows adjusted p.value; the x axis shows the Fold Change.

Drug target screening

Finally, we investigated whether COVID-19 meta-interactome nodes were suitable targets for drug repositioning. With this aim, we interrogated the Scalable Precision Medicine Oriented Knowledge Engine (SPOKE)²⁸, a continuously updated collection of biomedical data acquired from 25 databases, including CHEMBL, DrugBank and LINCS. We focused our search on the meta-hubs, i.e. the nodes with the highest Betweenness Centrality score in the whole meta-interactome and in the above analyzed subnetworks (Figure 2 and Figure 7), as their effects extend widely to the network, making them good candidates for therapeutic purposes. We only considered compounds in advanced stage of development (phase ≥3) and approved drugs, aiming at clinical actionability.

At first, we looked for drugs interacting directly with the protein derived from each gene and retrieved 10 proteins targeted by 66 different drugs (Extended data¹³) Among these appeared the JAK inhibitors Tofacitinib, approved to treat Rheumatoid Arthritis, and Fedratinib, developed for myeloproliferative diseases. This class of compounds is already under clinical investigation for COVID-19, owing to the capacity to interfere with SARS-CoV-2 infection of cells and to reduce cytokine hyperproduction³³. Anti-TNF agents, included in our results, follow the same track by blocking a pivotal molecule of the hyperinflammatory response. Among less clinically explored compounds, is Niclosamide that proved to target both LYN and STAT3 in the meta-interactome: it is an approved anti-helminthic drug, which has been recently demonstrated to have very potent antiviral activity against SARS-CoV-2³⁴.

Considering the effects of SARS-CoV-2 infection on COVID-19 meta-interactome gene expression, we extended the screening in SPOKE to drugs capable of reverting the pattern for each top node. Among these we found anti-diabetic drugs (sitagliptin, pioglitazone, rosiglitazone, canaglifozin), statins (atorvastatin, fluvastatin, simvastatin) and an angiotensin-receptor-blocker (valsartan): these results enhance the pathological relation of COVID-19 with cardiovascular diseases and diabetes, with shared pathways that may benefit from the same interventions.

Results also included antiviral agents, such as ribavirin, which is receiving positive results in COVID-19 management³⁵, and nelfinavir, which has been demonstrated to inhibit SARS-CoV-2 infection and replication both in silico and in vitro^36,37; along with these, anti-inflammatory, immunomodulatory and immunosuppressive agents also emerged. Of interest is the presence of progesterone, estradiol, selective estrogen receptor modulators (SERMs: Tamoxifen, Bazedoxifene, Raloxifene) and anti-androgen medications (finasteride, flutamide): these results suggest that the effects of diverse steroidal hormones at the transcriptional level may underlie the sex discrepancy seen in COVID-19 severity.

In conclusion, we found multiple repurposable molecules targeting SARS-CoV-2-host interactions: to organize the results, we implemented a “pleiotropy ranking” of the whole list of compounds, on the basis of the multiplicity of positive effects each drug is expected to have on COVID-19 meta-interactome both at the protein and at the gene expression level (Extended data¹³).

Underlying data

All data underlying the results are available as part of the article and no additional source data are required.

Extended data

Open Science Framework: SARS-CoV-2 meta-interactome suggests disease-specific, autoimmune pathophysiologies and therapeutic targets: extended data. https://doi.org/10.17605/OSF.IO/82T5Q¹³

Data are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).