A multidimensional coding architecture of the vagal interoceptive system.

Mice

All animal husbandry and procedures were performed in compliance with Yale University’s Institutional Animal Care and Use Committee and National Institute of Health (NIH) guidelines. All mice were age- and gender-matched adults (older than 8 weeks) and no differences between sexes were observed. Mouse lines Wild-type C57BL/6J (000664), Glp1r-ires-Cre (029283), Gpr65-ires-Cre (029282), Npy2r-ires-Cre (029285), P2ry1-ires-Cre (029284), Agtr1a-Cre (030553), Calb2-ires-Cre (010774), Nts-Cre (017525), Vglut2-ires-Cre (016963), Sst-ires-Cre (013044), Twist2-Cre (008712), Vip-ires-Cre (010908), Piezo2-GFP-ires-Cre (027719), Trpv1-Cre (017769), Pvalb-Cre (017320), Vglut1-ires2-Cre (023527), Chat-ires-Cre (031661), lox-ChR2 (024109), lox-tdTomato (007914), and Snap25-2A-GCaMP6s-D (025111) were from the Jackson Laboratory. Drd2-Cre (032108-UCD) mice were from the Mutant Mouse Resource and Research Center (MMRRC). lox-L10-GFP mice were described before 11 , 12 .

Generation of Projection-seq AAVs

UPB sequences were cloned from coding sequences of hChR2 (Addgene, plasmid 28017), hM3Dq (Addgene, plasmid 50474) and hM4Di (Addgene, plasmid 50475) and inserted right before the SV40 poly (A) of AAVrg-CAG-tdTomato-WPRE-SV40 (Addgene, plasmid 59462) using the In-Fusion HD Cloning Kit (Takara, 638909). Projection-seq AAVs (at titre > 10 12 –10 13 viral genomes per ml) were generated at the UNC Vector Core. Plasmids have been deposited to Addgene. UPB sequences (5′–3′): UPB-oesophagus (UPB5, Addgene plasmid 180783): ACAGCACCATCCTCAACTCCACCAAGTTACCCTCATCGGACAACCTGCAGGTGCCTGAGGAGGAGCTGGGGATGGTGGACTTGGAGAGGAAAGCCGACAAGCTGCAGGCCCAGAAGAGCGTGGACGATGGAGGCAGTTTTCCAAAAAGCTTCTCCAAGCTTCCCATCCAGCTAGAGTCAGCCGTGGACACAGCTAAGACTTCTGACGTCAACTCCTCAGTGGGTAAGAGCACGGCCACTCTACCTCTGTCCTTCAAGGAAGCCACTCTGGCCAAGAGGTTTGCTCTGAAGACCAGAAGTCAGATCACTAAGCGGA. UPB-stomach (UPB1, Addgene plasmid 180784): ATGGACTATGGCGGCGCTTTGTCTGCCGTCGGACGCGAACTTTTGTTCGTTACTAATCCTGTGGTGGTGAACGGGTCCGTCCTGGTCCCTGAGGATCAATGTTACTGTGCCGGATGGATTGAATCTCGCGGCACGAACGGCGCTCAGACCGCGTCAAATGTCCTGCAGTGGCTTGCAGCAGGATTCAGCATTTTGCTGCTGATGTTCTATGCCTACCAAACCTGGAAATCTACATGCGGCTGGGAGGAGATCTATGTGTGCGCCATTGAAATGGTTAAGGTGATTCTCGAGTTCTTTTTTGAGTTTAAGAATCCCTCTATGCTCTACCTT. UPB-duodenum (UPB6, Addgene plasmid 180785): AATGGCAGCTCGGGCAATCAGTCCGTGCGCCTGGTCACGTCATCATCCCACAATCGCTATGAGACGGTGGAAATGGTCTTCATTGCCACAGTGACAGGCTCCCTGAGCCTGGTGACTGTCGTGGGCAACATCCTGGTGATGCTGTCCATCAAGGTCAACAGGCAGCTGCAGACAGTCAACAACTACTTCCTCTTCAGCCTGGCGTGTGCTGATCTCATCATAGGCGCCTTCTCCATGAACCTCTACACCGTGTACATCATCAAGGGCTACTGGCCCCTGGGCGCCGTGGTCTGCGACCTGTGGCTGGCCCTGGACTGCGTGGTGAGCAACGCCTCCGTCATGAACCTTCTCATCATCAGCTTTGACCGCTACTTCTGCGTCA. UPB-colon (UPB4, Addgene plasmid 180786): ATCGATGGGCCTTAGGGAACTTGGCCTGTGACCTCTGGCTTGCCATTGACTGCGTAGCCAGCAATGCCTCTGTTATGAATCTTCTGGTCATCAGCTTTGACAGATACTTTTCCATCACGAGGCCGCTCACGTACCGAGCCAAACGAACAACAAAGAGAGCCGGTGTGATGATCGGTCTGGCTTGGGTCATCTCCTTTGTCCTTTGGGCTCCTGCCATCTTGTTCTGGCAATACTTTGTTGGAAAGAGAACTGTGCCTCCGGGAGAGTGCTTCATTCAGTTCCTCAGTGAGCCCACCATTACTTTTGGCACAGCCATCGCTGGTT. UPB-pancreas (UPB7, Addgene plasmid 180787): AAGATGGCAGGCCTCATGATTGCTGCTGCCTGGGTACTGTCCTTCGTGCTCTGGGCGCCTGCCATCTTGTTCTGGCAGTTTGTGGTGGGTAAGCGGACGGTGCCCGACAACCAGTGCTTCATCCAGTTCCTGTCCAACCCAGCAGTGACCTTTGGCACAGCCATTGCTGGCTTCTACCTGCCTGTGGTCATCATGACGGTGCTGTACATCCACATCTCCCTGGCCAGTCGCAGCCGAGTCCACAAGCACCGGCCCGAGGGCCCGAAGGAGAAGAAAGCCAAGACGCTGGCCTTCCTCAAGAGCCCACTAATGAAGCAGA. UPB-lung (UPB2, Addgene plasmid 180788): ACAGGACACCGGGTGCAGTGGCTGCGCTATGCAGAGTGGCTGCTCACTTGTCCTGTCATCCTTATCCGCCTGAGCAACCTCACCGGCCTGAGCAACGACTACAGCAGGAGAACCATGGGACTCCTTGTCTCAGACATCGGGACTATCGTGTGGGGGGCTACCAGCGCCATGGCAACCGGCTATGTTAAAGTCATCTTCTTTTGTCTTGGATTGTGCTATGGCGCGAACACATTTTTTCACGCCGCCAAAGCATATATCGAGGGTTATCATACTGTGCCA. UPB-heart (UPB8, Addgene plasmid 180789): AGGACACTTCCAATGAGTCCAGCTCAGGCAGTGCCACCCAGAACACCAAGGAACGCCCAGCCACAGAGCTGTCCACCACAGAGGCCACCACGCCCGCCATGCCCGCCCCTCCCCTGCAGCCGCGGGCCCTCAACCCAGCCTCCAGATGGTCCAAGATCCAGATTGTGACGAAGCAGACAGGCAATGAGTGTGTGACAGCCATTGAGATTGTGCCTGCCACGCCGGCTGGCATGCGCCCTGCGGCCAACGTGGCCCGCAAGTTCGCCAGCATCGCTCGCAACCAGGTGCGCAAGAAGCGGCAGATGGCGGCCCGGGAGCGCAAAGTGACACGAACGATCTTTGCCATTCTGCTGGCCTTCATCCT.

