key: cord-354096-x2skguz8
authors: Ray, Pradipta R.; Wangzhou, Andi; Ghneim, Nizar; Yousuf, Muhammad S.; Paige, Candler; Tavares-Ferreira, Diana; Mwirigi, Juliet M.; Shiers, Stephanie; Sankaranarayanan, Ishwarya; McFarland, Amelia J.; Neerukonda, Sanjay V.; Davidson, Steve; Dussor, Gregory; Burton, Michael D.; Price, Theodore J.
title: A pharmacological interactome between COVID-19 patient samples and human sensory neurons reveals potential drivers of neurogenic pulmonary dysfunction
date: 2020-06-01
journal: Brain Behav Immun
DOI: 10.1016/j.bbi.2020.05.078
sha: 
doc_id: 354096
cord_uid: x2skguz8

The SARS-CoV-2 virus infects cells of the airway and lungs in humans causing the disease COVID-19. This disease is characterized by cough, shortness of breath, and in severe cases causes pneumonia and acute respiratory distress syndrome (ARDS) which can be fatal. Bronchial alveolar lavage fluid (BALF) and plasma from mild and severe cases of COVID-19 have been profiled using protein measurements and bulk and single cell RNA sequencing. Onset of pneumonia and ARDS can be rapid in COVID-19, suggesting a potential neuronal involvement in pathology and mortality. We hypothesized that SARS-CoV-2 infection drives changes in immune cell-derived factors that then interact with receptors expressed by the sensory neuronal innervation of the lung to further promote important aspects of disease severity, including ARDS. We sought to quantify how immune cells might interact with sensory innervation of the lung in COVID-19 using published data from patients, existing RNA sequencing datasets from human dorsal root ganglion neurons and other sources, and a genome-wide ligand-receptor pair database curated for pharmacological interactions relevant for neuro-immune interactions. Our findings reveal a landscape of ligand-receptor interactions in the lung caused by SARS-CoV-2 viral infection and point to potential interventions to reduce the burden of neurogenic inflammation in COVID-19 pulmonary disease. In particular, our work highlights opportunities for clinical trials with existing or under development rheumatoid arthritis and other (e.g. CCL2, CCR5 or EGFR inhibitors) drugs to treat high risk or severe COVID-19 cases.

Caption. Our workflow, showing the different stages of RNA-sequencing, differential gene expression analysis, interactome prediction and identification of putative druggable targets, using healthy and COVID-19 BALF and healthy DRG samples.

Using an interactome-based framework we have described previously (Wangzhou et al., 2020) to find high-value pharmacologically relevant targets, we identify new, potential intervention points to reduce disease burden in patients with existing or under-development drugs. A key finding emerging from the data is that certain interventions used or under development for rheumatic or neuropathic pain might be useful against COVID-19. This finding is consistent with independent data showing that lung inflammation causes a neuropathic phenotype in neurons innervating the lung (Kaelberer et al., 2020) .

Supplementary tables can be found at our online repository, as Excel sheets: https://www.utdallas.edu/~prr105020/covid19/supp_tables.zip .

