Transcriptomic evidence for dense peptidergic networks within forebrains of four widely divergent tetrapods

The primary function common to every neuron is communication with other neurons. Such cell-cell signaling can take numerous forms, including fast synaptic transmission and slower neuromodulation via secreted messengers, such as neuropeptides, dopamine, and many other diffusible small molecules. Individual neurons are quite diverse, however, in all particulars of both synaptic and neuromodulatory communication. Neuron classification schemes have therefore proven very useful in exploring the emergence of network function, behavior, and cognition from the communication functions of individual neurons. Recently published single-cell mRNA sequencing data and corresponding transcriptomic neuron classifications from turtle, songbird, mouse, and human provide evidence for a long evolutionary history and adaptive significance of localized peptidergic signaling. Across all four species, sets of at least twenty orthologous cognate pairs of neuropeptide precursor protein and receptor genes are expressed in individually sparse but heavily overlapping patterns suggesting that all forebrain neuron types are densely interconnected by local peptidergic signals.


Introduction
The primary function of a neuron can be defined most simply in terms of direct communication with some modest number of other neurons and neuroglial cells.
Such communication may involve transmission at chemical or electrical synapses, ephaptic coupling via extracellular ionic changes or current flows, and modulation of synapses, membrane excitability and other aspects of neuronal function via secreted neuropeptides (NPs), monoamines or other small molecules. Small subsets of neurons may also transduce sensory stimuli (e.g. light, touch, chemical) or control muscle or gland cells. An animal's sensory, cognitive, and motor functions emerge as communication between individual neurons iterates through the extended networks of neurons and glia that compose entire nervous systems [1]. One major challenge central to systems neuroscience is thus to delineate synaptic and neuromodulatory network connectivity graphs. Here, we will explore new perspectives on neuromodulatory network architectures offered by the recent introduction of sequencing methods for single-cell assay of gene expression at the level of messenger RNA (mRNA).
The neurons composing most animal nervous systems are highly heterogeneous. Neuronal heterogeneity was first glimpsed by nineteenth-century microscopists and has become ever more fascinating with each advance in methods for anatomical, physiological, molecular, and computational analysis of individual neurons [2e8]. Unfortunately, because nervous systems commonly comprise formidably vast numbers of diverse individual neurons and connections, the prospects for cataloging functional diversity at the level of individual neurons seem very remote indeed, except perhaps for the very simplest animals. Neuron-type classification thus looms as the more practical approach, many classification schemes based on neuronal form, firing and chemical composition have been proposed over the decades. Such 'neurotaxonomies' are at the heart of most ongoing attempts to develop useful neuronal network connectivity graphs.
Single-cell mRNA sequencing (scRNA-Seq) methods are now driving the development of neurotaxonomies of unprecedented depth, precision, and analytical utility [2,4,7,9e13]. By capturing genome-wide patterns of differential gene expression at the level of single cells, scRNA-Seq data and taxonomies are now revealing molecular mechanisms that regulate neuronal differentiation and network development [8,14e17] and enabling powerful new network physiological analyses based on abilities to probe and perturb specific neuron types [2,18]. The construction of modulatory network graphs is one area where single-cell transcriptomes and neurotaxonomies are having an especially direct impact: the scRNA-Seq data themselves provide strong predictions regarding the differential expression of the 'effector' molecules (e.g. ligands and receptors) that determine the form, effects and typespecific targeting of neuromodulatory signaling [13,19e23]. Though critical anatomically diffuse modulatory actions of both neuropeptides (e.g. hypothalamic peptides) and monoamines (e.g. dopamine and serotonin) across broad brain areas are widely recognized [24], the new single-cell data suggest an unexpectedly pervasive presence of much more localized yet highly neuron-type-specific forms of peptidergic signaling. The present article will elaborate upon this prospect, focusing especially on emerging transcriptomic evidence for highly neuron-typeespecific peptidergic signaling localized to tightly circumscribed brain areas.

