Transcriptional neural-like signaling contributes to an immune-suppressive tumor microenvironment.

Tumor innervation has recently been documented and characterized in various settings and tumor types. However, the role that nerves innervating tumors play in the pathogenesis of cancer has not been clarified. In this study, we searched for neural signaling from bulk RNA sequencing from The Cancer Genome Atlas (TCGA) dataset and looked for patterns of interactions between different cell types within the tumor environment. Using a presynapse signature (PSS) as a probe, we showed that multiple stromal cell types crosstalk and/or contribute to neural signals. Based on the correlation and linear regression, we hypothesized that neural signals contribute to an immune-suppressive tumor microenvironment (TME). To test this hypothesis, we performed in vitro dorsal root ganglion (DRG)/macrophage coculture experiments. Compared to the M2 macrophage monoculture, the DRG/M2 macrophage coculture prevented anti-inflammatory M2 to pro-inflammatory M1 polarization by LPS stimulation. Finally, a survey of different TCGA tumor types indicated that higher RNA neural signature is predictive of poor patient outcomes in multiple tumor types.

INTRODUCTION

Cancer is the second most prevalent cause of death (following heart diseases) and responsible for ~20% of all deaths in the United States in 2018. Breakthroughs in cancer treatment and patient care led to a 27% reduction of cancer death rate in last 10 years (https://www.cdc.gov/cancer/dcpc/research/update‐on‐cancer‐deaths/index.htm). While the discovery of new biology and targeted therapies were the main drivers of this reduction, an improved understanding of the pathogenesis of tumor could be a force to further reduce the death rate. To accelerate discoveries, various big data cancer consortiums have been established, and multiple databases have been made publicly available for testing hypotheses and mining for novel therapeutic targets. Among them, the Cancer Genome Atlas (TCGA) is one of the biggest cancer omics programs that compiles over 20 thousand primary cancer and matched normal samples from 11 thousand patients spanning 33 cancer types. This comprehensive atlas contains clinical outcomes, biospecimens, pathology reports, single nucleotide polymorphisms (SNP) microarray, whole exome sequencing, diagnostic images, DNA methylation, bulk RNA sequencing, and reverse‐phase protein arrays (RPPA) of all samples. Since its availability, many high‐profile bioinformatics papers have been published based on TCGA dataset analysis. Hypotheses have been tested, pathways surveyed, and targets identified using this database. Combinations of TCGA with other databases and deconvolution algorithms such as xCell, Cibersort, , and ESTIMATE have further brought new ways to deconvolute the cell type fractions in the bulk tumor samples based on cell type‐specific gene signatures, which enables certain studies on tumor microenvironment.

Long before different cell types could be estimated through the tumor transcriptomics analysis, it had been widely accepted that tumor tissues contain not only cancer cells, but also nonmalignant cells such as cancer‐associated fibroblasts (CAFs), tumor vasculature endothelial cells, adipocytes, glia cells, and immune cells that either reside in or proximal to tumors. All these cancer and non‐cancer cells in the pathological tumor site are part of the tumor microenvironment (TME). The composition of the cell types in the TME not only could impact the efficacies of cancer drug responses, but could also be utilized to treat cancers such as in the case of anti‐VEGF and immunotherapies. ,

In recent years, an emerging field of cancer neuroscience or exoneural biology more broadly has shown that nerves are frequently present in tumor sites and are not only a passive structure in the TME, but also play active roles in promoting or inhibiting cancer progression. , In the central nervous system, it was found that glutamatergic neurons increase glioma cell proliferation by secreting growth factors as well as electrically communicate with the glioma cells. , , This mechanism is likely to be utilized not only by brain‐originated cancer cells, but also the other cancer types that have metastasized to the brain. Outside of the blood–brain barrier, the peripheral nervous system coordinates angiogenesis and perineural invasion. , Recently, new mechanisms of direct nerve‐tumor interactions have also been discovered by using various neuronal technologies such as neural modulations with toxins, surgical denervation, optogenetics, and chemogenetics. , , , , Evidence of peripheral nerve system (PNS) direct contributions to tumors has been found in pancreatic, prostate, , breast, head and neck, lung, cervical, , and gastric/colon cancers and could involve somatosensory, autonomic as well as vagal nerves.

Neurons use synapses (pre‐ and post‐) to transmit signals with other cells. In this study, we used the presynapse signature (PSS) as a representative signature (PSS is highly correlated with neuronal communication in TCGA, data not shown) to study the neural signal in TME. We started by checking correlations between PSS and postsynaptic signature across various tumors and validated the signature relevance. We then looked at the genes and pathways correlated to the PSS in TCGA data. Given that the percentage of neurons in tumors is relatively rare in most tumor types, we also looked at other stromal cells which can utilize or “mimic” neural signals. Aside from tumor cells and stromal cells, immune cells are a major component in the TME. We looked for possible effects of neural signals on different types of immune cells and confirmed our findings in in vitro dorsal root ganglion (DRG, a peripheral sensory neuron)/macrophage coculture experiments. Finally, we surveyed all tumor types included in TCGA to identify the types in which neural signals may affect patient outcomes.

MATERIALS AND METHODS

Data used for analysis

TCGA data

TCGA RNA‐seq data and clinical (phenotype) data were downloaded from the University of California Santa Cruz website (https://xenabrowser.net/datapages/). RNA‐seq data were scaled by 75% percentile expression level. Samples with scale factor <4 (low quality) or sample type code >= 10 (non‐tumor tissues) were excluded from the analysis. Scaled RNA‐seq data were log2 transformed (log2(RNAseq + 0.001) + log2(75), where the constant log2(75) is to make the TCGA data at the similar values of other RNAseq log2(TPM) values). This resulted in 9,805 tumor samples total. For the PAAD (pancreatic adenocarcinoma), we separated another subtype with high neuroendocrine signature, and designated them as PNET (14 samples).

Generation of brain enriched presynaptic gene set

Presynaptic gene set was obtained from Syngo portal (https://syngoportal.org/ontology.html). SynGo ontology is an evidence‐based, expert‐curated knowledge base for the synapse. There were originally 482 unique presynaptic genes downloaded from the portal. The mean Reads Per Kilobase of transcript per Million mapped reads (RPKM) of each gene in the brain tissue as well as the other tissues was found in GTEx v7 brain tissue database, also included in the supplementary data from the SynGo publication. The fold‐change enrichment in brain compared to other tissue was calculated, and an arbitrary enrichment of 1.2‐folds (20% more than the rest of the tissues) was used as a filter to select brain enriched presynaptic gene set, which led to a list of 256 genes. The expression levels of the 256 genes in the neuropil–a dense network of interwoven nerve fibers and their branches and synapses, together with glial filaments–were sorted based on the neuropil mRNA RNA‐seq result, the top 161 genes were kept for future analysis.

Deconvolution of TCGA data

Two algorithms were run on TCGA data to deconvolute cell types from the bulk RNA‐seq data. xCell (https://xcell.ucsf.edu/) was used to deconvolute 64 cell types, and we mainly used the xCell results for stroma cells. Cibersort (https://cibersort.stanford.edu/) was used to deconvolute immune cells. Immune cells and their components were used for the neural/immune interaction analyses. ESTIMATE (https://bioinformatics.mdanderson.org/estimate/rpackage.html) was used to estimate total immune, stromal, and tumor content.

Correlation calculation

R function “cor” and “cor.test” from “stats” package were used to compute correlations in the analyses.

General linear model for neural signals and neural immune interactions

The glm function from R stats package was used to fit presynapse score (PSS) with stromal or immune cells:

A constant 0.01 was added to estimated cell fractions in log2 transformation to avoid log2(0) problem.

The R function summary (fit) was used to extract the t‐ and p‐values from the fitting. R function fitted.values(fit) was used to compute correlations with PSS, and estimate the percentage of variances explained by the model.

DRG neuron and macrophage coculture assay

All animal studies were conducted in accordance with the American Veterinary Medical Association (AVMA) Guidelines for the Euthanasia of Animals, Guide for the Care and Use of Laboratory Animals, and the City of Cambridge. The Institutional Animal Care and Use Committee at Cygnal Therapeutics reviewed and approved all guidelines and procedures for rodent euthanasia and dissection of dorsal root ganglia. Policies are listed under standard operating procedure SOP513 and tissue collection protocol CYG2019‐003. Mice were sacrificed by carbon dioxide euthanasia.

C57/BL6 mouse DRG neurons were dissected and isolated following published protocols by Sleigh et al. In brief: mouse DRGs were extracted from adult mice following IACUC procedure. DRGs were dissociated by a mixture of collagenase and dispase for 30 min. Dissociated neurons were washed and plated on 50 µg/ml PDL (Sigma‐Aldrich, P7405‐5MG) and Matrigel (Corning Matrigel Basement Membrane Matrix, LDEV‐free, 10 ml. Product Number 354234, 1:30 diluted in PBS) coated wells in Neurobasal Medium‐A (Thermo Fisher Scientific, 10888022) containing B‐27 supplement (ThermoFisher, A3482801), and 25 ng/ml mouse NGF (Sigma, n6009). Final concentration of 5 µM anti‐mitotic reagent Ara‐C (Cytarabine, Sigma, C1768‐100 mg) was added to the culture medium to inhibit the proliferation of non‐neuronal cells such as fibroblasts. DRGs were cultured for 8 days before coculture was setup. DRG culture media was changed every 3 days with Neurobasal‐based media supplemented with B‐27, same concentrations of NGF and anti‐mitotic reagent Ara‐C.

Seven days before coculture, bone marrow cells were isolated from mouse leg bones. Cells were plated and cultured in DMEM supplemented with 10% FBS, 50 μM β‐mercaptoethanol (ThermoFisher, 21985023), and 20 ng/ml recombinant murine M‐CSF (Peprotech, 315‐02). Twenty‐four hours before coculture, M2 polarization was initiated by adding fresh medium containing 10 ng/ml mouse M‐CSF, 20 ng/ml mouse IL4 (Peprotech, 214‐14), and 10 ng/ml IL‐13 (Peprotech, 210‐13). On the day of coculture, DRG media was aspirated. M2 macrophages were lifted and added to DRG with or without 50 ng/ml LPS (Invivogen, tlrl‐eblps) in DMEM supplemented with 10% FBS and 50 μM β‐mercaptoethanol. Twenty‐four hours after coculture, cells were lifted from the plate by Accutase (Innovative Cell Technologies, AT104) and stained with corresponding biomarkers and analyzed by flow cytometry. For RNAseq, cells were lifted from the plate, total RNA was extracted for transcripts analysis.

RESULTS

PSS in tumors

To generate a PSS list representing the neural signaling in the tumor, we first obtained a raw list of presynaptic genes in an evidence‐based, expert‐curated database: SynGo (https://www.syngoportal.org/). We excluded genes that are primarily expressed in the cell soma because the vast majority of peripheral tumors do not contain the cell bodies of neurons. Neurons have complicated morphologies which allow the nerve processes, such as dendrites and axons, to extend far away from the nuclei. To accommodate these specialized structures, mechanisms have been developed to allow local protein translations at these nerve fibers, also called neuropils. Various labs have carried out neuropil transcriptomics to understand such mechanisms. Here, we resorted to neuropil (neural process) transcriptomics to filter the presynaptic genes that have transcripts in the distal nerve terminals. , That led to 161 PSS genes (see Section 2, Table S1). We computed a PSS score for each TCGA sample by averaging the log2 (expression level) of each signature gene. Figure 1A shows the boxplot of PSS scores across all 33 tumor types from TCGA (in fact, it is 34 tumor types since we split PAAD (pancreatic adenocarcinoma) into PNET (pancreatic neuroendocrine tumor) and PAAD, see Section 2). The scores are highest in tumors of the nervous system (GBM, glioblastoma multiforme, LGG, brain lower grade glioma, and PCPG, pheochromocytoma, and paraganglioma, see https://gdc.cancer.gov/resources‐tcga‐users/tcga‐code‐tables/tcga‐study‐abbreviations for tumor type abbreviations), but low in most of the solid tumors, demonstrating that the PSS signature contains genes that are specific to nervous system.

FIGURE 1

Presynapse signature (PSS) is highly expressed in neuron tumors and highly correlated with postsynaptic signature in TCGA. (A) Boxplot of PSS score in various TCGA tumor types. (B) Correlation of presynapse and postsynapse signature across all non‐neuronal TCGA tumor samples. (C–E) Correlation of presynapse and postsynapse signatures in TCGA neuronal tumors. GBM, Glioblastoma multiforme; LGG, Brain Lower Grade Glioma; PCPG, Pheochromocytoma and Paraganglioma

We then looked at the relationship between presynapse and postsynapse (also from SynGo database) signatures in various tumors in TCGA. Overlapping genes that could be either presynaptic or postsynaptic were removed from the postsynaptic signature for this analysis. These two signatures are highly correlated (correlation = 0.94, p = 0) in all non‐neuronal TCGA tumors (Figure 1B). We also checked the relationship in each of neuronal tumors individually (GBM: Figure 1C, LGG: 1D, and PCPG: Figure 1E), and they all show very tight correlations, in agreement with closely related biological processes. As a result, we believe the PSS and postsynaptic signature help corroborate each other.

Pathway enrichment analysis of PSS‐correlated genes in TCGA

To get an overview of the possible roles of neurons or neuronal signaling in the TME, we first looked at the PSS‐correlated genes and pathways in TCGA data. To avoid the dominant and confounding influence of neuronal tumors (GBM, LGG, and PCPG), these three tumor types were excluded from the analysis. We generated a PSS score for each tumor sample by averaging the log2 transformed expression levels and correlated them with each gene expression level across tumor samples (in 8934 non‐neuronal tumor samples). We found 2826 genes correlated with PSS with a correlation >0.3, and most highly enriched pathways were “neuron system,” “axon guidance,” “neuroactive ligand‐receptor interaction,” “collagen biosynthesis and modifying enzymes,” “potassium channels,” etc. (Figure 2A), and most of the ion channel genes (calcium, potassium, and sodium) are correlated with PSS (Figure 2B).

FIGURE 2

Genes and pathways correlated with presynapse signature (PSS) in TCGA. (A) pathways enriched in genes correlated with presynapse signature in all TCGA non‐neuronal tumors (all tumors except GBM, LGG, and PCPG). (B) Boxplot of correlation between individual genes and PSS in all TCGA non‐neuronal tumors. Most of ion channel genes (calcium, potassium, and sodium) are correlated with PSS. (C) Pathways enriched in correlated genes with PSS in STAD (stomach adenocarcinoma). (D) Pathways enriched in negatively correlated genes with PSS in STAD. (E) Pathways enriched in correlated genes with PSS in PAAD (pancreatic adenocarcinoma). (F) Pathways enriched in negatively correlated genes with PSS in PAAD. All results generated using Enrichr tool (https://maayanlab.cloud/Enrichr/)

A few tumor types had relatively higher average PSS scores than the other tumor types. For example, in STAD (414 tumors), 2110 genes had a correlation >0.5, and 459 genes had a correlation <−0.4. Pathway analysis showed that the neuron‐related pathways were highly enriched for the positively correlated genes (Figure 2C), whereas the proliferation and cell cycle‐related pathways were highly enriched in negatively correlated genes (Figure 2D). Similar pathway enrichment (Figure 2E,F) was observed in PAAD (165 tumor samples excluding PNET samples). These results further validate PSS as a reliable surrogate for neural signal in tumors despite lacking of direct representatives.

Stromal cells contributed to neuronal signal

To find out the origin of the PSS signals in TCGA, an xCell algorithm was used to deconvolute cell types from TCGA RNA‐seq data, and indeed the neuron was one of the cell types identified. The fraction of neurons estimated by xCell was very low in bulk tumors (Figure 3A), as neural cell bodies are typically not part of a tumor biopsy. Given the strong PSS signal, we thus asked the question if other cell types could contribute to neural signals, or “mimic” neurons.

FIGURE 3

Presynapse signature versus stroma cells. (A) Boxplot of deconvoluted fraction of neurons in each tumor type (by xCell). (B) Scatter plot of presynapse signature versus xCell deconvoluted stromal cells for six top ranked cell types: B1: neurons; B2: chondrocytes; B3: astrocytes; B4: fibroblasts; B5: HSC (hematopoietic stem cells); B6: mesangial cells. t‐ and p‐values are estimated by glm fit. (C) Heatmap of correlations between PSS and cell types in each of the TCGA tumors. Cell types are ranked by average correlation with PSS across all tumor types. (D) Modeling of presyanpse signature score by top six stroma cells plus tumor type in all TCGA non‐neuronal tumors. The model can explain 69% of the PSS variation. (E) Modeling of presyanpse signature score by top five stroma cells (excluding neuron) plus tumor type in all TCGA non‐neuronal tumors. The model can explain 58% of the PSS variation

To address this question, we deployed a two‐component general linear model for each cell type in the xCell, with one component being the estimated cell fraction (in log2 scale), and another component being the tumor type (see Materials and Methods). The t‐value from the general linear fit was used to rank the cell types. The top six cell types that could contribute to neural signals included neurons (t = 51), chondrocytes (t = 43), astrocytes (t = 43), fibroblasts (t = 40), HSC (hematopoietic stem cells) (t = 37), and mesangial cells (t = 30) (Figure 3B). The fact that PSS has the most significant association with deconvoluted neurons again confirms the relevance of this signature. Figure 3C shows the hierarchical clustering of correlation between PSS and all cell types in 34 tumor types, which shows that the PSS is highly correlated with the top six cell types in the vast majority of tumor types (top seven cell types are: neurons, fibroblasts, chondrocytes, astrocytes, HSC, total stroma score, and mesangial cells). A multi‐component general linear model was also used to assess the fraction of variations in the PSS score that could be explained by simple linear combinations of these stromal cells (Methods and Materials). The top six cell types (including neuron) can account for 69% of PSS variations across all non‐neuronal tumors (correlation 0.83, Figure 3D, and for individual tumor types, Figure S1). Excluding the neuron component, the remaining five cell types still accounted for 58% of PSS variation (correlation 0.76, Figure 3E). These results suggested that apart from neurons, stroma cells could be the other source of neural signals. Interestingly, tumor cells were negatively correlated with PSS (Figure S2A), which might simply be due to the fact that the fraction of tumor cells and the fraction of stroma cells are negatively correlated (Figure S2B), and that the tumor is unlikely a major source of neuronal signals. Correlations of several other fractions of cell types can be found in Figure S2C,D.

Neuronal signals and immune cells

In the bulk RNA sequencing, immune cells are a major component across almost all tumor types, and have a great impact on tumor cells. Next, we looked at the possible effect of neural signals on immune cells. To this end, we utilized the online Cibersort portal to deconvolute the various immune cells from TCGA RNA‐seq data, as the deconvolution algorithm focused on immune cells and provided detailed immune subtypes and polarization status. Using a similar general linear model as above, we calculated the t‐values of PSS score versus deconvoluted immune subtypes in each of the non‐neuronal tumor samples. The t‐value distributions of paired immune subtypes, including B cells, T cells, macrophages, and mast cells, are shown in scatter plots with colored tumor types in Figure 4. For all these pairs, the PSS was more associated with suppressive/less active subtypes than pro‐inflammatory/active subtypes (for example, more correlated to macrophages M2 than M1). This association was observed in a majority of tumor types (Figure 4E).

FIGURE 4

Presynapse signature versus immune cells (deconvoluted by Cibersort). (A) Scatter plot of presynapse signature versus B cells (left: naïve, right: memory), t‐ and p‐values are estimated by glm fit in all TCGA non‐neuronal tumors. (B) Scatter plot of presynapse signature versus CD4 memory T cells (left: resting, right: activated). (C) Scatter plot of presynapse signature versus macrophages (left: M2, right: M1). (D) Scatter plot of presynapse signature versus mast cells (left: resting, right: activated). (E) Heatmap of correlation between presynapse signature and pair of immune cells in each TCGA tumor types. (F) Scatter plot of xCell estimated neuron fraction versus total infiltrating immune cell fraction. Colors in scatter plots indicate different tumor types

We subsequently assessed whether the neural signals are correlated with overall immune cell fractions in solid tumors in general. By investigating the correlation of neural cell fraction (by xCell) with total immune cell fraction (by ESTIMATE) in each tumor sample, we found a negative correlation between these two components, suggesting that tumors with more neural fractions tend to have less immune cell fractions, resembling the so called “cold tumors” (Figure 4F).

A correlation only suggests a possible relationship; to investigate whether neuronal signals indeed drive a suppressive immune phenotype, we conducted in vitro validation experiments.

DRG neurons lock macrophages in an immunosuppressive phenotype

Macrophages represented one of the major immune populations in most solid tumors. The diversity and plasticity of macrophages can be shaped by the different tissue local microenvironments and in response to various physiological and pathological stimuli. , In response to varying types of activation, macrophages can be classified into classically activated macrophages (“M1”) and alternatively activated macrophages (“M2”). “M1” macrophages are immunostimulatory and are characterized by the production of high levels of pro‐inflammatory cytokines, while “M2” macrophages are immunosuppressive in nature and accompanied by the secretion of anti‐inflammatory cytokines. Macrophages usually display an anti‐inflammatory M2‐like status in tumors. ,

To decipher the crosstalk between neurons and anti‐inflammatory tumor‐associated macrophages, we developed an in vitro coculture system utilizing primary murine DRGs (peripheral sensory neurons) and primary “M2”‐like macrophages differentiated from mouse bone marrow.

To investigate the role of sensory neurons during “M2” into “M1” repolarization, “M2” macrophages were re‐polarized into “M1”‐like phenotype by lipopolysaccharide (LPS) stimulation during the coculture experiment. “M2” into “M1” repolarization of macrophages was evaluated by measuring the expression of classical “M2” marker Arginase1 and “M1” marker iNOS via flow cytometry (Figure 5A). In monoculture of M2 macrophages, repolarization with LPS for 24 h led to a drastic increase of the M1 macrophage marker iNOS protein as well as mRNA, while it led to only a small increase of M2 macrophage marker ARG1 protein or mRNA (Figure 5B,C).

FIGURE 5

DRG facilitates M2 polarization of macrophages. (A) Schematic of DRG‐macrophage coculture and polarization. (B) Representative FACS plots of iNOS and ARG1 expression in monoculture and coculture macrophages. (C) ARG1 and iNOS protein level in macrophages measured by FACS. (D) Arg1 and Nos2 (encoding iNOS) gene expression levels by RNAseq

Compared to monoculture macrophages, the iNOS protein level was reduced by over fivefold under LPS repolarization in macrophages cocultured with DRGs, suggesting a role of DRG in inhibiting M2 to M1 repolarization. In coculture M2 macrophages, the protein level of ARG1 was increased by twofold without repolarization compared with monoculture. After 24 h of LPS repolarization, the ARG1 protein level was higher by threefold in coculture compared with monoculture (Figure 5B,C). Besides the protein expression levels, we also observed the mRNA levels of Arg1 and Nos2 (encoding iNOS) were changed accordingly (Figure 5D).

In accordance with the in vitro results, a correlation analysis of 161 PSS genes with either M2 or M1 macrophages showed the majority of PSS genes were positively correlated with M2 while negatively correlated with M1 (Figure S3). In conclusion, our work suggests that neuron signals promote the immunosuppressive phenotype of macrophages, which confirms the observation of PSS association with M2 but less of M1 as shown in Figure 4C.

Neuronal signals versus patient outcome

Observing the PSS correlations with immune cells, we analyzed how these correlations could be translated to patient outcome and therefore performed log‐rank test over TCGA tumors. For each tumor type, patients were evenly divided into three groups based on the PSS score (top, middle, and bottom one‐third), and log‐rank tests were performed over three clinical endpoints: overall survival (OS), disease‐specific survival (DSS), and progression‐free interval (PFI). Among 33 tumor types (PNET excluded due to low sample number), we found that high PSS scores were associated with good outcomes (OS, DSS, PFI, Figure S4A–C) in LGG, but poor outcomes in stomach adenocarcinoma (STAD; OS, DSS, PFI, Figure 6A–C), bladder urothelial carcinoma (BLCA; OS, DSS, PFI, Figure 6D–F), uterine corpus endometrial carcinoma (UCEC; OS, DSS, Figure S4D–F), and colon adenocarcinoma (COAD; DSS, Figure S4G–I). Worth noting that when a correlation between PSS and outcome was observed, a high PSS score was generally correlated with poor prognosis in non‐neuronal tumors. Even though PAAD had higher PSS scores among non‐neuronal tumors, however, the PSS score was not prognostic in this tumor type (Figure S4 J–L).

FIGURE 6

Presynapse signature versus patient outcomes. Kaplan–Meier curves of STAD (stomach adenocarcinoma) and BLCA (bladder urothelial carcinoma) patients grouped by PSS scores, red, green, and blue represent low, mid, and top 1/3rd patients. P‐value is estimated by log‐rank test. (A) STAD overall survival (OS), (B) STAD disease‐specific survival (DSS), (C) STAD progression‐free interval (PFI). (D) BLCA overall survival (OS), (E) BLCA disease‐specific survival (DSS), (F) BLCA progression‐free interval (PFI)

DISCUSSION

To carry out meaningful gene expression correlation analyses, the selection of the PSS gene list is critical. We used peer‐reviewed publications for well‐accepted genes that are enriched in presynapses. SynGo PSS gene list included some genes that were not exclusively presynaptically localized, but they had been indicated in the literature for modulating presynaptic functions. We further narrowed down the PSS gene list by focusing on genes with enriched transcripts detectable in the neuropil analysis, since most neuronal fragments in tumors come from the neurites. It is of note that the neurite‐specific transcript publications were mostly based on rodent brain tissues, and thus were likely an approximation of what is happening in human tumor‐innervating nerves. If the tumor‐innervating nerves are very different from those of central nervous system (CNS), we would miss these genes. However, neural signaling pathways are relatively conserved across species. When we looked at the correlation of the PSS with neuronal tumor types as well as postsynaptic signatures, they had extremely high correlations. This suggests like neurons, the TME components may both “emitting” and “receiving” neural signals. On the other hand, this also suggests that the resolution of the current analysis was limited, and we could not differentiate different types of neural signals. Future experiments that detect mRNAs in the DRG or nodose ganglia neurites will provide valuable resources for a refined gene signature for tumor‐innervating nerves. Using the same methodology, other gene lists such as glial/Schwann cell markers could also be used to further validate the correlations between nerves and glia, as well as glia with the immune cells.

The analyses in this study were mostly based on correlations between gene signatures and deconvoluted cell types in a large reference dataset. While correlations may suggest possible relationships, thus providing hypotheses to be further validated or tested, they can occasionally be influenced by hidden variables (confounding factors). For convincing causal relations, the hypotheses need to be validated by carefully designed in vitro or in vivo perturbation experiments. When the PSS gene set was highly correlated with a certain cell type (deconvoluted cell types by xCell), it could either suggests that these cell types coexist with innervation, or that they could resemble nerve signaling, such as neuroendocrine cells and adrenal gland. Neural signal was negatively correlated with proliferation signature in many tumor types. It could either be that neural signals were inhibiting tumor cell proliferation or equally possible that tumor content was inversely proportional to the number of stromal cells, thus creating an artificial negative correlation. In fact, when we cocultured various tumor cell lines with various neuronal cells in vitro, we saw both promoting and inhibition of tumor proliferation (for example, DRG cells promote proliferation in small‐cell lung cancer and prostate cancer cell lines, but inhibit proliferation in pancreatic cell lines, data not shown).

When looking at the relationship between PSS and immune subtypes, correlation analysis showed that the PSS was more associated with suppressive/less active subtypes than pro‐inflammatory/active subtypes. Furthermore, tumors with higher neuron fractions tend to have less immune cell fractions.

How could the neural signaling modulate immune cells in TME? It was found that the nervous system forms close associations with immune cells such as T cells, macrophages, neutrophils, NK cells, mast cells, and dendritic cells. , , Neuropeptides, transmitters, cytokines, chemokines, and purines could modulate the immune cells from various perspectives, such as propagation, induction, recruitment, migration, maturation, T‐cell priming, and macrophage phagocytosis. Neurons have been shown to play a role in multiple diseases such as arthritis, asthma, diabetes (pancreatitis), inflammatory bowel disease, psoriasis, contact dermatitis as well as inflammatory diseases based on neural manipulation/modulation experiments in vitro as well as in vivo. , Pro‐ and anti‐inflammatory effects from the nervous system have been observed, suggesting a dynamic interaction between the nervous system and the immune system. , , Though these mechanisms have not been studied in the context of tumor, they might share similar mechanisms like PSS association with immunosuppressive states in various tumor types.

Whether other neural signal contributing cell types like fibroblasts can interact with macrophages similarly as bona fide nerves requires future dedicated experiments. We reason it is likely based on (a) the TCGA PSS and M2 macrophage correlation was observed across many solid tumor types, even though some of the tumor types had very little deconvoluted neurons, and (b) more than 60% of the presynapse signature genes are expressed (TPM > 1) in TME myofibroblasts (for example, LRRC15+ fibroblasts in PADC, NSCLC, Figure S5), where the fibroblast expression level was derived from single cell RNA‐seq data by Buechler et al. This suggests fibroblasts have needed molecules to perform the neural function.

Given the observation that neural signal is associated with immune‐suppressive component, we looked the neuropeptides and neuropeptide receptors in terms their correlation with immune cells in TME, and found the neuropeptide receptors are significantly more correlated to M2 than M1 macrophages (t‐test p = 2.3 × 10−11, Figure S6). Note that in the model system of mouse DRG and M2 coculture, significant number of neuropeptides and receptors are indeed expressed with TPM > 1. To identify the particular group of neural signals affecting immune polarization would require more investigations.

In summary, our analysis based on TCGA dataset investigated a neural gene expression signature in tumor settings and observed that (1) multiple stromal cell types may utilize neural signals as communication tools, (2) neural signals were associated with immune suppressive or inactive TME, which was further confirmed by in vitro DRG‐macrophage coculture experiments, and (3) the PSS was predictive of patient outcome in multiple, non‐neuronal tumor types. To the best of our knowledge, this is the first time that the neural signaling has been revealed to correlate with suppressive immune subtypes in TME across many tumor types.

AUTHOR CONTRIBUTIONS

H. Dai and S. Lou performed the analysis and drafted the paper; Y. Zhang, C. J. Shu, and S. Carden contributed to the in vitro DRG/M2 experiment; M. Thanawala compiled presynaptic/postsynaptic signatures; KC. Huang and L. Ji contributed to TCGA data analysis; T. Liao and M. Abbassi performed mouse DRG extraction; A. Lantermann, M. Sadaghiani, D. Blom, J. Wagner, and P. Huang contributed to analysis discussions, manuscript editing.

Supplementary Material

Click here for additional data file.