We first re-analyzed bulk RNA-seq data from published studies by mapping and quantifying relative gene abundance (in counts per million) from BALF (National Genomics Data Center accession number: PRJCA002326, https://bigd.big.ac.cn) of severe COVID-19 patients (Xiong et al., 2020b) , compared to BALF (NCBI SRAdb project SRP230751, https://www.ncbi.nlm.nih.gov) of healthy controls. The BALF method samples fluid from the lower respiratory tract and is commonly used as a diagnostic tool in ARDS. The method removes fluid, mucus and cells from the lung. This allows for RNA sequencing of immune and other cells found in the lung. Our healthy control samples had highly consistent RNA profiles (Pearsons' R = 0.9263). Differences between the two COVID-19 BALF samples can likely be attributed to variability of the immune response from patient to patient, as well as the lower sequencing depth of the samples, as even after pooling technical replicates for the COVID-19 patient samples we had ~ 12M and 7M sequenced reads. The two COVID-19 BALF samples were also found to contain SARS-CoV2 viral RNA reads confirming the presence of disease (Xiong et al., 2020b) .

Out of 31,069 genes in the reference genome, we detected 18,507 and 18,855 genes in the control samples, and 13,973 and 13,545 genes in the COVID-19 patient samples; differences being primarily attributable to a decrease in sequencing depth for COVID-19 samples. Due to this and the fact that our goal was to identify potentially targetable protein interaction pathways in COVID-19 patients, we performed downstream analysis only on 1372 genes that were upregulated in the COVID-19 BALF samples (Supplementary Table 1) .

Genes upregulated in the COVID-19 patient samples include a multitude of genes that recapitulate clinical characteristics, such as increased cytokine signaling causing a "cytokine storm" (including CCL2/3/4/7/8, CXCL1/2/6/8/28, CCL3L3, CCL4L2) , hypoxia (HIF1A, HLF), and inflammasome and sepsis-related genes (IL1R1/2, IL5RA, IL33, IL31RA). BALF typically contains RNA profiles of immune cells in alveolar fluid, primarily lymphocytes and macrophages. Additionally, we found that upregulation of transcription factor genes in COVID-19 samples identifies transcription factors associated with alveolar cell types (EHF, PAX9, ELF3, GHRL2) and immune cells (RFX3, SOX5, TP63, HOPX) with functions including regulation of antiviral pathways (NR3C2), based on ARCHS4 database (Lachmann et al., 2018) and the Enrichr gene set enrichment analysis tool (Kuleshov et al., 2016) (Supplementary Table 2 ). This suggests that cellular distress and virus-driven cell lysis causes lung cell mRNA content to be detected in COVID-19 BALF samples. Immune cell markers like CD4, CD8, and CD68 were detected in COVID-19 samples, but were not consistently upregulated. Since we only quantified relative (and not absolute) abundance, this is possibly due to increased proportion of alveolar cell mRNA content in disease samples. However, lymphocyte markers CD24 and CD38 were upregulated in the COVID-19 BALF samples, which is also consistent with reports of lymphocyte induced apoptosis in COVID-19 (Huang et al., 2020a) . These findings are congruent with reports of increased lactate dehydrogenase as a biomarker of COVID-19 severity wherein lactate dehydrogenase is a sign of pyroptosis -inflammasome driven programmed cell death (Rayamajhi et al., 2013; Chen et al., 2020; Han et al., 2020) .

We were primarily interested in ligand to receptor signaling from immune cells in the lung to lung-innervating sensory neurons. For this purpose, we used our interactome framework (Ramilowski et al., 2015) to elucidate possible pharmacological interactions that may be associated with disease in COVID-19 samples. We focused first on ligands upregulated in the COVID-19 samples and receptors expressed in human thoracic DRG samples, which contain neurons that innervate the lungs (Springall et al., 1987; Kummer et al., 1992) . We found that many receptors for cytokines identified as upregulated in COVID-19 patients including CCR2, CXCR2, CCR5, and CXCR10, and IL15RA were also expressed in hDRG suggesting a potential direct connection between these cytokines and sensory neuron activation in the lung ( Table 1) . We also noted increased EREG expression in COVID-19 samples, which is known to signal via the epidermal growth factor receptor (EGFR) to sensitize nociceptors (Martin et al., 2017) . Interestingly, there were a number of Eph ligand genes (EFNA1 and EFNA5) whose gene products are known to act via Eph receptors to modulate the activity of neurons as well as regulating connections between neurons (Henderson and Dalva, 2018) . These Eph ligands were expressed de novo in disease samples. A complete list of the interactome featuring differentially increased ligands in BALF from COVID-19 samples can be seen in Supplementary Table 3 , and expression profiles of the mouse orthologs in sensory neurons can be found in Supplementary Table 4 .

Previous work has shown that nociceptor-derived CGRP plays a key role in dampening the immune response to bacterial infection in the lung (Baral et al., 2018) . Therefore, we also explored possible interaction points between sensory neurons that are likely to innervate the lung and receptors expressed in BALF samples from patients with COVID-19. Here, we identified a striking number of interactions between Eph ligands found in hDRG and receptors upregulated in COVID-19 BALF samples ( Table 2) . This is meaningful as there is emerging literature on an important role of the ephrin -Eph system in immunity (Darling and Lamb, 2019) . Our data suggests potential bidirectional interactions between these receptors and ligands on sensory nerve endings and immune cells in the lung of COVID-19 patients. Another prominent interaction was with EGFR ligands and EGFR itself, again suggesting bidirectional interactions between neurons and immune cells in the lung in COVID-19. Finally, we also noted potential interactions between neurexins and neuroligins, which are also involved in junctions between neurons (Levinson and El-Husseini, 2005) , suggesting that multiple types of mediators may be involved in remodeling nerve ending morphology within the lung driven by immune cell activation in COVID-19. Collectively, these neuron-derived ligands could potentially exacerbate lung inflammation in COVID-19 contributing to a positive feedback loop. A complete list of the interactome featuring differentially increased receptors in BALF from COVID-19 samples can be seen in Supplementary Table 5 and expression profiles of the mouse orthologs in sensory neurons can be found in Supplementary Table 6 .

The immune response is tuned by both harnessing specific cell types (by chemotaxis and other mechanisms), and by modulating molecular profiles of these cell types (by transcriptional, post-transcriptional, translational and post-translational mechanisms). While our interactome identification from the BALF bulk RNA-seq is a starting point for identifying neuro-immune interactions in COVID-19 patient lungs, bulk RNA sequencing has limitations and other approaches are also being used to characterize the immune response in this disease. We therefore studied the emerging COVID-19 literature for studies where raw datasets were not publicly available, but the studies provided gene or protein sets that were associated with clinical outcomes or pathologies. In addition to the Xiong et al study (Xiong et al., 2020b) , we integrated gene and protein lists from sources that included BALF single-cell RNA-seq (scRNA-seq), blood protein levels and immune cell profiles from COVID-19 patients (Liao et al., 2020; Xiong et al., 2020b; Zhou et al., 2020a; Zhou et al., 2020b) . This effort identified more than 150 genes. Since a large proportion of these were ligands, we constructed a ligand-receptor interactome from the ligand identified in the studies cited above to hDRG expressed receptors based on our previous RNA-seq experiments . These findings were consistent with the interactome generated from BALF samples from COVID19 patients with prominent interactions for cytokines, chemokines and other inflammatory molecules (Supplementary Table 7 ). These included CCL2/3/4/5/7/10, IL1 and TNF as well as anti-inflammatory molecules IL10 and TGFB1 and TGFB2. These finding implicate the potential off-label utility of anti-rheumatoid arthritis drugs that target TNF (infliximab or etanercept) or IL1 (anakinra) as well as EGFR inhibitors such as cetuximab for interfering with neuro-immune interactions in COVID-19.

Using only deeply sequenced thoracic hDRG bulk RNA-sequencing as a proxy for human nociceptor expression can generate both false positives and false negatives. Bulk hDRG sequencing libraries incorporate mRNA from not just sensory neurons, but also from glial, immune and vascular cells, leading to the possibility of false positives. Additionally, many lung innervating nociceptors are from the nodose or jugular ganglia so our dataset from thoracic hDRG may not represent the full possibility of receptors or ligands that could contribute to lung physiology from sensory afferents. Unfortunately, no sequencing data is available from the human nodose or jugular ganglia. To predict which of the hDRG-expressed ligands or receptors are likely present in human nociceptors, we integrated several mouse datasets. We used FACS sorted neuronally enriched mouse DRG cell pool RNA-seq (Liang et al., 2019) , and mouse DRG scRNA-seq data (Usoskin et al., 2015) . We additionally included two datasets for gene expression in the mouse jugular-nodose complex (JNC), including the Kupari et al scRNA-seq dataset (Kupari et al., 2019 ) and a study using bulk RNA-seq to study gene expression changes in the JNC following lung inflammation with Lipopolysaccharide (LPS) (Kaelberer et al., 2020) . From this analysis (Supplementary Tables 3-6) we found clear evidence for nociceptor expression of EFNA1, EFNB1 and EFNB2 suggesting the ephrin signaling is likely to influence lung nociceptors. The same was true for neuroligin genes NLGN2 and NLGN3. Likewise, we found strong evidence for cytokine receptor gene expression in nociceptors, most notably CCR2, LIF (which was specific to the nonpeptidergic population) and IL15RA but also ITGAM and ITGA9. As mentioned above, EGFR ligands were strongly induced in COVID-19 patients. The ERBB2 gene was strongly expressed in single cell and sorted neuron samples but the EGFR gene was sparsely expressed or not detected. However, EGFR was increased in nerve injury samples suggesting that this receptor is induced in nociceptors after injury. Since ERBB2 encodes a receptor that acts in concert with EGFR, more work is needed to better understand how this receptor may function to modulate the activity of lung nociceptors.

A major focus of current research into SARS-CoV-2 and COVID-19 is identification of potential therapeutics. From that perspective, a recent paper described a protein-protein interactome of SARS-CoV-2 proteins with the human proteome and identified druggable interaction points (Gordon et al., 2020) . A target identified from that paper was the RNA binding protein, eIF4E, the phosphorylation of which has been implicated in coronavirus replication (Banerjee et al., 2002; Mizutani et al., 2004) and the suppression of anti-viral responses (Herdy et al., 2012) . Previous work suggests that eIF4E phosphorylation by MNK, a druggable target identified by Gordon et al. (2020) , plays a key role in regulating the excitability of nociceptors (Moy et al., 2017; Megat et al., 2019; Shiers et al., 2020) and also regulates the translation of many cytokines and chemokines (Furic et al., 2010; Herdy et al., 2012; Amorim et al., 2018) . Based on this, we intersected differentially expressed BALF ligand and receptor genes in COVID-19 samples with mRNAs known to be regulated by eIF4E phosphorylation from previous studies (Furic et al., 2010; Aguilar-Valles et al., 2018; Amorim et al., 2018) . To do this we compiled a list of genes shown to be significantly changed in COVID-19 patients from the studies we cited and described above (Supplementary Tables 3, 5 and 7) . We also compiled a list of mRNAs whose translation has been shown to be regulated by eIF4E phosphorylation based on the literature (Furic et al., 2010; Aguilar-Valles et al., 2018; Amorim et al., 2018) . We then looked for genes that were common to these two lists.

We found that 8 inflammatory mediators upregulated in the BALF of COVID-19 patients were reduced in vitro and/or in the brain of animals lacking eIF4E phosphorylation (Table 3) . CCL2, which is among the most highly upregulated chemokines in COVID-19 patients, is one example. CCL2 impairs virus clearance, reduces macrophage maturation, and increases lethality after coronavirus infection (Held et al., 2004; Trujillo et al., 2013) . Furthermore, CCL2 levels in the nasal aspirate of asthmatic individuals are correlated with severe respiratory symptoms following viral infection (Lewis et al., 2012) . Upregulation of CCL2 following SARS-CoV (the virus that caused the first SARS epidemic) infection is mediated by Ras-ERK pathway that is activated by ACE2 signaling (Chen et al., 2010) . A downstream target of the Ras-ERK pathway is MNK, which in turn phosphorylates eIF4E to influence mRNA translation. It is possible that the novel coronavirus, SARS-CoV-2, which also enters cells through ACE2, enhances the production of CCL2 via MNK-eIF4E signaling. Hence, eFT508 (also known as Tomivosertib), a MNK1/2 inhibitor currently undergoing clinical trials (Reich et al., 2018) , may prove beneficial against COVID-19. The expected mechanistic outcome of eFT508 treatment would be decreased translation of many chemokines and cytokines that are predicted to act on sensory innervation of the lung. Therefore, this particular therapeutic approach could decrease pathological immune-neuronal interactions without completely sequestering individual immune mediators.

We also identified potential therapeutics for COVID-19 by intersecting the interactomes generated here with the DGIdb database (Cotto et al., 2018) . We identified 144 gene products with characterized antagonists, inhibitors, blockers or modulators in this curated database (a selection of these are shown in Table 4 , full dataset in Supplementary Table 8 ). This database highlights the potential for targeting cytokines and chemokines, such as CCL2, the CCR2 receptor and/or the EGFR receptor in COVID-19. It also reveals several drugs that interact with ephrin -EphB signaling that may have promise if structural remodeling of lung afferent endings contributes to COVID-19 pathology. Interestingly, the NMDA receptor 2B (GRIN2B) subunit was revealed in this list, which can be targeted with ketamine and other drugs. The corresponding gene GRIN2B was strongly induced in BALF samples of COVID-19 samples. The NMDA receptor is activated by glutamate, which is released by sensory neurons at their peripheral terminals. NMDA receptor activity and localization is strongly influenced by EphB receptors (Hanamura et al., 2017; Henderson and Dalva, 2018) , which were also dramatically increased in COVID-19 BALF samples. These findings reveal a potential new pathway for sensory modulation of lung pathophysiology involving glutamate from sensory afferents acting on NMDA receptors expressed by resident cells in the lung.

We used a computational framework to test the hypothesis that immune cells in the lungs of severe COVID-19 patients produce factors that are likely to interact with the sensory neuronal innervation of the lung. Our work finds support for this hypothesis and identifies potential interactions that may be critical drivers of disease severity in COVID-19. Some of these, such as CCL2, epiregulin acting via the EGFR receptor and TNF can be inhibited with drugs in late stage development (e.g. CCL2) or existing therapeutics (e.g. EGFR or TNF) (summarized in Figure 2 ). All 3 of these targets have been implicated in increasing nociceptor excitability (Hensellek et al., 2007; Belkouch et al., 2011; Martin et al., 2017) . It may be necessary to target multiple pathways to overcome the severe respiratory phenotype produced by this disease. To that end, MNK inhibitors, which were also identified by Gordon et al. (Gordon et al., 2020) may be particularly useful because they have an effect on coronavirus replication in murine models (Banerjee et al., 2002) , interfere with the translation of chemokine and cytokines implicated in COVID-19 and reduce excitability of nociceptors (Moy et al., 2017; Megat et al., 2019) that may be driving neurogenic inflammation in the disease. Clinical trials are obviously needed to test these hypotheses. Endpoints for such trials should likely include oxygen saturation, perceived shortness of breath and pneumonia severity. Many of the candidate therapeutics we have identified could have potential impacts of immune status (for instance TNF sequestering therapeutics are immune suppressant) so this must be considered from a safety standpoint.

There are important caveats to our work. First, the interactomes we have built are based largely on data from a limited number of patients. Similar results have been found across studies (Huang et al., 2020b; Liao et al., 2020; Xiong et al., 2020b; Zhou et al., 2020a; Zhou et al., 2020b) , but there is an obvious need for additional sequencing on lung immune cells from patients with mild and severe COVID-19. The availability of such data, either bulk RNA-seq or single cell RNAseq would greatly improve our ability to hone in on more specific targets.

Our work identifies that: 1) Several proinflammatory cytokines and chemokines are upregulated in the COVID-19 BALF samples that are known to be translationally-regulated via MNK-eIF4E signaling. This offers a unique opportunity to disrupt activity of many inflammatory proteins via MNK inhibitors. 2) Many upregulated COVID-19 inflammatory mediators interact with receptors potentially expressed on thoracic sensory neurons that innervate the lung. Activation of these sensory neurons may cause them to release neuropeptides back into the lung environment to cause vasodilation, immune cell recruitment, neurogenic inflammation, and potentially even pain upon breathing.