Modulation of synaptic network function by neuropeptides
'Neuropeptide' refers to a broad set of neuroactive peptide fragments derived from large numbers of distinct neuropeptide precursor proteins (NPPs) [25,26] that are fragmented, covalently modified, and stored within secretory vesicles, then secreted in response to neuronal electrical activity. Although definitions of 'neuropeptide' vary, one of the strictest requires that the ligand in question act via a G-proteinecoupled receptor (GPCR) [27]. Large numbers and the high conservation of genes encoding NPPs and neuropeptide-selective GPCRs (NP-GPCRs) across all modern eumetazoan phyla speaks powerfully to their fundamental importance to nervous system function [28e31]. Though GPCR genes comprise a very large family (on the order of 1000 members), the encoded receptors exert their effects by way of a very limited number of direct signal transduction processes, primarily via members of the eponymous heterotrimeric G protein family [32,33]. Though certainly an oversimplification (e.g. the study by Jung et al. [34]), actions of neuropeptides can be classified to a first-order according to the G-protein alpha subunit (G a ) preference of the 'cognate' NP-GPCR selective for that particular NP. For NP-GPCRs expressed in the brain areas considered here, primary subunit preference falls into just three classes, G a i/o (inhibits cAMP production), G a q/11 (augments calcium signaling via phospholipase C activation), and G a s (stimulates cAMP production).
The impacts of NP signaling upon neuronal function are profound yet highly neuron-type-dependent [22]. Both cAMP and calcium, the two main second messengers regulated by NP-GPCR activation, are known to potently regulate ion channel function via protein phosphorylation, while ion channels are, in turn, the major substrates of neuronal excitability and major determinants of synaptic transmission [24]. Though just three classes of direct G a effects upon cAMP and calcium signaling are recognized, the downstream coupling between those second messengers and membrane ion channels is multi-faceted and far more diverse, depending upon neuron-typeedependent differential expression of dozens of protein kinase genes and hundreds of ion channel genes. The powerful effects of a given NP thus depend not only upon the corresponding NP-GPCR's G a subunit preference but also upon highly neuron-type-dependent profiles of protein kinase and ion channel gene expression.
Ancestors of all modern eumetazoans likely accomplished their cellecell signaling primarily via diffusional 'broadcast' of one cell's secreted chemical messenger (often a peptide) over varying distances to GPCRs on a second, target cell [30,32,35e37]. There are even strong phylogenomic hints that such 'paracrine' signaling predated the evolutionary invention of neurons with their axons, dendrites, and focal synaptic connections [7,21,28,29,37e40]. Abundant expression of many different NPP and NP-GPCR genes has now been reported in most or all neurons studied by single-cell transcriptomics across a wide range of eumetazoan phyla [21, 36,41]. Given that most or all neurons of present-day animals also use at least one small molecule neurotransmitter (e.g. glutamate or GABA), it follows that most neurons probably signal via both neuropeptides and a small molecule transmitter [27,42]. Individual neurons, moreover, often express multiple NPPs and multiple NP-GPCRs, suggesting that a combinatorial NP code may serve to support very dense yet directed and specific addressing of NP signals from one neuron (or neuron type) to another [20,21].
Though fast synaptic networks are clearly fundamental to the fast reactions and behaviors of contemporary eumetazoans, critical parameters such as synaptic strength, dendritic integration, and membrane excitability must be continually and critically adjusted to support the homeostasis, fine-tuning, flexibility, and plasticity of nervous system function. Such adjustment is the province of neuromodulation, adaptive nervous system function clearly depends on a well-coordinated interplay of fast synaptic and slower neuromodulatory cellecell signaling [28,43e45]. From this recognition and results reviewed here, there emerges a view of the neuronal network as a superposition of many synaptic and modulatory networks, integrated by individual neurons as common network nodes. To explore evolutionary pressures that may have shaped and preserved such intimately mingled synaptic and modulatory architectures, this commentary will examine results from focused analysis of published Table 1 The present study is based on published single-cell mRNA-Seq data acquired from the four tetrapod species and brain areas tabulated here, all using droplet-based RNA-Seq platforms [68]. The listed publications describe the corresponding transcriptomic and taxonomic methods in detail and link to the indicated publicly available datasets upon which the present analysis is based. Bakken et al. [49] https://portal.brain-map.org/atlases-and-data/rnaseq/human-m1-10x Table 2 The thirty-two neuropeptide signaling genes (11 NPPs and 21 NP-GPCRs) tabulated here are transcribed at high levels in most or all of the four tetrapod species and forebrain areas considered. Each row lists one NPP gene, active neuropeptide product(s) derived from the encoded neuropeptide precursor protein, and one to five cognate NP-GPCR genes (i.e. genes that encode receptors selective for one or more products of the encoded NPP). Each NPP/NP-GPCR pair defines a potential directed signaling link from an NPP-expressing source cell to an NP-GPCR-expressing target cell. Accordingly, this table represents a total of 24 potential NPP to NP-GPCR potential connections. Symbols Tac2 (mouse) and TAC3 (turtle, bird, human) refer to proper orthologs. Color shading of NP-GPCR symbols (bottom right) indicates the expected primary GPCR signal transduction mode.