AAV infection of visceral organs

Retrograde adeno-associated virus (AAVrg) 51 represents a powerful genetic tool in VSN studies 52 – 54 . AAVrg-tdTomato and AAVrg-GFP were purchased from the UNC vector core. AAVrg-FLEX-tdTomato (28306-AAVrg), AAVrg-CAG-FLEX-rc [Jaws-KGC-GFP-ER2] (84445-AAVrg) and pAAVrg-hSyn-Con/Fon hChR2(H134R)-EYFP (55645-AAVrg) were purchased from Addgene. All viruses contained 10 12 –10 13 viral genomes per ml, and 0.05% Fast Green FCF (Sigma-Aldrich, F7252-5G) was occasionally used to facilitate visualization. In all surgeries, mice were anaesthetized with 1–2% isoflurane on a heating pad, followed by subcutaneous injections of meloxicam (5 mg kg −1 ) and buprenorphine (1.5 mg kg −1 ). Lung: virus (5 μl diluted with 75 μl of saline) was injected through a tracheal catheter into the lung using a Hamilton syringe. Heart: the mouse was ventilated (tidal volume 0.21 ml, respiration rate 110 breaths per minute) using a mouse ventilator (SAR-1000, CWE) via an intubated angiocatheter. The heart was gently exposed via thoracotomy. Virus (5 μl) was injected intramurally (20 nl per second) using a Nanoject III injector at multiple sites covering most of the area of the heart. Stomach, duodenum, colon and pancreas: the target organ was gently exposed via an abdominal incision. Virus was injected intramurally at multiple sites covering most of the area (stomach: 5 μl; all subregions in both dorsal and ventral sides; duodenum: 2.5 μl, a length of 1.5 cm from the pyloric sphincter; transverse colon: 2.5 μl, a length of 1.5 cm; pancreas: 5 μl). Oesophagus: the cervical oesophagus underneath the trachea was surgically exposed via a neck incision. The abdominal oesophagus and the oesophageal sphincter were gently exposed via an abdominal incision. Virus was injected into oesophageal muscularis and serosa layers (cervical: 2.5 μl; abdominal: 1 μl, between the diaphragm and the oesophageal sphincter). To visualize and quantify the anatomical location of VSNs innervating various visceral organs in the nodose ganglion (Extended Data Fig. 1a–c ), the following AAVrgs were injected into visceral organs in wild-type mice for Extended Data Fig. 1a–c : heart (tdTomato)/lung (GFP), heart (tdTomato)/stomach (GFP), heart (tdTomato)/oesophagus (GFP), pancreas (tdTomato)/stomach (GFP), heart (tdTomato)/duodenum (GFP), and colon (tdTomato). To compare infection efficiencies of AAVrg and the conventional neural tracer cholera toxin B subunit (CTb), AAVrg-tdTomato (2.5 μl) and Alexa Fluor 647-conjugated CTb (1.0 mg ml −1 , 2.5 μl, Thermo Fisher Scientific, C34778 ) were mixed and co-injected into the stomach in wild-type mice (Extended Data Fig. 1d ). AAVrg and CTb labelled partially overlapping VSNs with similar efficiencies, raising the possibility that AAVrg might have some preferential tropism for specific VSN populations over CTb. However, it is worth noting that co-infection of VSNs with two similar AAVs (AAV9-FLEX-tdTomato, AAV9-GFP) via nodose injection also resulted in similar partially overlapping labelling 11 , suggesting that this likelihood is low. This possibility is further reduced by the extensive validation of Projection-seq results using anterograde tracing in 14 Cre mouse lines as described in Extended Data Figs. 4 – 8 . To verify the level of coverage of injected virus within the organ (Extended Data Fig. 1e–q, s–u ), AAVrg-tdTomato or Fast Green FCF (5% w/v, same volume as AAVrg) were injected into wild-type mice as described above. Organs were processed as described in ‘Histology and immunochemistry’. Injected dye or virus were able to cover most subregions and virus-infected cells were observed in most tissue layers including the pulmonary alveoli, stomach and intestinal villi, and the entire myocardial layer of the heart ventricles. We also noticed that the atria and bronchial walls were not well covered. To examine whether AAVrgs could label vagal motor neurons in the DMV, AAVrg-FLEX-tdTomato (2.5 μl) and AAVrg-GFP (2.5 μl) were mixed and co-injected into the stomach of Chat-ires-Cre mice (Extended Data Fig. 1r ). To examine the degree of viral spread (Extended Data Fig. 1v, w ), AAVrg-tdTomato and AAVrg-GFP were separately injected into indicated organs in the same mice. To determine the infection efficiency of Projection-seq AAVs (Extended Data Fig. 2a, b ), Projection-seq AAVs and AAVrg-GFP were injected into the stomach. For Projection-seq analysis (Fig. 1 ), Projection-seq AAVs with UPB-lung, UPB-heart, UPB-oesophagus, UPB-stomach, UPB-duodenum, UPB-colon and UPB-pancreas were injected into corresponding organs. The heart, cervical oesophagus and lung were sequentially injected on the first day, and abdominal oesophagus, stomach, duodenum, pancreas and colon were sequentially injected into the same mice on the second day. Approximately 30% of VSNs were labelled from the seven organs (Extended Data Fig. 2c ). For RNAscope HiPlex assays of retrogradely labelled VSNs (Extended Data Fig. 2p, q ), AAVrgs encoding reversed transcripts of fluorescent proteins, including AAVrg-FLEX-tdTomato (heart, stomach, duodenum), AAVrg-CAG-FLEX-rc [Jaws-KGC-GFP-ER2] (lung, oesophagus, pancreas) and AAVrg-hSyn Con/Fon hChR2(H134R)-EYFP (colon) were used to avoid potential contamination of RNAscope results from virus-introduced fluorescent signals. For tracing the central targets of diverse vagal pathways, AAVrg-tdTomato was injected into visceral organs of wild-type mice (Extended Data Fig. 13 ), and AAVrg-FLEX-tdTomato was injected into visceral organs of indicated Cre lines (Fig. 4 , Extended Data Fig. 14 ). Vagal ganglia were collected seven days after surgery and brains were collected two to three weeks after surgery. PCR analysis of VSNs Vagal ganglia were collected from control mice or mice with stomach injection of AAVrg-tdTomato-UPB2. RNA was extracted using the Monarch Total RNA Miniprep Kit (NEB, T2010S) and reverse-transcribed into cDNA using the SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific, 18091050). Primer sets used (Extended Data Fig. 2b ): UPB2 (5′–3′): CATCGATACCGTCGACACAGGACACCGGGTGCAGTG (forward); CTGCTCGAAGCGGCCGCTGGCACAGTATGATAACCCTCG (reverse). UPB4 (5′–3′): CATCGATACCGTCGACATCGATGGGCCTTAGGGAAC (forward); TGCTCGAAGCGGCCGCAACCAGCGATGGCTGTGCC (reverse). scRNA-seq and Projection-seq of VSNs innervating various visceral organs Neuron isolation and sequencing For control scRNA-seq, vagal ganglia (left and right) were collected from 40 age- and gender-matched C57BL/6J wild-type mice (10 mice per sample) and VSNs were acutely isolated and enriched using previously described methods 11 , 12 , 17 . Approximately 5,000–10,000 VSNs were loaded in each channel of the 10X microfluidic device to target 3,000–6,000 cells as an output of one sample. Single-cell cDNA libraries were prepared at the Yale Center for Genomic Analysis (YCGA) and sequenced using an Illumina NovaSeq S4 sequencer at 150–300 million reads to achieve a fine sequencing depth of 30,000–50,000 reads per cell. For Projection-seq, 30 age- and gender-matched wild-type mice (divided into four samples) were first injected with different Projection-seq AAVs into thoracic and abdominal organs (details in ‘AAV infection of visceral organs’). Vagal ganglia were collected seven days later. Of the examined VSNs, 29.5% (102/346) were tdTomato + . VSNs were sequenced as mentioned above.