Another important caveat is that the directionality of the influence of nociceptor-released factors on lung immunity and COVID-19 disease state is not currently clear. Some studies on sepsis, which occurs in severe COVID-19 patients, have shown that nociceptors promote mortality in sepsis (Bryant et al., 2003) but others show that TRPV1-positive nociceptors are protective against mortality in sepsis (Guptill et al., 2011) . In the case of lung infection by bacteria, CGRP, which is released by nociceptors, has been shown to play a protective role by limiting immune reaction within the lung (Baral et al., 2018) . Little work has been done to understand how viral infections and sensory neurons interact to promote or suppress the immune response. However, the existing literature on sensory innervation of the airway and lung with bacterial pathogens suggests that sensory afferents can have detrimental and beneficial effects depending on the context (Chiu et al., 2012; Chavan et al., 2018; Ruhl et al., 2020) . We have interpreted our interactome findings under the assumption that neurogenic responses may worsen the disease state in COVID-19, potentially leading to ARDS and fatality, but without perturbational animal model studies or additional human clinical data we cannot rule out the possibility that a subset of the interactome mitigates the disease state. However, recent findings that lung inflammation creates a neuropathic-like state in lung-innervating nociceptors (Kaelberer et al., 2020) suggests that our hypothesis of neurogenic response worsening the disease is supported by available data.