Single-cell transcriptomic neurotaxonomies and data resources
Data from four recent publications [46e49] offer novel transcriptomic perspectives on local peptidergic signaling within forebrains of four divergent tetrapods (see Table 1). Each publication describes the collection of a large single-cell mRNA-Seq dataset and the development of a corresponding neurotaxonomy. These four outstanding publications were chosen to represent a broad range of amniote tetrapod vertebrates (reptilian, avian and mammalian) with scRNA-Seq transcriptomes sampled by comparable methods, well-established gene orthology, and published neurotaxonomies. Results

Figure 1
Statistics of single-cell NPP and NP-GPCR mRNA counts by a subclass, brain region and species. (a) Neurotaxonomic classifications based on inference from single-cell RNA-Seq data, as described in the four Table 1 Table 1 and each gene as listed and color coded in Table 2. Dashed gray horizontal lines indicate that mRNA for that ortholog was not detected for that species and gene. Mean counts were normalized to the maximum for each gene/species row, represented using a color scale at the lower right. Numbers at the right of maps for each species indicated rank order of row-maximum mean counts. Lists of both NPP and NP-GPCR genes are    16 for mouse and 14 for human). I will refer to these categories as 'subclasses'. In the case of GABAergic neurons, it is also now possible to infer a cell's embryological neurogenesis region, or 'birthplace', from transcription factor expression patterns. I will refer to coarse subdivisions of the major GABA neuron class into 'minor classes' based on transcriptomically inferred brain region of last division, LGE (lateral ganglionic eminence), medial ganglionic eminence or CGE (caudal ganglionic eminence). Here, the terms 'type' and 'typespecific' are used in generic senses rather than in reference to any specific taxon.
Type-specific expression of cognate NPP and NP-GPCR genes in four divergent tetrapods Table 2 lists 11 NPP and 21 cognate NP-GPCR genes for which corresponding mRNA was detected at high levels in neurons of most or all Table 1 species and brain areas. Each row begins with a symbol for one of the NPP genes, then calls out the active neuropeptide product (or products) of that gene and one to five highly expressed cognate NP-GPCR genes, with color fills indicating G a subunit preference of each. The abundant expression of these 11 NPP and their 21 cognate NP-GPCR genes, all within tightly circumscribed brain areas of each species, suggests the likelihood of local, paracrine peptidergic signaling. Transcripts of additional NPP and NP-GPCR genes, beyond the 32 listed in Table 2, were detected in each species and area, but because they do not align as cognate pairs, they seem more likely to mediate signaling to or from more distant regions (e.g. endocrinelike actions of hypothalamic neuropeptides, or afferent or efferent projections from or to distant neurons). I focus here, therefore, on the 32 Table 2 genes, following a special interest in the relatively unexplored prospect of dense and localized paracrine signaling. Figure 1 displays neurotaxonomic expression maps representing scRNA-Seq data, averaged per subclass, for each of the 32 Table 2 genes and each of the Table 1 species. Though cross-species alignment of subclasslevel neurotaxonomies remains a work in progress (as discussed in each of the four parent publications and elsewhere [6e8,50]), alignment of major and minor classes, as in Figure 1, is straightforward (though the lateral ganglionic eminence GABA neuron minor class was detected only for turtle and bird samples). For most species studied so far, peak scRNA-Seq counts for numerous NPPs rank among the very highest across all neuronal protein-coding genes [19e21], as perhaps to be expected for genes whose products are secreted rather than incorporated into relatively stable cellular components. Seven of the listed NPP genes (NPY, SST, PENK, CCK, TAC1, VIP and ADCYAP1) show very high peaksubclass transcript counts across all four species, while expression levels of CRH, PDYN, TAC2/3 and PNOC are low or insignificant in one or more of the four. The expression levels of the Figure 1 NP-GPCR genes are much lower than those of most NPPs but nonetheless fall comfortably within the range of other genes encoding other neuronal membrane proteins known to be strongly functional. An extended version of Figure 1, which includes more numerical expression details and statistics is available at https://github.com/profsjsmith/tetrapods.
Perhaps the most striking feature common to all Figure 1 expression maps is the strong variation across classes and subclasses for each NPP and NP-GPCR gene and the fact that such variation is patterned differentially for each gene. Both similarities and differences in patterning amongst species are notable in Figure 1. For all four tetrapods, more NPP genes show higher and row vectors representing the expression of cognate NPP and NP-GPCR genes, respectively, calculated from single-cell RNA-Seq data and color scaled exactly as for individual gene rows in Figure 1. (a4) An adjacency matrix is calculated as the outer product of the NPP column vector (taken as a proxy for cellular NP secretion rate) and the NP-GPCR row vector (taken as a proxy for cellular NP responsivity). Predicted signaling weight from column to row neuron types, normalized to matrix maximum, as indicated by color scale at left. (a5) Overlaid lines demarcate major (magenta) and minor (cyan) classes as in Figure 1 and expression in GABA neurons than in glutamate neurons, while the opposite is true for NP-GPCR genes. Figure 1 also shows that subclass specificity of NPP expression is overall considerably sharper than that of NP-GPCR expression. While cross-species comparison of expression of particular genes at the level of subclasses is problematic as noted previously, Figure 1 shows that neurotaxonomic patterning for most of these 32 NP genes still shows numerous obvious similarities across species at the minor and major class levels. These similarities across widely divergent species argue for substantial adaptive value of the observed patterning of NPP and NP-GPCR gene expression.
Transcriptomic prediction of local peptidergic network architectures RNA-Seq data permit a loosely quantitative prediction of peptidergic signaling directed from a 'source' neuron secreting one particular neuropeptide to a cognate receptor-bearing 'target' neuron [19,22]. The logic here is premised upon (1) taking the source neuron's NPP transcript counts as a surrogate for the time-averaged rate at which that NPP's active neuropeptide product is secreted, (2) taking the target neuron's cognate NP-GPCR counts as a surrogate for that neuron's sensitivity to that neuropeptide, (3) assuming paracrine 'broadcast' diffusion of the neuropeptide from source to any nearby target, and then (4) multiplying the source neuron's NPP counts by the target neuron's NP-GPCR counts to predict signaling strength. Though this prediction logic involves numerous assumptions unlikely to hold with much precision and neglects many other factors certain to be of cardinal importance, the resulting predictions are likely nonetheless to be of value as experimentally testable working hypotheses [22] and bases for theoretical exploration [51].
The prediction logic outlined above is perhaps most compelling for the case where a GABAergic neuron is the NP source. Compared to glutamatergic neurons, GABA neurons are generally the more prolific in NPP gene expression, and most GABA neurons lack distant axonal projections (but see the study by Melzer and Monyer [52]). This fact and the limited lifetimes of NPs once secreted make distant actions of peptide release by GABA neurons unlikely [21]. The extraordinarily high and widespread abundance of many NPP gene transcripts in GABA neurons (and the fact that every transcript comes at a substantial metabolic cost) therefore argues that NP products must play important physiological signaling roles [19,21] because no other roles for NPP gene products are known. Nonetheless, similar d if just slightly less compelling d arguments still can be made on behalf of glutamate neurons as local peptidergic signal sources.
The prediction logic described above can be cast in neurotaxonomic matrix form to generate network-level peptidergic connectivity hypotheses. For a given cognate NPP/NP-GPCR pair and a given species, NPP counts for each neuron subclass can be cast as a column vector and the cognate NP-GPCR's counts as a row vector. (Each of these vectors is identical as displayed in one gene/species row in the Figure 1 expression map, except with transpose of the NPP row to form a column). An adjacency matrix expressing all subclass-tosubclass signaling predictions for any cognate pair can then be calculated as the outer product of those two vectors. Figure 2A lays out this prediction scheme in graphic form. The neurotaxonomic NP gene expression data of Figure 1 and the 24 cognate NPP/NP-GPCR pairings of Table 2 allow the prediction of many such adjacency matrices for all four species. Figures 2B-2E displays a sampling of six such adjacency matrix predictions for each species as a representative sample; the remainder can be accessed at https://github.com/ profsjsmith/tetrapods.
As might be expected from the sparsity of the Figure 1 expression maps, the individual adjacency matrix outer products represented in Figure 2 are themselves even more sparse. Though each predicted peptidergic network connects only sparse subsets of cell types, the complete sets of predicted NP networks nonetheless provide dense modulatory connectivity; they interlink all major and minor classes and subclasses. Strong similarities among the network predictions for the four widely divergent species are moreover clear from inspection of the matrices represented in Figure 2.
Matrix elements highlighted along the diagonals in Figure 2 indicate subclasses where the adjacency matrix predicts autocrine signaling, that is, where neurons classified within a given subclass express both a given NPP and its cognate NP-GPCR. Autocrine neuropeptide signaling is suggested for many other systems [21,38], so evidence of the existence of such signaling in the brain is certainly of interest. Cognate co-expression predicts autocrine signaling involving NP-GPCRs in all three G protein alpha subunit preference classes, so a wide range of negative and positive feedback linkages are possible, including links that could mediate cell-autonomous homeostasis [53] or regenerative amplification of peptide secretion [21,38]. Because there is good evidence for dendritic NP secretion [42], it is furthermore possible that autocrine NP actions could be localized to individual dendrites or dendritic branches. By contributing to dendrite-or branch-specific integration, plasticity, and/or homeostasis [54], autocrine NP signaling could thereby enhance the information processing capacity inherent in individual neurons [55]. The similarities of typespecific autocrine NP expression patterns across the four divergent species evident in Figure 2, moreover, suggest that autocrine NP signaling has very likely been of significant advantage to all four throughout their evolution.
The highly abundant and type-specific co-expression of numerous cognate NPP and NP-GPCR genes within confined brain regions suggests the existence and functional importance of dense local neuromodulatory networks. For each species and brain area considered here, the predicted peptidergic networks aggregate to densely interconnect all or very nearly all neuronal subclasses. Figure 2 furthermore highlights striking similarities in the type-specific patterning of such connectivity across widely divergent tetrapods, most obvious at the levels of major and minor class but also discernible for subclasses, despite imperfect alignment of subclasses across species. Regardless of whether these similarities reflect evolutionary conservation or convergence, they suggest that dense local peptidergic networks are probably of ancient yet still very substantial adaptive value.

