Neuron-specific spinal cord translatomes reveal a neuropeptide code for mouse dorsal horn excitatory neurons

The spinal dorsal horn harbors a sophisticated and heterogeneous network of excitatory and inhibitory neurons that process peripheral signals encoding different sensory modalities. Although it has long been recognized that this network is crucial both for the separation and the integration of sensory signals of different modalities, a systematic unbiased approach to the use of specific neuromodulatory systems is still missing. Here, we have used the translating ribosome affinity purification (TRAP) technique to map the translatomes of excitatory glutamatergic (vGluT2+) and inhibitory GABA and/or glycinergic (vGAT+ or Gad67+) neurons of the mouse spinal cord. Our analyses demonstrate that inhibitory and excitatory neurons are not only set apart, as expected, by the expression of genes related to the production, release or re-uptake of their principal neurotransmitters and by genes encoding for transcription factors, but also by a differential engagement of neuromodulator, especially neuropeptide, signaling pathways. Subsequent multiplex in situ hybridization revealed eleven neuropeptide genes that are strongly enriched in excitatory dorsal horn neurons and display largely non-overlapping expression patterns closely adhering to the laminar and presumably also functional organization of the spinal cord grey matter.

The ability to sense and discriminate different noxious and innocuous somatosensory stimuli is essential for all higher animals and humans in order to react adequately to external stimuli and internal conditions 1,2 . The spinal dorsal horn, i.e., the sensory part of the spinal cord, constitutes a key element in this process. It receives somatosensory signals from peripheral neurons and processes these signals together with other inputs descending from supraspinal sites in a complex network of inhibitory and excitatory interneurons before relaying them via projection neurons to supraspinal centers 3 . Projection neurons make up less than 10% of all dorsal horn neurons, while more than 90% of the neuronal population are interneurons of which, between 60 and 70% are excitatory glutamatergic neurons, and the rest is inhibitory (GABA and/or glycinergic) 3,4 . Besides these principal neurotransmitters, many neurons of the spinal dorsal horn use additional signaling molecules, especially neuropeptides, for neuronal communication, which together with the structural connectivity set the basis for modality specific sensory processing.
The spinal cord is organized in a laminar fashion, which has initially been proposed on the basis of differences in cell density and morphology 5,6 but has later been shown to also reflect a functional organization. This becomes obvious especially from the lamina specific innervation pattern by the different types of peripheral sensory neurons: unmyelinated C fibers, which mainly carry noxious and thermal information, terminate in the superficial dorsal horn (laminae I-II), while thickly myelinated Aβ fibers, which convey innocuous signals including touch and proprioceptive information, terminate in the deep dorsal horn (laminae III-V) 3,7,8 . A laminar organization of neuronal function is also supported by gene expression patterns that follow laminar patterns 9-13 and a recent publication indicating a topographic map for spinal sensorimotor reflexes 14 . Furthermore, optogenetic and chemogenetic experiments support a modality-specific processing by distinct genetically defined neuron populations 8,[15][16][17] . However, the basis of this modality specific processing beyond structural connectivity is only incompletely understood.
In the present study, we performed a genome-wide unbiased search for genes involved in neuronal communication in excitatory or inhibitory spinal cord neurons. To this end, we employed the translating ribosome affinity purification (TRAP) technology 18 to characterize the translatomes (i.e. the whole range of mRNAs that are not only transcribed but also translated). We generated vGluT2::bacTRAP (Slc17a6), vGAT::bacTRAP (Slc32a1) and Gad67::bacTRAP (Gad1) mice, which express the eGFP-tagged ribosomal subunit L10a (RPL10a) under the control of the respective gene regulatory elements. Gene ontology and pathway analyses revealed that, in addition to transcription factors and genes related to the neurotransmitter phenotype of these neurons, genes encoding elements of neuropeptide signaling constitute a functionally defined class of genes that is highly enriched in excitatory spinal neurons. About 25% of the most strongly enriched genes were neuropeptides. Subsequent multiplex in situ hybridization for neuropeptide encoding genes revealed largely non-overlapping expression patterns that follow the laminar organization of the spinal cord thus suggesting a role of neuropeptide signaling in segregation of sensory modalities.
From the spinal cord of these mice, polysome-bound mRNA can be isolated by co-immunoprecipitation of mRNA bound to the eGFP-L10a tagged ribosomes. We confirmed the correct expression of the eGFP-L10a transgene in the lumbar spinal cord of the different mouse strains by immunohistochemistry. NeuN was used as a neuronal marker, Pax2 as a ubiquitous marker for spinal inhibitory 19,20 and Lmx1b as a marker for the majority of dorsal spinal excitatory neurons ( Fig. 1D-I). In the vGluT2::bacTRAP line, the eGFP-L10a expression is restricted to Pax2-negative neurons (98.8 ± 0.5% Pax2 − NeuN + ; Fig. 1D,G,J), indicating its exclusive expression in excitatory spinal neurons. In the two mouse lines targeting inhibitory neurons, vGAT::bacTRAP and Gad67::bacTRAP, the expression was confined to Pax2-positive neurons (98.3 ± 0.8%, 97.4 ± 1.6%, respectively) and virtually absent from excitatory, Lmx1b-positive neurons of the dorsal horn (0.25 ± 0.22%, 0.19 ± 0.34%, for vGAT::bacTRAP and Gad67::bacTRAP mice, respectively; Fig. 1E,F,H,I,K,L). Next, we determined the fraction of vGluT2, vGat and Gad67 mRNA expressing neurons that co-express the eGFP-L10a (TRAP) transgene in the respective TRAP driver line. To this end, we performed multiplex in situ hybridizations on lumbar spinal cord sections of Gad67::bacTRAP, vGAT::bacTRAP and vGluT2::bacTRAP mice (Fig. S1). 96.5 ± 3% of Gad67 + neurons ( Fig. S1A-C), 93.8 ± 4% of vGAT neurons (Fig. S1D-F) and 95.4 ± 1.2% of vGluT2 neurons (Fig. S1G-I) expressed the TRAP transgene in the respective mouse line. To further validate the three bacTRAP lines we isolated the mRNA bound to the eGFP-L10a tagged ribosomes from each line and analyzed the enrichment of selected marker genes in the polysome pull down fraction versus the input using qPCR. In all cases, we found an enrichment of the transcript encoding the transgene (eGFP) (Fig. 1M-O). In addition, when comparing the pull down fraction versus the input fraction isolated from the vGluT2::bacTRAP line we found an enrichment of transcripts encoding excitatory markers (Lmx1b and vGluT2) and a depletion of the transcripts encoding for inhibitory markers (vGAT , GlyT2 and Gad1) (Fig. 1O). Conversely, in the two lines expressing the TRAP construct in inhibitory neurons we found an enrichment of vGAT, GlyT2 and Gad1 and a depletion of Lmx1b and vGluT2 (Fig. 1N,M).   www.nature.com/scientificreports/ General differences in gene expression in excitatory and inhibitory spinal neurons. Next, we used the bacTRAP mouse lines to identify key features and subtype specific markers of the two main subpopulations of spinal interneurons (i.e. excitatory or inhibitory dorsal horn neurons). To this end, we isolated and sequenced the cell-type-specific polysomal mRNA from the lumbar spinal cord of three animals from each mouse line. Out of 22,880 detectable genes, between 13,515 and 13,774 were found expressed in each mouse line ( Fig. 2A, Table S1), defined by a normalized count (number of reads per gene, normalized to the total number of reads and the gene length) of > 10. Table S1 thus presents the first source of a browsable genome wide atlas of translated mRNAs in excitatory and inhibitory spinal neurons. In order to identify those genes that show a regionally non-overlapping expression pattern, we focused our analyses on genes whose regional expression was accessible by in situ hybridization. To this end, we assumed that genes with normalized counts > 100 should be reliably detectable with in situ hybridization. This threshold is three times lower than the expression of Grpr, a gene that we previously found expressed at low levels in a small subset of dorsal horn neurons. In all three mouse lines, approximately 11,200 genes were expressed at levels exceeding this threshold ( Fig. 2A). Differential gene expression analyses (DGEAs) revealed 800 genes that were differentially expressed (with a false discovery rate (FDR) ≤ 0.05) in vGluT2::bacTRAP versus vGAT::bacTRAP mice, and 1,790 genes in vGluT2::bacTRAP versus Gad67::bacTRAP mice (Fig. 2B). The approximately two-fold lower number of significantly enriched genes in the comparison of vGluT2::bacTRAP versus vGAT::bacTRAP mice probably results from the lower expression of the eGFP-L10a transgene in the vGAT::bacTRAP line and hence a lower mRNA yield and higher variability in read numbers. When we restricted the analysis to genes with an at least two-fold enrichment, comparison of the translatomes of vGluT2 positive and vGAT positive neurons revealed 438 differentially expressed genes. Of those, 256 genes were enriched in excitatory neurons, and 182 in inhibitory neurons. When the Gad67-TRAP line was used instead of the vGAT-TRAP line, a higher number of differentially expressed genes were detected (624, of which 300 and 324 were enriched in excitatory and inhibitory neurons, respectively) (Fig. 2B). Between 70 and 80% of the differentially expressed genes identified with the vGAT-TRAP line were also found in the comparison with the Gad67-TRAP line (Fig. 2C), supporting the robustness of our approach. Subsequent cluster analyses ( 21 (http://fgcz-shiny .uzh.ch/fgcz_heatm ap_app/)) over the nine different samples revealed, as expected, that the three expression analyses made for each mouse line clustered together. As expected, the next hierarchical segregation separated the inhibitory (vGAT and Gad67) from the excitatory (vGluT2) samples. Within the inhibitory cluster vGAT samples segregated from Gad67 samples (Fig. 2D,E). The cluster analysis therefore identified the vGAT + , Gad67 + and vGluT2 + samples as three populations with distinct translatome profiles.
To further validate our bacTRAP approach we compared our gene expression results to recently published single cell sequencing data. We chose the 50 genes that were most enriched in one of the three bacTRAP lines employed in our study to the single cell sequencing data of Haring et al. 9 (Tables S2-S5). We found that many of our hits were also present in the single cell data (60-70% of the 50 most enriched genes). Fourteen to 20 genes identified per comparison in our study were either not detected or only at a level that was too low to assign them to a subpopulation. Among those genes that were not detected by the single cell analysis were many transcription factors such as Vsx2, Lmx1b, Pax5, Barhl1 or Phox2a but also other genes (e.g. Ighg3). To verify the expression of some of these genes, we selected three for further investigation with in situ hybridization. We verified expression of Barhl1 and Phox2a, which we found enriched in excitatory cells and the immunoglobulin heavy chain gene Ighg3, which we surprisingly found enriched in inhibitory neurons (Fig. S2).
Finally, we were interested in the correlation of the genes highly enriched in our study with the molecular groups that had been defined by Haring et al. 9 . We therefore analyzed the distribution of the 50 most enriched genes among the 30 neuronal subpopulations identified by Haring et al. (Tables S2-S5). Among the genes found enriched in excitatory neurons in our study (Tables S2, S4), three genes (Nrn1, Cacna2d1 and Slc17a6) displayed expression in all 15 excitatory subtypes of Haring et al. Three other genes (Ucn3, Cpne4, Aldh1a2) have been attributed to only one excitatory subtype, while all other genes (44) were either attributed to more than one subtype or not reported in the single cell analysis. Among the genes enriched in inhibitory neurons (Tables S3,  S5), seven genes displayed pan-inhibitory expression, seven were assigned to only one subtype and the remaining genes were assigned to more than one subtype or not reported. All genes that we found to be enriched in either inhibitory or excitatory neurons were (if detected) also exclusively or preferentially ascribed to the same category by Haring et al. (http://linna rsson lab.org/dorsa lhorn ). This comparison suggests first that the BAC-TRAP approach has a lower detection limit than in the single cell RNA sequencing approach and, second, that most of the subpopulation specific genes detected with single cell sequencing are not only transcribed but also translated. type of genes are enriched in inhibitory or excitatory spinal neurons we employed the Ingenuity Pathway Analysis (IPA, QIAGEN) tool, which assigns a unique subcellular location ("location") and coarse function ("type(s)") to the enriched genes ( Fig. 2F,G). Twenty-four percent of the genes enriched in vGluT2 over vGAT encoded for cytoplasmic proteins, 24% for extracellular proteins, 22% for proteins in the plasma membrane, 20% were nuclear proteins and the remaining 11% were located elsewhere. Similar localizations were found when genes enriched in the inhibitory populations were analyzed (Fig. 2F). One hundred-thirty out of 256 genes enriched in vGluT2 over vGAT and 109 of 182 genes enriched in vGAT over vGluT2 could be assigned to a coarse function (Fig. 2G). The largest group of genes that distinguished inhibitory from excitatory neurons in the adult spinal cord were transcription regulators (28% and 39% of the genes enriched in vGAT and vGluT2, respectively). Two other large molecular groups present in the set of enriched genes were enzymes and G protein coupled receptors (Fig. 2G).
The Ingenuity "pathway analysis" tool confirmed some the results of the GO term analyses but also identified additional receptor signaling pathways overrepresented in the different neuron populations (Fig. 4A-D).
Glutamate and GABA receptor signaling as well as serotonin receptor signaling were overrepresented in inhibitory neurons (Fig. 4B,D). As expected, the enrichment of GABAergic signaling was mainly due to genes involved in the production (Gad1, Gad2) and transport (vGAT , Gat1) of GABA (Tables S1, S2-S5), while GABA receptors showed relatively similar expression levels in the different populations (Table S8). Overrepresentation of Glutamate receptor signaling was due to the aforementioned enrichment of kainate receptors as well as to an enrichment of AMPA (Gria1), NMDA (Grin2c), and metabotropic (Grm2, Grm3) glutamate receptor subunits. Overrepresentation of serotonin receptor signaling (Fig. 4B) was due to the enrichment of serotonin receptors (Htr1a, Htr1d, Htr3a) in inhibitory neurons. Serotonin receptor 3a (Htr3a), encoding ionotropic serotonin receptors, was highly enriched in inhibitory spinal interneurons (Table S13) and displayed an expression pattern that was restricted to a population in the deep dorsal horn (Fig. S3). Tables S6-S15 illustrate relative expression levels and enrichment of glutamate, GABA, glycine, acetylcholine, serotonin, adrenergic, dopaminergic and histamine receptors. In summary, several receptor-signaling pathways are enriched within the differentially expressed genes, indicating that the respective neurotransmitter signaling pathways may engage inhibitory and excitatory spinal neurons differentially.
A neuropeptide code for glutamatergic spinal dorsal horn neurons. Neuropeptide signaling related terms (e.g. neuropeptide hormone activity, neuropeptide receptor activity and neuropeptide receptor binding, 4.6-26.3-fold overrepresented as compared to all expressed genes, FDR between 0.00016 and 0.0488) were highly overrepresented among the molecular function terms. These signaling terms were predominantly enriched in excitatory neurons (Fig. 3A,C). A comparison between the translatomes of the vGluT2 and Gad67 neuron populations demonstrated that 12 of the 50 genes most strongly enriched in excitatory neurons (Tables S2-S5)   www.nature.com/scientificreports/ three out of the 50 most enriched genes, neuropeptide enrichment was far less common in inhibitory neurons than in excitatory populations. In order to more deeply investigate the expression patterns of the neuropeptideencoding genes enriched in excitatory neurons, we employed fluorescent in situ hybridization (RNAscope). Of the 13 identified genes, two genes, adenylate cyclase activating polypeptide (Adcyap1) and corticotropin-releasing hormone (Crh), could not be detected reliably with in situ hybridization experiments. For the remaining 11 genes, we performed multiplex RNAscope experiments, combining up to three different neuropeptides at a time. Figure 5 (Fig. 5A-O) displays five of these combinations, in which every neuropeptide is represented at least once. Quantification of the co-expression of the depicted neuropeptides (Fig. 5C,F,I,L,O) suggested that most are expressed in different neuronal subtypes, yet one combination (Tac2 and Nmu) displayed a high degree of co-expression (Fig. 5F). To systematically determine the overlap between the different neuropeptide-expressing populations, we quantified for each pair of neuropeptides the percentage of co-expressing neurons. As some neuropeptides displayed a wide range of expression levels, we performed separate analyses either including low expressors (threshold near to background levels) or taking only cells into account with an expression clearly above background (threshold > 20 dots per cell). We created two expression matrices summarizing the percent of co-expressing neurons either taking only medium to high expression levels into account (Fig. 6) or all neuropeptide expressing neurons irrespective of the expression level (Fig. S4). Restricting the analysis to the medium to high expressors led to more pronounced segregations, accentuating the absence or presence of co-expression. When focusing on medium to high-level expression, our analyses suggested that at least eight excitatory subpopulations are marked by either a single or a combination of several neuropeptides: The expression patterns of several neuropeptides described above suggest a spatially restricted expression limited to only one or two spinal laminae (Fig. 5). To address the spatial distribution of each neuropeptide in more detail we imaged spinal cords at high resolution at the L3 and L4 levels. High resolution tiles were stitched, and the resulting images where aligned to reference sections from the spinal cord 26 (Fig. S5). We then determined the number of cells expressing the respective neuropeptide in each lamina. Nine out of the eleven neuropeptides were predominantly expressed in only one or two spinal laminae (Fig. 7A). Npff neurons were concentrated in lamina I, Grp, Nmu, Tac1 and Tac2 neurons were present mainly in lamina II, whereas Nts, Trh and Ucn3 neurons predominated in lamina III. Cck neurons were mainly located in the deep dorsal horn (laminae III-V). In case of the remaining two neuropeptides, Nmb or Nms, expression patterns were weak and dispersed over most of the dorsal horn and were therefore not quantified. These data indicate that the different laminae of the spinal dorsal horn are patterned by the expression of unique combinations of only partially overlapping neuropeptide genes (Figs. 6, 7B).

Discussion
In this study we have characterized the mRNA translation profiles (translatomes) of inhibitory (GABA/glycinergic) and excitatory (glutamatergic) spinal cord neurons. With this data we provide the first searchable dataset of neuron specific translatomes for the spinal cord. Pathway and Gene Ontology analyses of the translatomes of these two populations revealed that subtypes of neurotransmitter receptors and of down-stream signalling cascades differ between excitatory and inhibitory neurons and that the heterogeneity within the large population of excitatory dorsal horn neurons and their lamina organization is impressively mirrored by the expression of 11 neuropeptides or neuropeptide precursors.

Correlation of transcriptome and translatome in spinal neurons. mRNA translation profiles
obtained by the TRAP approach provide information about cell type-specific gene expression but without the cellular resolution obtained with single cell transcriptome profiling. However, the fraction of polysome-bound mRNA should more accurately reflect protein expression than approaches based on total or especially nuclear RNA analyses 18 . This is also supported by our study. The great majority of differentially expressed genes, revealed in our translatome comparisons, have previously been identified as marker genes of neuron subpopulations in studies using single cell RNAseq. However, we also found genes enriched in inhibitory and excitatory neurons that were either not detected with single cell RNAseq or detected only at levels that were too low to reliably assign them to a subpopulation (e.g. Phox2a). Other genes may have been filtered out to reduce artefacts (e.g. cell stress responses) potentially introduced by the necessary dissociation of the tissue (e.g. Ighg3). Hence, the sensitivity of the TRAP approach appears to be at least equal to single cell RNA sequencing. At the same time, our comparison of the most enriched genes versus data from single cell RNA sequencing also suggests a relatively high congruence between the transcriptome and the translatome. A comprehensive comparison of our translatome data with a recently developed harmonized transcriptional atlas of the spinal cord 27 may provide additional insights into the degree to which transcribed mRNAs are also translated.

Differential use of neurotransmitter and neuromodulator signaling pathways. Spinal inhibi-
tory and excitatory neurons receive synaptic inputs from a wide variety of peripheral, spinal and supraspinal neurons, including glutamatergic, GABAergic, glycinergic, serotonergic, cholinergic, and adrenergic neurons. We found that both inhibitory and excitatory spinal neurons express receptors for all these neurotransmitters (Tables S6-15). However, in some cases, neurotransmitter signaling appears to occur via different receptor sub-  Kainate receptors display different modes of signaling and a metabotropic action of kainate receptors has been suggested to be involved in presynaptic regulation and inhibition of GABA release 28,29 . A preferential expression of kainate receptors on inhibitory spinal neurons might therefore be involved in a glutamate-mediated inhibition of GABA (and potentially glycine) release. Even more striking is the enrichment of several serotonin receptors on inhibitory neurons with the most pronounced example being Htr3a, one of the two ionotropic serotonin receptor subtypes. Expression of Htr3a is more than eightfold enriched in inhibitory neurons and in situ hybridization indicated expression restricted to the deep dorsal horn. It is therefore conceivable that serotonin released in the spinal cord from descending supraspinal neurons has a differential effect on excitatory and inhibitory neurons.

Lamina specific expression of neuropeptides suggests a neuropeptide code involved in modal-
ity specific sensory processing. Our data confirm and extend previous single cell studies that identified neuropeptides as markers of excitatory neuronal subtypes 9,10 . Our multiplex in situ hybridization experiments demonstrate a lamina-specific distribution of many of the enriched neuropeptides with rather limited co-expression between different neuropeptides. In total, we find that the 11 different neuropeptides analysed in this study span dorsal horn laminae I-V, thus potentially setting up a map for a neuropeptide sensory signalling code. This is in good agreement with recent literature, as a number of these genes have either already been suggested to define specific subsets of spinal neurons 11,[30][31][32] or have already been used as driver genes for functional manipulation of spinal neuron subtypes 14,[33][34][35][36][37][38][39] . These studies demonstrate that neuronal subpopulations marked by the expression of specific neuropeptides exert modality specific functions in sensory processing. Grp+ and Ucn3+ neurons are required for the relay of pruriceptive information 33,37,39 , Tac1+ neurons are required for pain associated coping behaviours, Nts+ neurons are suggested to contribute to temporal summation of nociceptive input and "wide-up" 34 and Cck expressing neurons are crucially involved in the development of mechanical allodynia 35,38 . More importantly, emerging evidence suggests that neuropeptides do not only label specific subtypes of spinal neurons but are themselves crucial for modality specific sensory processing. The release of the neuropeptides Cck and Nts has been associated with pronociceptive or antinociceptive effects [40][41][42] . Signal- www.nature.com/scientificreports/ ling via the gastrin releasing peptide (Grp), which is concentrated in lamina II, is critical for itch perception in mice 33,39,43,44 and required for successful signal transmission from second order Grp to third order Grp receptor (GRPR) positive neurons in lamina II of the spinal cord 45 . Another example of a neuropeptide modulating pruriceptive information at the spinal level is somatostatin (SST). We found SST to be enriched in spinal excitatory neurons but we also detected substantial expression in inhibitory neurons (Table S1). SST expression in inhibitory dorsal horn neurons has been previously detected by Proudlock et al. 46 and Haring et al. 9 (see http://linna rsson lab.org/dorsa lhorn /) but was not reported by others 47,48 . However, SST is primarily expressed in peripheral neurons and in excitatory lamina II dorsal horn neurons. After release of SST from peripheral sensory neurons or spinal dorsal horn neurons it inhibits inhibitory dynorphin (Dyn) neurons, which in turn inhibit Grp receptor expressing neurons 49 . SST release thereby facilitates pruriceptive signalling. In this context, it is interesting to note that while the majority of spinal SST expressing neurons are excitatory (glutamatergic), the effect of SST on Dyn neurons is inhibitory. These studies demonstrate that neuropeptides are able to add an additional layer of complexity to neuronal communication and may in certain cases dominate over the effect of fast-acting neurotransmitters such as glutamate, GABA and glycine. The lamina specific expression of many neuropeptides in dorsal horn neurons suggests that distinct groups of neuropeptide expressing cells are activated by a restricted subset of primary afferents. Hence, distinct neuropeptides might only be released in response to certain sensory stimuli.

Conclusions
The present study has utilized three newly generated lines of BAC TRAP transgenic mice that allow cell typespecific translatome analyses of excitatory glutamatergic and inhibitory GABAergic/glycinergic neurons. In the present study, these mice allowed us to identify several signalling pathways enriched in distinct neuron populations of the dorsal horn. Particularly interesting to us appeared an enrichment of neuropeptide signalling in different populations of excitatory neurons with several neuropeptides exhibiting a rather lamina specific distribution suggesting an important role in modality specific sensory processing. While such a specific involvement is already known for some of these neuropeptides, the function of others is still obscure. The mouse lines described in this report and the accompanying translatome database should hence be well suited to reveal hitherto unknown mechanisms of normal sensory processing and their changes in disease states.

Materials and methods
Animals. Experiments were performed on 6-12-week-old mice kept at a 12:12 h light/dark cycle that received food and water ad libitum and were kept in a 12 h light/dark cycle. All methods were carried out in accordance with relevant guidelines and regulations. Genotyping PCRs. The presence of the eGFP-L10a transgene in the three different mouse lines was determined using a specific PCR for the transgene or a generic PCR against eGFP with one of the following primer pairs:  For immunofluorescence staining the slides were washed for 5 min in PBS to remove the embedding medium, followed by blocking with 5% normal donkey serum (NDS, AbD Serotec, RRID:SCR_008898) in 0.1% Triton X-100-PBS for at least 30 min at room temperature (RT). Sections were incubated with primary antibodies (see resource table) in the blocking solution at 4 °C overnight, followed by three washes in PBS for 5 min and incubation with secondary antibodies (fluorophore-coupled donkey antibodies, Jackson ImmunoResearch) in blocking solution at RT for 30-60 min. Subsequently, sections were washed in PBS and incubated with 4′,6-diamidino-2-phenylindole (DAPI, 500 ng/ml in PBS) for 10 min and again washed in PBS. Finally, sections were covered with DAKO fluorescent mounting medium (Dako, RRID:SCR_013530) and coverslips.
In situ hybridization (ISH). For ISHs, 6 weeks old, naïve male C57BL/6J mice were used. For spinal cord preparation, the vertebral column containing the lumbar part of the spinal cord were dissected, immediately after decapitation of the mouse. A syringe filled with ice-cold, RNAse-free PBS was fitted tightly against the caudal opening of the spinal canal and pressure was applied to press the spinal cord out, which was timed to the lumbar part. After dissection, spinal cords were snap frozen in liquid nitrogen and stored at − 80 °C until embedding in NEG50 frozen section medium (Richard-Allen Scientific). Frozen blocks were again stored at − 80 °C until cutting into 16-25 μm thick sections and mounted as described for IHC.
For single-plex chromogenic ISH the manual RNAscope 2.5 HD BROWN Assay (ACD, Bio-Techne, catalog no. 322300) was used. The manufacturer's pretreatment protocol for fresh frozen tissue (document no. 320536-TN, rev. date 11112016) and detection protocol (document no. 322310-USM, rev. date 11052015) were followed. Hybridization with Amp 5 was increased from 30 to 60 min for all probes in order to increase the signal. Counterstaining in 50% hematoxylin was reduced from 2 min to 30 s. Probes are listed in the resource table.  For fluorescent imaging, the pinhole was set to 1 airy unit for every channel, which were scanned sequentially to avoid overlapping emission spectra or with a combination of the ultraviolet and infrared channel in one track, where emission spectra overlap is minimal.
For the validation of the correct expression of the eGFP-L10a transgene, the 40× objectives on the LSM 700 or 710 were used with a zoom of 0.6× to acquire a tile scan over the complete lumbar dorsal horn (Pixel size: × 0.346 μm, y 0.346 μm, z 1.628/2.487 μm). A small z-stack (usually, 3.3 μm within 3 optical sections) at a bit depth of 8 was acquired. For the analysis, the cell counter plugin of ImageJ was used. Manual counting was done in the middle optical section; the other sections were used as an aid to distinguish between neighboring cells. Three hemi-sections from three animals of each mouse line were analyzed. Ratios were calculated per animal and then averaged.
For the quantification of the co-expression of the neuropeptides with multiplex FISH, the 40× objective on the LSM 800 was used with a zoom of 2× and an image size of 1024 × 1024 px (pixel size: × 0.078 μm, y 0.078 μm). An array of non-overlapping images that covers the complete dorsal horn was scanned at a bit depth of 16. Analysis was done with CellProfiler 2.2.0 software 52 . The cell profiler pipeline can be provided upon request. Briefly, nuclei were detected using the IdentifyPrimaryObjects module and expanded with the IdentifySecondaryObjects module to define the cell soma (from now on referred to as "cell"). The signal dots were detected using the IdentifyPri-maryObjects module, adapting the "threshold correction factor" and the "lower and upper bounds on threshold" www.nature.com/scientificreports/ for every slide (neuropeptide combination) to minimize false positive and false negative detection. The signal dots of the different channels were assigned to the cells they lie in using the RelateObjects module. The number of cells and related signal dots (per channel) were exported using the ExportToSpreadsheet module. Images with artefacts or high background were manually excluded from the automatic analysis in order to avoid false positive signal. Per slide (neuropeptide combination), three animals with two hemi-sections per animal were analyzed. Additionally, per combination one hemi-section hybridized with the 3-plex negative control probe and amplified with the same Amp 4 Alt was imaged with the same microscope setting and analyzed with the same CellProfiler setting in order to determine the background. Data analysis was conducted with R. To determine the signal background, the 3-plex negative control was analyzed. For every combination and channel, the 90th percentile of dots per cell in the negative control was determined. This value was set as the threshold above which a cell is counted as positive for the respective signal. 20 dots per cell were set as a threshold for high expression. The number of single, double and triple-positive cells was counted. Ratios were calculated per animal and then averaged.
Cell-type-specific polysomal mRNA isolation from excitatory and inhibitory spinal cord neurons. To specifically isolate polysomal mRNA from genetically targeted cells, the three bacTRAP mouse lines described above were used and the protocol published by Heiman, et al. 53 was followed. In brief, solutions and the affinity matrix were prepared as described. For the affinity matrix of Gad67::bacTRAP 300 μL and of vGluT2::bacTRAP and vGAT::bacTRAP 150 μL Streptavidin MyOne T1 Dynabeads and the corresponding volumes of biotinylated protein L and GFP antibodies 19C8 and 19F7 were used in a final volume of 200 μL.
Per mouse line three male mice per mouse line, aged between 8 and 10 weeks were used. The vertebral column containing the lumbar part of the spinal cord was dissected, immediately after decapitation. A syringe filled with ice-cold dissection buffer was fitted tightly against the caudal opening of the spinal canal and pressure was applied to press the spinal cord out, which was trimmed to the lumbar part. For the neuropathic pain analysis, the lumbar spinal cord was bisected along the midline and the left side (ipsilateral to the peripheral nerve surgery) was further processed. The tissue was washed, homogenized, centrifuged, lysed and again centrifuged as described in the protocol by Heiman et al. 53 .
Immunopurification was done as described by Heiman et al. 53 . RNA purification was done using the RNeasy Plus Micro Kit (Qiagen, catalog no. 74034). RNA was eluted from beads with 350 μl Buffer RLT Plus from the kit, supplemented with 40 μM dithiothreitol (DTT, Thermo Fisher Scientific, catalog no. R0861) and then processed according to the manufacturer's protocol. RNA concentration was determined using the RiboGreen fluorescence assay (Thermo Fisher Scientific, catalog no. R11491) on a NanoDrop 3300 microvolume fluorospectrometer. RNA quality was assessed by high sensitivity RNA screen tape on an Agilent 2200 TapeStation. RNA samples with RINe values of 6.8-7.8 were used for sequencing library construction.
Sequencing library construction, next-generation sequencing (NGS), mapping and differential gene expression analysis (DGEA). The sequencing library was prepared using the Ovation Mouse RNA-Seq System 1-16 kit (NuGEN, catalog no. 0348-32). The manufacturer's protocol was followed. 13-20 ng RNA were used as input material for the inhibitory-excitatory comparison and 30 ng RNA for the neuropathic pain analysis. The optional cDNA fragmentation step was performed on a Covaris S220 focused-ultrasonicator (Covaris) to obtain a median fragment size of approx. 300 bp. As recommended, the optimal cycle number for amplification was determined by qPCR (cycle number within the exponential phase of the amplification plot) and turned out to be 18.
Sequencing of the library, mapping to the reference genome, as well as DGEA were performed by the Functional Genomics Center Zurich (FGCZ). NGS was done on two lanes of a HiSeq 2500 system (illumina) with single-end sequencing of 125 bp. Mapping was done in R using the STAR package, version 2.4.2a [68]. The reference genome build Mus_musculus Ensembl GRCm38 release 80 was used for alignment. Reads mapped to exons were counted using countOverlaps in GenomicRanges version 1.20.5 (parameters: countNonredundant = TRUE, minFeatureOverlap = 10, minMapQuality = 10, keepMultiHits = TRUE) 54 . DGEA was done in R using either the package edgeR 3.10.2 or DESeq2 1.8.1 [69-71]. Bioinformatic analysis was performed using the R package ezRun (https ://githu b.com/uzh/ezRun ) within the data analysis framework SUSHI 55 . Quality checkpoints 56 , such as quality control of the alignment and count results, were implemented in ezRun (https ://githu b.com/uzh/ezRun ) and applied throughout the analysis workflow to ensure correct data interpretation. In details, read alignments were comprehensively evaluated using the RnaBamStatsApp (https ://githu b.com/uzh/ezRun /blob/maste r/R/app-RnaBa mStat s.R) in the ezRun package, in terms of different aspects of RNA-seq experiments, such as sequence quality, gDNA and rRNA contamination, GC/PCR/sequence bias, sequencing depth, strand specificity, coverage uniformity and read distribution over the genome annotation. Gene counts were comprehensively evaluated using the CountQCApp (https ://githu b.com/uzh/ezRun /blob/maste r/R/app-count QC.R) in the same R package, including sample correlation, sample clustering, as well as clustering of high variance features.
Gene Ontology (GO) term and pathway enrichment analyses. GO term enrichment analysis was done using the online PANTHER tool (http://panth erdb.org/). Genes significantly (FDR ≤ 0.05) and ≥ 2-fold enriched in excitatory (vGluT2) over inhibitory (vGAT or Gad67) and vice versa were analyzed. In the PAN-THER overrepresentation test, the sets of enriched genes were compared to the set of all expressed genes as a reference. GO molecular function (MF) complete was used as the annotation data set. The significantly overrepresented (FDR ≤ 0.05) MFs were analyzed. When several GO terms of a hierarchy were significantly overrepresented, only the most specific GO term was displayed in the graphical representation. For the analysis of Scientific Reports | (2021) 11:5232 | https://doi.org/10.1038/s41598-021-84667-y www.nature.com/scientificreports/ Gad67 vs. vGluT2, two redundant GO terms were not displayed in the graphical representation in order to allow comprehensibility. Pathway analysis was done with Ingenuity Pathway Analysis (IPA, QIAGEN). Genes significantly (FDR ≤ 0.05) and ≥ 2-fold enriched in excitatory (vGluT2) over inhibitory (vGAT or Gad67) and vice versa were analyzed. To get a coarse classification of the enriched genes into subcellular location and function ("type"), these parameters were exported from the software.