A third caveat is the fact that we have predicted neuro-immune interactions based on the assumption that lung and airway-innervating sensory neuronal transcriptomes remain mostly unchanged after infection. It is formally possible that sensory neurons could express ACE2 allowing viral infection of these cells, but available datasets are unable to clarify this question. A bigger potential issue is that lung inflammation clearly alters the transcriptome of neurons that innervate the lung in rodents (Kaelberer et al., 2020) , but the time-course and magnitude of this effect at the single cell level has not been studied previously, and little is known about evolutionary conservation of the transcriptional response between humans and rodents. These are obviously important areas of research that need to be explored in future studies.

In conclusion, our work identifies several potential therapeutics that target neuroimmune interactions in the lung that may be useful for treatment of pulmonary pathology such as ARDS in COVID-19 patients. Based on the currently available evidence, we suggest further investigation of MNK and EGFR inhibitors, TNF sequestering treatments and CCL2 inhibitors, alone or in combination, for the treatment of pulmonary pathology in COVID-19. There are safety concerns that must be considered as some of these therapeutics (in particular TNF) decrease immunity and may have a negative influence on patient outcomes. (Furic et al., 2010; Huang et al., 2020a; Liao et al., 2020; Xiong et al., 2020b) CCL7 C-C motif chemokine ligand 7 ↑ ↓ (Furic et al., 2010; Liao et al., 2020; Xiong et al., 2020b) Differentially expressed genes in BALF of COVID-19 patients that are known to be regulated by MNK -eIF4E signaling based on studies in eIF4E S209A -knock-in mutant cells and animals. These knock-in animals and cells are null mutants for eIF4E phosphorylation. Downregulation or upregulation with respect to eiF4E signaling refers to genes that have decreased or increased translation respectively in the absence of MNK-eIF4E signaling. Table 4 Gene product Drug (mechanism)