Bioinformatic processing

Transcriptomic data were aligned to the mm10 mouse genome reference (control scRNA-seq) or a custom mouse genome reference with additional sequence information of UPBs (Projection-seq) using the Cell Ranger software v.3.0.2 (10X Genomics). The following quality control metrics were applied to filter low-quality cells in scRNA-seq: number of genes per cell > 500; number of genes per cell < 8,000; percentage of mitochondria genes < 10%. A total of 56,575 cells were sequenced (31,182 control scRNA-seq samples: 7,842, 7,580, 7,939, 7,821; 25,393 Projection-seq samples: 6,403, 6,381, 5,803, 6,806). Control scRNA-seq and Projection-seq data were then integrated and processed using the R package Seurat v.3 55 , and 42 cell clusters identified using the top 30 principal components (PCs) were visualized using UMAP 56 (Extended Data Fig. 2f ). A total of 32,558 neurons selected from 25 Slc17a6 + clusters (16,267 control scRNA-seq samples: 4,961, 3,894, 4,288, 3,124; 16,291 Projection-seq samples: 3,863, 3,431, 4,449, 4,548) were re-clustered into 67 populations with the top 100 PCs (Extended Data Fig. 2g , top). A total of 27,800 placode-derived neurons selected from 52 Phox2b + clusters (13,210 control scRNA-seq samples: 3,835, 3,091, 3,452, 2,832; 14,590 Projection-seq samples: 3,443, 2,973, 4,019, 4,155) were re-clustered into 52 clusters with the top 100 PCs and visualized with UMAP separately (Fig. 1b for Projection-seq; Extended Data Fig. 2g , bottom for control scRNA-seq) or together (Extended Data Fig. 2h ). DEGs for the 52 clusters were identified using the Wilcoxon rank-sum test implemented in Seurat from the Projection-seq dataset (Extended Data Fig. 2d ). VSN clusters were then manually grouped into 12 subpopulations on the basis of expression of DEGs and their locations on the UMAP plot. E-VSNs selectively express markers for damaged sensory neurons, such as Sprr1a and Ecel1 16 , 57 , 58 (Extended Data Fig. 3a ). The percentage of E-VSNs increased after the Projection-seq process (control: n = 627, 4.7%; Projection-seq: n = 1,528, 10.5%), and consistently more Sprr1a + VSNs were observed after AAVrg injection (Extended Data Fig. 3b–d ), demonstrating that neurons damaged during both cell dissociation and Projection-seq procedures can be easily identified. Thus, E-VSNs were removed from further analysis.

Projection-seq analysis