Outlook
The above analysis casts single-cell mRNA-Seq data from circumscribed forebrain areas of divergent tetrapods as local modulatory network predictions and reflects upon possible implications of these predictions for brain function and evolution. A multitude of questions remains, of course, regarding the many underlying assumptions and the correspondence of these predictions to actual nervous system physiology. Certainly, limitations of present scRNA-Seq methods must be considered and quantitative connections between transcript abundance, rates of translation, protein abundance, processing, and trafficking all must be investigated much more deeply. To explore and understand the actual Spatio-temporal dynamics of local neuropeptide signaling, it will also be essential to explore the diffusion and enzymatic inactivation of NP ligands in brain interstitial spaces, via both experiments and theory [21]. Moreover, given the fundamentally recurrent, overlapping, and redundant nature of neuronal networks, the rational design and interpretation of experiments to fathom neuropeptide impacts upon network function remain as formidable challenges. Fortunately, new physiological tools (e.g. the study by Smith et al. [22], Melzer et al. [45], Tejeda et al. [56], and Vu et al. [57]) and new theoretical and computational modeling approaches (e.g. the study by Jekely [21], Nusbaum et al. [44], Liu et al. [51], and Reddy et al. [58]) seem ideally suited to address many of these hard outstanding challenges are now at hand.
Further refinement of brain neurotaxonomies and their more precise and reliable alignment across divergent species, brain regions, and data acquisition modalities remain critically important broad goals for neuroscience. It will be of great interest to see what more advanced and tightly reconciled neurotaxonomies reveal both conservation and drift in the evolution of peptidergic network architectures. It also will be of extraordinary interest to see what innovations in multi-modal data acquisition and neurotaxonomy reveal about the alignment of neuromodulatory and synaptic connectivity graphs and how such new insights advance future theoretical explorations of adaptive brain function.
Similarities of highly distinctive patterns of NPP and NP-GPCR gene expression across four widely divergent tetrapods argue for the deep functional importance of local peptidergic signaling. The full cladistics of NPP and NP-GPCR genes further suggest that local peptidergic modulation may govern the function of all eumetazoan nervous systems [16,21,28,29,40,59e61]. Such premises surely will be further tested, modified, and elaborated as parts of ongoing experimental and theoretical explorations of single-cell transcriptomics in species ever more broadly and deeply spanning the metazoans (e.g. the study by Taylor et al. [20], Jekely [21], Crocker et al. [41], Sebe-Pedros et al. [62], Croset et al. [63], Deng et al. [64], and Kozma et al. [65]). Finally, it seems highly likely that such broadened and deepened exploration of neuropeptide-related gene expression will illuminate presently unknown general principles of circuit homeostasis, modulation, and plasticity. Such light may, in turn, bring us closer to finally fulfilling the long-awaited promise of neuropeptides and their receptors as targets for better neurological and psychiatric therapeutics [66,67].

Conflict of interest statement
Nothing declared.