Carlumab ( 

IRB approval for RNA sequencing from human DRG samples was provided by University of Texas at Dallas as described previously . Other datasets used here were based on analysis from datasets described in published studies from other groups.

Human thoracic DRG data was described in . Data for the curated COVID-19 interactome were obtained from (Huang et al., 2020b; Liao et al., 2020; Xiong et al., 2020b; Zhou et al., 2020a) , with BALF RNA-seq data from COVID-19 patients (National Genomics Data Center sample ids CRR119894-7, https://bigd.big.ac.cn/bioproject/browse/PRJCA002326 ) obtained from (Xiong et al., 2020b) . BALF control data was obtained from control RNA-seq samples (Michalovich et al., 2019) in healthy patients (SRAdb sample ids SRR10571724, and SRR10571732, https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=SRP230751).

Mouse DRG single cell sequencing data were obtained from mouse DRG (Usoskin et al., 2015) and JNC (Kupari et al., 2019) . Mouse bulk RNA-sequencing datasets were obtained from sorted neuronal cell pools from the DRG (Liang et al., 2019) and whole tissue RNA-seq from the JNC complex before and after LPS injury (Kaelberer et al., 2020) .

Three previous studies of eIF4E phosphorylation (Furic et al., 2010; Aguilar-Valles et al., 2018; Amorim et al., 2018) were used to identify the set of genes which can be potentially affected by MNK inhibition. The DGIdb database was used to identify known drugs (Cotto et al., 2018) that target proteins in our identified interactomes.

The interactome database is described in detail in our previous work (Wangzhou et al., 2020) .

BALF RNA-seq data (Xiong et al., 2020b) was compared with control RNA-seq by mapping RNA-seq reads to the reference genome Refseq hg38 and reference transcriptome Gencode v33 (Frankish et al., 2019) using the STAR aligner (Dobin et al., 2013) . Read counts per gene was then summarized for each sample by HT-Seq Count (Anders et al., 2015) . Replicates from each COVID-19 patient were pooled to increase sequencing depth. Due to the higher sequencing depth of the control samples (over 40 million reads), compared to the COVID-19 samples (~ 7 and 12 million reads after pooling (Xiong et al., 2020b) ), and due to specific interest in finding therapeutic targets, only genes identified as differentially expressed and increased in relative abundance in COVID-19 samples were further analyzed. Low replicate numbers cause difficulty in estimation of dispersion, and lower sequencing depth increases sampling variance in the data. To conservatively identify differentially expressed genes, analysis was performed using two well established differential expression analysis tools -edgeR (Robinson et al., 2010) and DEseq2 (Love et al., 2014) , both of which model gene counts as negative binomial distributions and estimate mean and dispersion parameters from the data. edgeR moderates estimated gene count dispersions towards a common sample-specific estimate for all genes. DEseq2 models dispersion as a function of the expression level. Genes that have an adjusted p-value of < 0.001 and log2 fold change > 2.0 were used to identify differentially expressed gene sets for each tool. Differentially expressed gene sets from both tools were intersected to identify genes that are differentially expressed in both models, and only genes with read frequency > 1 in 200,000 genic reads for COVID-19 samples were retained to reduce potential effects of sampling variance. For differentially expressed genes, we provide relative abundance in the form of counts per million mapped genic reads (CPM) for each sample, as well as p-values and multiple testing adjusted pvalues for both tools (Supplementary Table 1) . Gene set enrichment analysis was performed using the Enrichr framework, and limited to gene -transcription factor co-expression modules with 50 or more genes (Supplementary Table 2) .

Ramilowski et al. (Ramilowski et al., 2015) described 2557 pairs of ligand-receptor interactions that were used to populate the initial interactome database. However, hundreds of known ligands and receptors were missing from this resource. In order to curate a more complete ligand-receptor list, we collected gene lists from gene family and ontology databases HUGO (Yates et al., 2017) and GO (Carbon et al., 2009) . Genes corresponding to the GO terms of cell adhesion molecules (CAM), G-protein coupled receptor (GPCR), growth factor, ion channels, neuropeptides, nuclear receptors, and receptor kinases were tabulated. Based on a database (Szklarczyk et al., 2019) and the literature, corresponding receptors or ligands for these additional genes were identified and added to the list of ligand-receptor pairs. A total of 3098 pairs of ligandreceptor interactions were used in the interactome analysis for all possible ligand-receptor interactions (Wangzhou et al., 2020) .

Not all ligands and receptors are encoded as genes, so we included some additional records in our database, such as enzymes known to synthesize ligands and paired with the corresponding receptor. Additionally, some other interactions were included that do not have a traditional ligand -receptor relationship. In these cases, the gene for the upstream signaling molecule was analyzed in lieu of the ligand, and the gene for the downstream signaling molecule in lieu of the receptor. The entire interactome database is described in detail in (Wangzhou et al., 2020) .

Based on gene expression in BALF and human DRG samples, the interactome extends to thousands of candidate ligand-receptor pairs. While many of these signaling interactions are likely involved in the immune response of COVID-19, we chose to focus on pairs where the ligand or receptor was shown to be differentially increased or discriminative marker for cell types with differentially increased proportions in the immune response since these are more likely to be successful points of therapeutic intervention. Thus, we filtered the candidate interactome based on the list of differentially expressed genes we obtained from the literature and re-analysis of COVID-19 BALF data (Supplementary Tables 3, 5 and 7) . Potentially relevant gene product -drug interactions (inhibitor, antagonist or blocker) based on these interactomes were mined from the dgiDB database (Supplementary Table 8) .

Finally, we decided to tabulate these interactions with the ability to rank them using four metrics reported in Supplementary Tables. Two of these help rank immune genes : mean expression level in COVID-19 samples, and degree of differential expression (fold change) in COVID-19 versus control BALF (Supplementary Tables 3 and 5) . These are likely to identify signaling pathways with the greatest frequency of ligand-receptor interaction, or maximal changes in the frequency or intensity of such interactions. Two other metrics help rank neuronal genes: gene expression levels in sorted mouse DRG neuron datasets (Liang et al., 2019) and in mouse DRG neuron scRNA-seq datasets (Usoskin et al., 2015) , and degree of differential expression (log fold change) in LPS based injury models of mouse airway-innervating neurons (Kaelberer et al., 2020) (Supplementary Tables 4 and 6) . These are likely to identify neuronally expressed genes, as well as genes involved in ARDS in mammalian lungs. The directionality of changes in the LPS mouse injury model may not necessarily be consistent with changes in COVID-19 molecular changes, but are likely to be indicative of genes involved in neuroinflammation in mammalian lungs.

In order to identify genes that were stably expressed in the hDRG, only protein interactions where the corresponding neuronal gene expression in the hDRG was >= 0.15 TPM (as reported in North et al., 2019) were retained. To quantify relative gene abundances in (Liang et al., 2019) , reported coding gene expressions were renormalized to TPMs across all samples, and samples having a clear bimodal distribution of TPMs were retained for further analysis and quantile normalized, of which the FACS-sorted neuronal datasets are used in our study for integrative analysis. Reported gene abundances in (Kaelberer et al., 2020) were renormalized to a million. In the (Usoskin et al., 2015) dataset, the reported fraction of cells for each neuronal subpopulation (where a gene is detected) was used. To calculate a similar metric in the (Kupari et al., 2019) dataset, two steps were performed. Based on the reported clusterings of the cell types in the JNC, the four most distinct groups of neuronal subpopulations (Kupari et al., 2019) , comprising jugular ganglial neuronal subpopulations 1-3 and 4-6, and nodose ganglial neuronal subpopulations 1-11 and 12-18 were used. To calculate the fraction of cells where genes were detected in each of these four groups, cells in the lowest quartile (bottom 25 th percentile) of the cell population with respect to library diversity (number of unique molecular identifiers or UMI < 23,385) were discarded.

Translational control of depression-like behavior via phosphorylation of eukaryotic translation initiation factor 4E

Loss of eIF4E Phosphorylation Engenders Depressionlike Behaviors via Selective mRNA Translation

HTSeq-a Python framework to work with high-throughput sequencing data

A key role for gp130 expressed on peripheral sensory nerves in pathological pain

Murine coronavirus replication-induced p38 mitogen-activated protein kinase activation promotes interleukin-6 production and virus replication in cultured cells

Nociceptor sensory neurons suppress neutrophil and gammadelta T cell responses in bacterial lung infections and lethal pneumonia

The chemokine CCL2 increases Nav1.8 sodium channel activity in primary sensory neurons through a Gbetagamma-dependent mechanism

Transient receptor potential cation channel, subfamily V, member 4 and airway sensory afferent activation: Role of adenosine triphosphate

Capsaicin-sensitive nerves regulate the metabolic response to abdominal sepsis

Neurology of allergic inflammation and rhinitis

Functional implications of the multiple afferent pathways regulating cough

Neural regulation of airway smooth muscle tone

Sensory nerves and airway irritability

AmiGO: online access to ontology and annotation data

Neuro-immune interactions in inflammation and host defense: Implications for transplantation

Clinical and immunologic features in severe and moderate Coronavirus Disease

Upregulation of the Chemokine (C-C Motif) Ligand 2 via a Severe Acute Respiratory Syndrome Coronavirus Spike-ACE2 Signaling Pathway

Neurogenic inflammation and the peripheral nervous system in host defense and immunopathology

DGIdb 3.0: a redesign and expansion of the drug-gene interaction database

Emerging Roles for Eph Receptors and Ephrin Ligands in Immunity

Role of the transient receptor potential vanilloid 1 in inflammation and sepsis

STAR: ultrafast universal RNA-seq aligner

Nociceptors: the sensors of the pain pathway

GENCODE reference annotation for the human and mouse genomes

phosphorylation promotes tumorigenesis and is associated with prostate cancer progression

BLU-5937: A selective P2X3 antagonist with potent anti-tussive effect and no taste alteration

A SARS-CoV-2-Human Protein-Protein Interaction Map Reveals Drug Targets and Potential Drug-Repurposing

Disruption of the transient receptor potential vanilloid 1 can affect survival, bacterial clearance, and cytokine gene expression during murine sepsis

Sensory nerve terminal mitochondrial dysfunction induces hyperexcitability in airway nociceptors via protein kinase C

Lactate dehydrogenase, a Risk Factor of Severe COVID-19 Patients. medRxiv:2020

Extracellular phosphorylation of a receptor tyrosine kinase controls synaptic localization of NMDA receptors and regulates pathological pain

Differential roles of CCL2 and CCR2 in host defense to coronavirus infection

EphBs and ephrin-Bs: Trans-synaptic organizers of synapse development and function

The cytokine TNFalpha increases the proportion of DRG neurones expressing the TRPV1 receptor via the TNFR1 receptor and ERK activation

Translational control of the activation of transcription factor NF-kappaB and production of type I interferon by phosphorylation of the translation factor eIF4E

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Blood single cell immune profiling reveals the interferon-MAPK pathway mediated adaptive immune response for

Activation of a nerve injury transcriptional signature in airway-innervating sensory neurons after Lipopolysaccharide induced lung inflammation

Enrichr: a comprehensive gene set enrichment analysis web server 2016 update

The sensory and sympathetic innervation of guinea-pig lung and trachea as studied by retrograde neuronal tracing and doublelabelling immunohistochemistry

An atlas of vagal sensory neurons and their molecular specialization

Race to find COVID-19 treatments accelerates

Massive mining of publicly available RNA-seq data from human and mouse

Building excitatory and inhibitory synapses: balancing neuroligin partnerships

Nasal cytokine responses to natural colds in asthmatic children

A transcriptional toolbox for exploring peripheral neuro-immune interactions

The landscape of lung bronchoalveolar immune cells in COVID-19 revealed by single-cell RNA sequencing

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Epiregulin and EGFR interactions are involved in pain processing

Nociceptor Translational Profiling Reveals the Ragulator-Rag GTPase Complex as a Critical Generator of Neuropathic Pain

Obesity and disease severity magnify disturbed microbiome-immune interactions in asthma patients

Phosphorylation of p38 MAPK and its downstream targets in SARS coronavirus-infected cells

The MNK-eIF4E Signaling Axis Contributes to Injury-Induced Nociceptive Plasticity and the Development of Chronic Pain

Afferent neural pathways mediating cough in animals and humans

Electrophysiological and transcriptomic correlates of neuropathic pain in human dorsal root ganglion neurons

A draft network of ligand-receptor-mediated multicellular signalling in human

Comparative transcriptome profiling of the human and mouse dorsal root ganglia: an RNA-seq-based resource for pain and sensory neuroscience research

Detection of pyroptosis by measuring released lactate dehydrogenase activity

Structure-based Design of Pyridone-Aminal eFT508 Targeting Dysregulated Translation by Selective Mitogen-activated Protein Kinase Interacting Kinases 1 and 2 (MNK1/2) Inhibition

edgeR: a Bioconductor package for differential expression analysis of digital gene expression data

Mycobacterium tuberculosis Sulfolipid-1 Activates Nociceptive Neurons and Induces Cough

Efferent functions of C-fiber nociceptors

Reversal of peripheral nerve injury-induced neuropathic pain and cognitive dysfunction via genetic and tomivosertib targeting of MNK

Retrograde tracing shows that CGRP-immunoreactive nerves of rat trachea and lung originate from vagal and dorsal root ganglia

STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets

Silencing Nociceptor Neurons Reduces Allergic Airway Inflammation

Potent binding of 2019 novel coronavirus spike protein by a SARS coronavirus-specific human monoclonal antibody

Transgenic CCL2 Expression in the Central Nervous System Results in a Dysregulated Immune Response and Enhanced Lethality after Coronavirus Infection

Unbiased classification of sensory neuron types by large-scale single-cell RNA sequencing

Calcitonin Gene-Related Peptide Negatively Regulates Alarmin-Driven Type 2 Innate Lymphoid Cell Responses

Receptor Recognition by the Novel Coronavirus from Wuhan: an Analysis Based on Decade-Long Structural Studies of SARS Coronavirus

A pharmacological interactome platform for discovery of pain mechanisms and targets

Nociceptors--noxious stimulus detectors

Transcriptomic Characteristics of Bronchoalveolar Lavage Fluid and Peripheral Blood Mononuclear Cells in COVID-19 Patients

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients

Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2

Genenames.org: the HGNC and VGNC resources in 2017

Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study

Highlights for : A pharmacological interactome between COVID-19 patient samples and human sensory neurons reveals potential drivers of neurogenic pulmonary dysfunction

 Our paper identifies key neuro-immune interactions in the lungs of COVID-19 patients that potentially drive neurogenic pulmonary dysfunction  Our framework uses an interactome prediction model to identify relevant ligandreceptor signaling between human sensory neurons and bronchial alveolar lavage fluid from COVID-19 patients  Our findings reveal a landscape of ligand-receptor interactions in the lung caused by SARS-CoV-2 viral infection and point to potential interventions to reduce the burden of neurogenic inflammation in COVID-19 disease  Our work highlights opportunities for clinical trials with existing or under development rheumatoid arthritis and other (e.g. CCL2, CCR5 or EGFR inhibitors) drugs to treat high risk or severe COVID-19 cases