Further analysis was performed for the Projection-seq dataset. A total of 42 out of 1,701 Prdm12 + (2.5%) neural-crest-derived jugular VSNs and 4,791 out of 14,590 (32.8%) Phox2b + placode-derived nodose VSNs 18 , 59 were recognized as UPB-positive (expression level > 0.8), suggesting that VSNs retrogradely labelled from the seven major visceral organs in our study mainly originate from the nodose but not the jugular ganglia. We therefore focused on nodose VSNs. After removing E-VSNs, 3,539 out of 4,609 UPB-marked nodose VSNs (76.8%) expressed a single barcode (lung-UPB, 855; oesophagus-UPB, 595; heart-UPB, 177; stomach-UPB, 1,166; duodenum-UPB, 110; pancreas-UPB, 356; colon-UPB, 280), and 740 out of 4,609 (16.1%) VSNs expressed were marked by two UPBs. A correlation matrix (Extended Data Fig. 2e ) was calculated on the basis of the numbers of single-UPB and dual-UPB marked VSNs across the seven examined visceral organs and plotted using the R pheatmap package. Correlation analysis indicates that the two UPBs were predominantly from physically adjacent rather than random organs (Extended Data Fig. 2e ). This observation was consistent with retrograde tracing results (Extended Data Fig. 1b, v, w ), demonstrating that (1) UPB leakage was minimal during sample preparation; and (2) dual-labelled VSNs innervate regions close to both organs. Single UPB-labelled VSNs were designated as organ-specific VSNs, VSNs dual labelled with UPB-stomach and UPB-oesophagus were designated as oesophageal-sphincter projecting neurons and VSNs dual labelled with UPB-stomach and UPB-duodenum were designated as pyloric-sphincter projecting neurons. To examine VSNs that innervate different physiological systems, VSNs singly marked by oesophagus, stomach, duodenum and colon-UPBs were combined as gut VSNs (Fig. 1c, d ). A 3D UMAP plot (Fig. 1c , bottom) was generated using a previously published method 60 ( https://github.com/Dragonmasterx87/Interactive-3D-Plotting-in-Seurat-3.0.0 ). DEGs for organ-specific VSNs, or between thoracic (combination of oesophagus, heart and lung-UPB labelled) and abdominal (combination of stomach, duodenum, colon and pancreas-UPB labelled) VSNs were identified using the Wilcoxon rank-sum test (Extended Data Fig. 3e ). Mouse transcription factors were identified (Fig. 1d , Extended Data Fig. 3g ) by comparing DEGs with the AnimalTFDB 3.0 mouse database 61 ( http://bioinfo.life.hust.edu.cn/AnimalTFDB/#!/ ). For example, Pou4f1 , which is essential for DRG neuron specification 32 , is preferentially expressed in lung VSNs (Extended Data Fig. 3g ). Regulatory networks in UPB-labelled VSNs (Extended Data Fig. 3h ) showing the expression of Pou4f1 downstream regulators including Runx1 and Isl2 were predicted by IPA analysis (Qiagen) using the Upstream Analysis (Upstream Regulators) module. For cell–cell interaction analysis between organ-innervating VSNs and various organ cell types (Extended Data Fig. 3i ), UPB + VSNs labelled from the indicated organ were first extracted from the Projection-seq data and then integrated with a previously published scRNA-seq dataset containing various cell types from that corresponding organ using the R package Seurat v.3 55 . The following datasets were used: heart 62 (CM, cardiomyocyte; EDC, endothelial cell; EP, epicardial cell; FB, fibroblast), lung 63 (ATI, alveolar epithelial type I cell; ATII, alveolar epithelial type II cell; B, B cell; C&S, ciliated and secretory cell; DC, dendritic cell; EDC, endothelial cell; FB, fibroblast; MO, monocyte; Mac, macrophage; NK, natural killer cell; Neutro, neurophil; Peri, pericyte; T, T cell), colon 64 (Endo, endothelial cell; Immu, immune cell; EN, enteric neuron; Glia, glial cell; Entero, enteroendocrine cell; Mus/Fb, muscle cell and fibroblast), duodenum 65 (duodenal enteric neuron subtypes: EXMN, excitatory motor neuron; INMN, inhibitory motor neuron; IN, inter neuron; IPAN, intrinsic primary afferent neuron) and pancreas 60 (ISL, islet cell; ACI, acinar cell; DUCT, duct cell; MES, mesenchymal cell; IMVS, immune and vascular cell).

Cell–cell interactions between organ-innervating UPB +

VSNs and indicated organ cell types (or duodenal enteric neuron subtypes) were then analysed using the CellPhoneDB 66 ( https://github.com/Teichlab/cellphonedb , v.2.0.0) Python package. GO pathway analyses of DEGs in heart, lung, gut and pancreas VSNs, and along the tissue trajectory (Extended Data Fig. 8s ), were performed using the Gene Ontology Resource GO Enrichment Analysis tool 67 – 69 ( http://geneontology.org ).

Calculation of trajectory score

Distances from the nodose ganglia to the beginning and the end of various organs were measured and normalized to the body length from the neck to the rectum (Fig. 1e–g ). The mean distances were calculated as Position organ (lung: 0.241 ± 0.005; heart: 0.280 ± 0.005; oesophagus: 0.228 ± 0.005; stomach: 0.469 ± 0.006; duodenum: 0.567 ± 0.006; transverse colon: 0.589 ± 0.005; pancreas: 0.545 ± 0.007; n = 4 mice). The organ position score for VSN clusters (Fig. 1e ), indicating their target preference along the body’s rostral–caudal axis, was calculated as the weighted average of Position organ using the percentage of organ-specific VSNs in the cluster ( P organ-cluster = Number of UPB organ + VSNs in the target cluster/Number of all UPB organ + VSNs) as the weight value for each organ, expressed as Σ( P organ-cluster × Position organ )/Σ P organ-cluster . The organ trajectory score of an organ-specific VSN (Fig. 1f, g ) was measured as its projection position along the organ trajectory on the UMAP plot (as shown in Fig. 1e ). Tissue layer trajectory was identified using Slingshot 70 , and DEGs along this trajectory were discovered using tradeSeq 71 .

Tissue layer trajectory scores of DEG +

VSNs were measured as their relative positions along the tissue layer trajectory by Slingslot using the pseudotime function. On the basis of the inner–outer position, each tissue type was given an index score (mucosa or inner epithelium, 0; muscle, 1; connective tissue, 2), and the tissue layer index for DEG + VSNs (Fig. 2a–c ) in each organ was calculated based on the percentage of DEG + VSN endings in the target tissue layer and the tissue layer index across multiple tissue layers as Σ((Number of DEG + VSN endings in the target tissue layer/Number of DEG + VSN endings in all tissue layers) × Index tissue ) (Fig. 2a–c ).

AAV infection of vagal ganglia

Vagal ganglia injection was performed as previously described 11 , 12 . In brief, mice were anaesthetized with 1–2% isoflurane and maintained on a heating pad. Both left and right vagal ganglia were surgically exposed. A virus mix containing a 1:1 dilution of AAV9-FLEX-tdTomato and AAV5-CAG-GFP with 0.05% (w/v) Fast Green FCF was injected (160 nl for each side, 20 nl per second) using a Nanoject III injector (Drummond). Mice were euthanized four weeks after surgery for tissue collection (see ‘Histology and immunochemistry’). RNAscope HiPlex assay RNAscope HiPlex assays were performed following the manufacturer’s protocol (Advanced Cell Diagnostics). Vagal ganglia were acutely dissected and freshly frozen in cryo-embedding medium (OCT). Cryosections (10 μm for vCatFISH; 14 μm for others) were cut using a cryostat (Thermo Fisher Scientific), mounted onto Superfrost Plus slides (Thermo Fisher Scientific), and stored at −80 °C until use. Slides were immediately immersed into fresh 4% paraformaldehyde (PFA) in RNase-free PBS for 60 min at room temperature, followed by dehydration with 50%, 70% and 100% ethanol. Samples were then digested with protease IV for 30 min at room temperature. After hybridization with designed probes (Supplementary Table 1 ) for 2 h at 40 °C, the sections were treated with HiPlex Amp 1–3, and then HiPlex Fluoro T1–T3, before counterstaining and mounting. The slides were imaged using a Leica SP8 confocal microscope equipped with a motorized stage, a PMT detector, a HyD SP detector, four laser lines (405 nm, 488 nm, 552 nm and 638 nm) and a 20× objective (HC PL APO 20×/0.75 CS2). Imaging loci were exactly recorded in the LAS X software in the first group for image registration, and then applied for imaging the same sections for the following groups. After each group, the fluorophores were cleaved by 10% cleaving solution (ACD, 324130). The sections were then hybridized with HiPlex Fluoro T4–T6 in group 2, HiPlex Fluoro T7–T9 in group 3 and HiPlex Fluoro T10–T12 in group 4. In total, 12 genes were detected on a single section. The following changes were made for vCatFISH: (1) slides were not mounted after HiPlex Fluoro and instead imaged in 4× SSC with a 16× immersion objective (HC FLUOTAR L 16×/0.8 IMM motCORR VISIR); (2) after the first-round imaging (4 groups, 12 probes), probes were removed using the HiPlexUp reagent, and sections were hybridized with another 12 probes (2 h, 40 °C) for the second-round analysis. In total, 22 genes were analysed (see ‘vCatFISH analysis’). Histology and immunochemistry Mice were anaesthetized with urethane (1.5 g per kg) and transcardially perfused with 15 ml cold PBS (pH 7.4) containing 10 U ml −1 heparin (Sigma-Aldrich, H4784), followed by 25 ml cold 4% PFA. Fast Green FCF injected organs (Extended Data Fig. 1e, g, j, s–u ) were dissected and imaged under a Stereo Microscope with a digital camera (AmScope). For all others, visceral organs and/or vagal ganglia were dissected, post-fixed in 4% PFA at 4 °C (overnight for visceral organs, 30 min for ganglia), and then kept in cold PBS at 4 °C. Brains were post-fixed (4% PFA, overnight, 4 °C), cryoprotected in 30% sucrose PBS solution for two days at 4 °C, frozen in OCT and then stored at −80 °C until cryosection. For brain samples, 40-μm cryosections were mounted onto Superfrost Plus slides. Cryosections were washed (3× PBS), permeabilized (0.1% Triton X-100, PBS), blocked (5% normal donkey serum, PBST (PBS, 0.05% Tween-20)) and incubated with primary antibodies (chicken anti-GFP, 1:1,000, Aves Labs; rabbit anti-RFP, 1:1,000, Rockland) diluted in blocking buffer for 2 h at room temperature. Then, the slides were washed (3× PBST), and incubated with fluorophore-conjugated secondary antibodies (Alexa Fluor 647-conjugated AffiniPure donkey anti-chicken IgY, 1:1,000; Alexa Fluor 594-conjugated AffiniPure donkey anti-rabbit IgG(H+L), 1:1,000, Jackson ImmunoResearch) diluted in blocking buffer for 2 h at room temperature. After incubation, the samples were washed (3× PBST), and mounted with Fluoromount-G with DAPI before imaging with the Leica SP8 confocal microscope. Visceral organs were cleared with the CUBIC method 72 and stained with the following protocol unless specifically mentioned. In brief, the dissected organ was immersed into 1/2-water-diluted reagent-1 (25 wt% urea, 25 wt% Quadrol, 15 wt% Triton X-100) with shaking at 37 °C for 3–6 h, followed by reagent-1 (R1) with shaking at 37 °C for 7 days. R1 was replaced fresh every two days. Next, the tissue was washed (3× PBS/0.01% NaN 3 ), blocked (2% normal donkey serum, 0.1% Triton X-100, PBS/0.01% NaN 3 ) and incubated with primary antibodies (chicken anti-GFP, 1:200; rabbit anti-RFP, 1:200) in blocking buffer with shaking for 7 days at room temperature. Samples were then washed (0.1% Triton X-100, PBS/0.01% NaN 3 ) and incubated with fluorophore-conjugated secondary antibodies diluted in blocking buffer with shaking for five days at room temperature. After antibody incubation, the samples were washed and immersed in 1/2-PBS-diluted reagent-2 (25 wt% urea, 50 wt% sucrose, 10 wt% triethanolamine) overnight at room temperature, and then reagent-2 (R2) at 37 °C for 2 days. The samples were finally immersed in oil for at least 1 h and flattened to approximately 500 μm in a custom-built imaging chamber and imaged using the Leica SP8 confocal microscope as described above, with a 10× objective (HC PL APO 10×/0.40 CS2, working distance: 2.1 mm) or a 40× objective (HC PL FLUOTAR L 40×/0.60 CORR, working distance: 3.3 mm). Some heart samples were similarly processed but not flattened and imaged with a LaVision Vltramicroscope II light-sheet microscope at the CNNR Imaging Core at Yale University (Extended Data Fig. 5a ) or the Leica SP8 confocal microscope with a 16× immersion objective (HC FLUOTAR L 16×/0.8 IMM motCORR VISIR, working distance: 8 mm). Heart slices (1 mm thickness) and gastrointestinal organs to determine viral spread (Extended Data Fig. 1f, v, w ) were imaged under a Leica M205FCA Fluorescent Stereo Microscope with a CoolLED pE300 white illumination, GFP and DSR Filter sets, 1× objective (Plan M-series) and a Leica DFC7000 T camera. The pancreas was excluded from anatomical analysis owing to multiple technical challenges (keeping its original structure and distinguishing various tissue layers after clearing, antibody penetration and imaging efficiency). A modified CUBIC protocol was used for clearing vagal ganglia: samples were cleared in R1 for one day, incubated in primary and secondary antibodies (chicken anti-GFP, 1:1,000; rabbit anti-RFP, 1:1,000; Alexa Fluor 647-conjugated AffiniPure donkey anti-chicken IgY, 1:1,000; Alexa Fluor 594-conjugated AffiniPure donkey anti-rabbit IgG(H+L), 1:1,000) overnight, respectively, and treated with R2 for one day. Cleared ganglia were imaged with the Leica SP8 confocal microscope. vCatFISH analysis Surgery, stimulus delivery and imaging Gpr65 tdT -GCaMP6s mice were continuously anaesthetized (1–2% isoflurane/oxygen) during the experiment. A tracheal tube was inserted for air injection. The upper oesophagus and pyloric sphincter were cannulated and flushed with saline multiple times to remove residual food particles. The duodenum (around 0.5 cm below the pyloric sphincter) was cannulated with a bundle of six PE-10 tubing for separate delivery of saline, water and glucose (1 M), Ensure (Ensure Original Vanilla Nutrition Shake), 10× PBS and 150 mM HCl (pH = 0.84). Left vagal ganglia were exposed and immobilized on a stable platform 12 . During calcium imaging, a series of stimuli were delivered to the same mice in the following sequence: (1) lung inflation for 20 s with 600 ml min −1 flow (oxygen through the tracheal tube), twice with a 2-min interval; (2) small intestine stretch via fast injection of 600 μl saline through the duodenal cannula; stomach stretch with 100 μl, 300 μl and 600 μl saline through the oesophageal cannula for 30 s (duration precisely controlled by closing or opening of the pyloric cannula); (4) small intestine infusion with 100 μl saline, water, 1 M glucose, Ensure, 10× PBS and 150 mM HCl in sequence with a 3-min interval between infusions. GCamp6s fluorescence was measured from two focal planes 15 μm apart using a two-photon microscope (920 nm excitation, Leica TCS SP8, Mai Tai laser from Spectra-Physics, HyD SP detector). The imaging frequency for each plane was 1.72 s per frame. At the end, electrical stimulation was applied to the vagus nerve and a z -stack of the ganglia was taken for both GCaMP and tdTomato (552 nm, one photon excitation) signals for cell registration. RNAscope and registration After GCaMP imaging, vagal ganglia were immediately embedded in OCT, frozen in situ, and cut into 10 um cryosections. RNAscope HiPlex assays for the following 22 genes were performed (Fig. 3b , Extended Data Fig. 9 ): Trpa1 (R1T1), Runx3 (R1T2), Uts2b (R1T3), Gabra1 (R1T4), Slit2 (R1T5), Kcng1 (R1T6), Piezo2 (R1T7), Ddc (R1T8), Vip (R1T9), Trpv1 (R1T10), tdTomato (R1T11), Gpr65 (R1T12), Chodl (R2T1), Glp1r (R2T2), Grm5 (R2T3), Slc17a7 (R2T4), P2ry1 (R2T5), Tmc3 (R2T6), Car8 (R2T7), Nts (R2T8), Cckar (R2T9) and Calca (R2T11). The following criteria were used to categorize VSN subpopulations (Fig. 3b , Extended Data Fig. 9 ): A-VSNs: Runx3 + and/or Piezo2 + / P2ry1 + , Tmc3 − ; B-VSNs: Gabra1 + ; C-VSNs: multiple hits for Slit2 + , Piezo2 + , Ddc + , Tmc3 + , Trpv1 − , P2ry1 − . D-VSNs: Tmc3 + , Trpv1 − ; F-VSNs: Gpr65 + ; G-VSNs: Trpv1 + , multiple hits for Uts2b + , Vip + , Glp1r + , Cckar + ; H-VSNs: Trpv1 + , Tmc3 − , Trpa1 + ; I-VSNs: multiple hits for Tmc3 + , Car8 + , Cckar + , Piezo2 − . J -VSNs: Trpv1 + , Calca + and some P2ry1 + , Piezo2 − ; K-VSNs: Trpa1 + , Kcng1 + , Trpv1 + , Calca + ; L-VSNs: P2ry1 + , Trpv1 − . tdTomato + neurons that serve as geographic landmarks for cell registration were registered first between in vivo GCaMP images and RNAscope images (Extended Data Fig. 9e ). Coordinates of registered tdTomato + cells were then used to calculate the transformation matrix between in vivo 3D images and RNAscope sections. We reason that tdTomato + cells from a given RNAscope section should also be located in the same plane in the transformed in vivo image stack, therefore we used a plane correction script simulating a virtual plane with minimum total projection distances for tdTomato + cells from multiple RNAscope sections. The transformation matrix was then applied to the in vivo image stack (3D extension/plug-in in LASX software) to generate successive in vivo imaging planes that resemble RNAscope sections. tdTomato − neurons were then registered on the basis of their relative distance and depth to their neighbouring tdTomato + cells (Extended Data Fig. 9e–g ). Registration for each sample was performed by at least two people independently. In total, 57.5% (349/607, 6 mice) of responsive VSNs were unambiguously registered, which is comparable to a similar approach developed in the trigeminal ganglion 73 . Cells that were not successfully registered were removed from further analysis.

Analysis of neural activity

Regions of interest were manually extracted from GCaMP images. The stimulus induction frame was set as 0 unless specifically mentioned. Baseline signal ( F ) was defined as the average GCaMP6s fluorescence over a 10-frame period (17.2 s, frame −20 to frame −10) before stimulus induction and neuronal activity was calculated as Δ F / F . Cells were coded as responsive to a given stimulus if the maximum GCaMP6s fluorescence was more than 100% above baseline during stimulus period (lung inflation and stomach stretch: between stretch on and off; intestine stretch: within 40 frames/68.8 s after injection; intestine infusion: within 90 frames/154.8 s after infusion). Peak response was identified as the maximum Δ F / F within the stimulus period. To compare adaptation rates in lung stretch-sensitive VSNs, GCaMP6s traces were aligned at activation frame (Extended Data Fig. 11c , arrow, set as 0), and the activation duration was calculated as the number of frames between the prior and post peak frames at which VSN activity (Δ F / F ) reached 10% of the peak response (Extended Data Fig. 11d ). For intestine stretch, activation frame was defined as the prior peak frame at which VSN activity (Δ F / F ) reached 10% of peak response. Projection-seq-guided anterograde tracing in visceral organs The percentage of gene + VSNs targeted in various Cre lines in individual clusters (Extended Data Fig. 4c ) was calculated as the number of gene + VSNs normalized by the number of all VSNs in each cluster using the control scRNA-seq dataset. Fold enrichment of gene + VSNs in each cluster (Extended Data Fig. 4d ) was calculated as the percentage of gene + VSNs in the target cluster normalized by the overall percentage of gene + VSNs in all clusters.

Identification of enriched clusters and DEGs

Primary oesophagus, stomach, duodenum, colon, heart and lung VSN clusters were identified as clusters containing more than 4% of corresponding UPB single-labelled VSNs (Extended Data Figs. 5c , 6 b, 7 c, 8a, e, k, r ), with the following additions: G2, G5-duodenum VSNs (both 3.64%) and I4-heart VSNs (3.95%). VSN clusters enriched for stomach regions 4, 6–8 (Extended Data Fig. 7e, h ) were defined using the following criteria: (1) the cluster contains more than 4% of dual-UPB (region 4: oesophagus/stomach; region 6: stomach/pancreas; region 7: stomach/colon;region 8: stomach/duodenum, respectively) labelled VSNs; and (2) in the cluster, the percentage of dual-UPB labelled VSNs is at least 5% higher than both the percentage of stomach-UPB labelled VSNs and the percentage of the other UPB (oesophagus, pancreas, colon, duodenum, respectively) labelled VSNs. VSN clusters enriched for stomach region 5 (fundus, S-only, Extended Data Fig. 7h ) were defined as: (1) the cluster contains more than 4% of stomach-UPB labelled VSNs; and (2) in the cluster, the percentage of stomach-UPB labelled VSNs is at least 5% higher than the percentages of stomach/oesophagus, stomach/pancreas, stomach/colon and stomach/duodenum dual-UPB labelled VSNs. Among the eight primary duodenum VSN clusters (F1, F4, G2, G5, H2, H3, I5 and J3), three (F1, H2, J3) were more enriched in stomach/duodenum dual-labelled VSNs, suggesting that they preferentially project to the pyloric sphincter over the duodenum. Therefore, the other five duodenum VSN clusters were focused on for DEG analysis (Extended Data Fig. 8a ). Fractions of DEGs (Extended Data Figs. 5 e, 6 c, 7 f, 8b, f, l ) were calculated as the number of DEG + VSNs in the indicated cluster divided by the number of DEG + VSNs in all enriched clusters indicated in the panel. Percentages of gene + VSNs in identified enriched clusters (Extended Data Figs. 5 d, 6 d, 8h ) were calculated as the number of gene + VSNs divided by the total number of VSNs in that cluster.

Histology analysis

Four types of gut VSN ending types, three types of heart VSN ending types and five types of lung VSN ending types were classified on the basis of their morphologies and locations using whole-mount preparations in Vglut2 tdT mice (Fig. 2d , Extended Data Figs. 4 b, 5 b, 6a ; n = 4–7). For quantitative analyses of MEs and IMAs in the indicated gastrointestinal regions; alveoli, longitudinal and patch-terminals in the lung; and varicose endings and IMAs in the heart (Extended Data Figs. 5 f, 6 l, 7 b, m, 8d, j, n ), the area covered by each sensory ending type was measured using the Leica Application Suite X software and divided by the total area of each sample to derive the innervation intensity. For IGLEs in the gastrointestinal tract, neuroepithelial body (NEB) endings in the lung and flower-spray endings in the heart, the number of terminal clusters was counted (Extended Data Figs. 5 f, 6l ). The number of IGLEs was divided by the total area of each sample to derive the innervation intensity. Normalized innervation intensity for gastrointestinal endings (Extended Data Figs. 7 b, m, 8d, j, n ) was calculated as the innervation intensity of the indicated sensory ending type formed in indicated Cre tdT mice divided by the innervation intensity of the indicated sensory ending type formed in Vglut2 tdT mice. Fold changes for various sensory ending types in indicated stomach regions (Extended Data Fig. 7d, j ) were calculated as the innervation intensity of the corresponding sensory ending type over the indicated stomach region divided by the innervation intensity of the corresponding sensory ending type over the entire stomach in Vglut2 tdT mice ( n = 4). Annotation of VSN clusters Heart VSNs were predominantly distributed in four clusters (Extended Data Fig. 5 ), with the following DEGs: D1 ( Piezo2 ), H4 ( Drd2 ) and I3/I4 ( Agtr1a ). In Piezo2 tdT mice, most cardiac afferents were varicose surface endings. In Drd2 tdT mice, cardiac afferents densely innervated myocardium with both branched and parallel IMAs. Agtr1a + VSNs predominantly formed flower-spray endings in the heart and the aortic arch. Our results thus reveal the identity of various VSN cardiac ending types (Extended Data Fig. 5h ). Although five types of lung VSN endings were characterized, only two large groups of primary lung VSN clusters were revealed (Extended Data Fig. 6 ). A3-VSN fibres labelled in Agtr1a tdT , Vglut1 tdT , Pvalb tdT and P2ry1 tdT mice travelled along segmental bronchi and terminated at airway bifurcations wrapping around NEBs. Of note, VSN fibres in Agtr1a tdT and P2ry1 tdT mice also formed similar endings wrapping around taste buds in the larynx 17 and the upper oesophagus (Extended Data Fig. 6h ), suggesting that VSN bud endings in different organs are likely to have similar genetic signatures. K1–3/L2-VSNs labelled in Twist2 tdT and P2ry1 tdT mice formed alveoli endings. The other three VSN lung endings that were not effectively marked by Projection-seq were all on bronchial airways, consistent with our observation that lumen-delivered virus did not effectively cover bronchial walls (Extended Data Fig. 1i ). Npy2r + afferents formed dense longitudinal endings wrapping around the bronchioles; both Npy2r + and P2ry1 + VSNs formed patchy endings on segmental bronchi or around bronchial bifurcations; Piezo2 + pulmonary afferents mainly formed two ending types: (1) bud ending wrapping around NEBs (Extended Data Fig. 6g ) as Piezo2 is partially in A3-VSNs; and (2) near the bronchial bifurcation with branch endings (Fig. 3e , left), largely from the C3 cluster (1.5% of lung-UPB-marked VSNs). Annotations for lung VSN types are summarized in Extended Data Fig. 6m ). Most stomach VSNs were distributed in nine clusters (Extended Data Fig. 7c ). We divided the stomach into five regions based on their proximation to other UPB-targeted organs (Extended Data Fig. 7a, b ) and took advantage of dual-UPB labelled VSNs that innervate regions close to both organs to decode each region individually (region 4: stomach/oesophagus-UPB dual-labelled VSNs; region 6: stomach/pancreas-UPB dual-labelled VSNs; region 7: stomach/colon-UPB dual-labelled VSNs; region 8: stomach/duodenum-UPB dual-labelled VSNs). Compared to other stomach regions, pIMAs and cIMAs but no other endings were enriched around the oesophageal sphincter (Extended Data Fig. 7d , top). Accordingly, three clusters (J1, J3 and I1) were enriched in stomach/oesophagus-UPB dual-labelled VSNs (Extended Data Fig. 7e, f , top). We then identified DEGs ( P2ry1 for J1/J3 and Calb2 for I1) and examined IMAs around pyloric sphincters in P2ry1 tdT and Calb2 tdT mice. Most (75.0%) of the P2ry1 + IMAs around the oesophageal sphincter were pIMAs and all Calb2 + IMAs were cIMAs (Extended Data Fig. 7g ), suggesting that J1/J3-VSNs form pIMAs and I1-VSNs form cIMAs. In additional to J1- and I1-VSNs, three clusters were enriched in S/D-VSNs with the following DEGs: F1 ( Sst , Gpr65 ), H2 ( Vip ) and C4 ( Glp1r , Piezo2 ), whereas all afferent types in Vglut2 tdT mice were enriched in region 8 around the pyloric sphincter. Both Sst + and Gpr65 + VSNs formed MEs on the stomach, as reported 12 , 16 , with Sst + endings mainly in the antrum and Gpr65 + endings more evenly distributed across the stomach, indicating that F1-VSNs form MEs. Both Glp1r + and Piezo2 + VSNs predominantly formed IGLEs, suggesting that C4-VSNs form stomach IGLEs. Vip + VSNs exhibited pIMA morphology in region 8, suggesting that H2-VSNs also form pIMAs (Extended Data Fig. 7d–g ). We also used anatomical tracing results to facilitate the annotation of other stomach VSN clusters. We first determined the relative innervation intensity (referred as anatomical fold change or AF afferent type-stomach region ) for each ending type (ME, pIMA, cIMA and IGLE) in each stomach region (4, 6, 7 and 8), calculated as the innervation intensity of corresponding ending types in various stomach regions normalized by the innervation intensity of the same ending type in the entire stomach in Vglut2 tdT mice (Extended Data Fig. 7j ). According to anatomical results (Extended Data Figs. 7 , 8 ), F- and G-VSNs formed mucosal endings (MEs); H2-, H4-, J1- and J3-VSNs formed pIMAs; I1-VSNs formed cIMAs; and C4-VSNs formed IGLEs. We then performed sensory ending type simulation for seven VSN clusters (I2, I4, I5, I6, I7, J2 and J4) in the stomach (Extended Data Fig. 7i ) using a MATLAB script, assuming that each could form one of the four ending types (ME, pIMA, cIMA and IGLE) independently. Therefore, in total 4 7 = 16,384 possibilities were tested. For each possibility, the relative innervation intensity derived from Projection-seq data (referred as Projection-seq fold change or SF afferent type-stomach region ) of each ending type (ME, pIMA, cIMA and IGLE) in each stomach region (4, 6, 7 and 8) was calculated as the percentage of corresponding dual-labelled VSNs forming this ending type normalized to the percentage of stomach-UPB single-labelled VSNs forming this ending type, expressed as Σ(percentage of dual-labelled VSNs for this stomach region in all VSN clusters that form the corresponding ending type)/Σ(percentage of stomach-UPB single-labelled VSNs in all VSN clusters that form the corresponding ending type). We then calculated the total variance between anatomically and Projection-seq-derived innervation intensities in all stomach regions across all ending types, expressed as Σ( n = region 4, 6, 7, 8) (SF ME- n − AF ME- n ) 2 + Σ( n = region 4, 6, 7, 8) (SF pIMA- n − AF pIMA- n ) 2 + Σ( n = region 4, 6, 7, 8) (SF cIMA- n − AF cIMA- n ) 2 + Σ( n = region 4, 6, 7, 8) (SF IGLE- n − AF IGLE- n ) 2 . The trial with the lowest variance was defined as the best fit. AF afferent type-stomach region and SF afferent type-stomach region for this condition were plotted together (Extended Data Fig. 7j ). Simulation indicated that J2, J4 and I7 account for cIMAs, and I2 and I4–I6 for IGLEs. This prediction was further supported by our data and previous findings 16 that both Agtr1a + and Oxtr + VSNs, mainly I-VSNs, formed IGLEs close to the stomach antrum (Extended Data Fig. 7k, l ). Five clusters enriched for duodenum VSNs over stomach/duodenum dual-labelled VSNs were characterized with the following DEGs: Gpr65 (F4), Vip (H3/G2/G5), Glp1r (G2/G5), and Agtr1a (I5) (Extended Data Fig. 8a ). Gpr65 + and Vip + VSNs predominantly formed indistinguishable MEs arborizing intestinal villi in the duodenum (the density of Vip + MEs was much lower), suggesting that both F- and G-VSN clusters form MEs. By contrast, Agtr1a + neurons representing I5-VSNs primarily formed IGLEs in the duodenum (Extended Data Fig. 8b–d ). Similarly, of the five primary colon VSN clusters, Agtr1a is a DEG for I-VSNs and Agtr1a + VSNs mainly formed IGLEs in the colon (Extended Data Fig. 8e–g ). Trpv1 is highly expressed in both I- and H4- but only partially in F3-VSNs, and the percentage of IGLE, IMA and ME endings formed by Trpv1 + VSNs correlated well with the percentage of Trpv1 + VSNs in I, H4 and F3 clusters (Extended Data Fig. 8h–j ). Extensive Trpv1 + IMAs were observed in the oesophagus, and Trpv1 is highly expressed in J3 among the five primary oesophagus VSN clusters, consistent with the notion that J3 VSNs form IMAs in the stomach. Piezo2 + vagal oesophageal afferents—predominantly from C5–8 clusters—formed dense IMAs, suggesting that unlike C4-VSNs (stomach IGLEs), C5–C8 VSNs also form oesophageal IMAs. Probably owing to low infection efficiency, no apparent oesophageal ME or IGLE clusters were revealed. Notably, among all Cre lines examined, afferents in NTS tdT mice preferentially formed oesophageal IGLEs (Extended Data Fig. 8n ) and were later used to study the central projections of oesophageal IGLEs. Nevertheless, oesophageal MEs were predominantly formed by Trpv1 + and Gpr65 + VSNs, suggesting they might have similar genetic signatures to MEs in other gastrointestinal organs (Extended Data Fig. 8n, p ). Annotations for gut VSN types are summarized in Extended Data Fig. 8q .

Correlation analysis

Correlations among various VSN characteristics, including 7 visceral organs (lung, heart, oesophagus, stomach, pancreas, duodenum and colon), 4 tissue layer types (epithelium, specialized epithelial cells, muscle and elastic connective tissue), 6 VSN ending types (epithelial, budding, pIMA, cIMA, plate of puncta and varicose), 11 VSN subpopulations and 4 response patterns (mechanical sustained, mechanical transient, polymodal and chemical), were plotted in Fig. 3f . Correlation indexes between pairs of VSN characteristics (Fig. 3g ) were calculated on the basis of the number and pattern of connections between variables in the two characteristics, with two assumptions: (1) correlation is negatively related to the total number of connections: two characteristics are perfectly correlated (correlation index = 1) if every variable of one characteristic is only connected with one variable of the other characteristic (meaning a one-to-one correlation with the minimum number of possible connections), and two characteristics are completely uncorrelated (correlation index = 0) if every variable of one characteristic is connected to every variable of the other characteristic (meaning the maximum number of possible connections); and (2) with the same number of connections, correlation is stronger if the connections are more evenly distributed, meaning a smaller statistical variance of the number of connections among all variables within both characteristics. We first counted the number of variables within each characteristic (variables with no connections were removed from the analysis) and the number of connections between each pair of characteristics. We then calculated the normalized number of connections ( C ), expressed as (number of connections − minimum number of possible connections)/(maximum number of possible connections − minimum number of possible connections), for each pair of characteristics. We next calculated the statistical variance of connections per variable for both characteristics (v1 and v2). As both variances contribute to the correlation equally, the correlation index was finally calculated as (1 − C)/((v1 + 1) × (v2 + 1)). Projection-seq-guided retrograde tracing in the brainstem Mouse line and organ combinations were selected on the basis of Projection-seq anterograde tracing results: lung alveoli ending (Twist2-ires-Cre, lung), lung NEB ending (Agtr1a-Cre, lung), oesophageal IGLE (Nts-Cre, oesophagus), oesophageal IMA (Piezo2-GFP-ires-Cre, oesophagus), heart IMA (Drd2-Cre, heart), stomach IGLE (Glp1r-ires-Cre, stomach), stomach mucosal ending (Gpr65-ires-Cre, stomach), stomach IMA (Vip-ires-Cre, stomach), duodenal IGLE (Agtr1a-Cre, duodenum), duodenal mucosal ending (Gpr65-ires-Cre, duodenum) and colon IGLE (Agtr1a-Cre, colon). Organ injection of AAVrgs is described in ‘AAV infection of visceral organs’. Brain processing and imaging is described in ‘Histology and immunochemistry’.

Quantitative analyses of vagal central projections in the brainstem

Area innervated by vagal afferents retrogradely labelled from various visceral organs at different Bregma levels (−7.20, −7.32, −7.48, −7.56, −7.76, −7.92, −8.0 mm) were measured using Fiji (ImageJ) (Extended Data Fig. 13c ). The percentage innervation of indicated brainstem sub-nucleus at certain Bregma level (Fig. 4c , Extended Data Fig. 13d, e ) was calculated with the following steps: (1) the fluorescence intensity in each sub-nucleus (FI Bregma-subnucleus ) was calculated as average fluorescence in the sub-nucleus minus background fluorescence measured in a region in the sub-nucleus with no fluorescence-labelled vagal fibres; (2) the area of each sub-nucleus (A Bregma-subnucleus ) was measured; (3) the total fluorescence (TF Bregma ) was calculated as Σ(Bregma, subnucleus) (FI Bregma-subnucleus × A Bregma-subnucleus ); and (4) the percentage innervation (PI) of a sub-nucleus at a certain Bregma was calculated as FI Bregma-subnucleus × A Bregma-subnucleus /TF × 100. The correlation variance between VSNs labelled from various visceral organs (Extended Data Fig. 13f ) was calculated as Σ(all subnuclei, all Bregma level) (PI pathway 1 − PI pathway 2 ) 2 and a phylogenetic tree was generated based on the correlation variance matrix using the seqlinkage function in MATLAB.

Statistics and reproducibility

All statistical analyses were performed using GraphPad Prism 8. All data are reported as mean ± s.e.m. unless specifically mentioned. Significance for pair comparisons was determined and P values were reported using a two-tailed Student’s t test. Significance for multiple comparisons was first determined using a one-way ANOVA and adjusted P value was then reported using Tukey’s multiple comparisons test. The sample sizes represent the number of mice used for experiments and data analysis, which is determined according to the consistence of the results. No significant inter-individual variability was observed in our results, suggesting that the sample sizes were sufficient to demonstrate the findings. The precise number of mice is reported in the figure legends and Methods. The study did not report contrasts between treatment and control groups, so randomization was not applicable. For sequencing and tracing experiments, there were no treatment and control groups, and no hypothesis was tested regarding molecular identities. Thus, blinding was not applicable. Representative images and experiments were repeated independently in multiple mice with similar results: Fig. 2d , Plate of puncta: stomach ( n = 4), heart ( n = 3); varicose: heart ( n = 5), aorta ( n = 6); pIMA: colon ( n = 3), heart ( n = 6); cIMA: stomach ( n = 4), oesophagus ( n = 3); epithelial: stomach ( n = 4), lung ( n = 5); bud: lung ( n = 5), oesophagus ( n = 3); Extended Data Fig. 1a , n = 3; Extended Data Fig. 1f , n = 2; Extended Data Fig. 1h, i , n = 3; Extended Data Fig. 1k–q , n = 3; Extended Data Fig. 1r , n = 2; Extended Data Fig. 1s, t , n = 3; Extended Data Fig. 1v , n = 2; Extended Data Fig. 1w , n = 4; Extended Data Fig. 2a , n = 5; Extended Data Fig. 2b , n = 2; Extended Data Fig. 2k, l , n = 4; Extended Data Fig. 2p , n = 4; Extended Data Fig. 3c , n = 4; Extended Data Fig. 4b , oesophagus ( n = 3), stomach ( n = 4), duodenum ( n = 3), colon ( n = 3), aorta ( n = 6); Extended Data Fig. 4e , n = 3; Extended Data Fig. 5a , n = 4; Extended Data Fig. 5g , Piezo2 ( n = 9), Drd2 ( n = 6), Agtr1a ( n = 3), Npr2r ( n = 2); Extended Data Fig. 6e , n = 4; Extended Data Fig. 6g–k , Npy2r ( n = 3), Trpv1 ( n = 5), P2ry1 ( n = 6), Agtr1a ( n = 4), Piezo2 ( n = 9), Vglut1 ( n = 4), Pvalb ( n = 3), Twist2 ( n = 3); Extended Data Fig. 7g, l , P2ry1 ( n = 4), Calb2 ( n = 3), Glp1r ( n = 5), Vip ( n = 5), Gpr65 ( n = 7), Sst ( n = 5), Agtr1a ( n = 12); Extended Data Fig. 8c , Glp1r ( n = 3), Agtr1a ( n = 4); Extended Data Fig. 8g , Trpv1 ( n = 4), Agtr1a ( n = 3); Extended Data Fig. 8o, p , Nts ( n = 5), Trpv1 ( n = 3), Agtr1a ( n = 4), Gpr65 ( n = 3); Extended Data Fig. 9a, e, f , n = 5; Extended Data Fig. 10a , n = 5; Extended Data Fig. 13a , lung ( n = 3), heart ( n = 3), oesophagus ( n = 2), stomach ( n = 3), duodenum ( n = 3), colon ( n = 3), pancreas ( n = 3); Extended Data Fig. 14 , n = 3 for all groups. Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Online content Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-022-04515-5.

Supplementary information Supplementary Figure 1 Original source image for electrophoresis, related to Extended Data Fig. 2b. Red box indicates how the gel was cropped for the final figure. Reporting Summary Supplementary Table 1 Detailed information for RNAscope probes.