Lymph nodes are innervated by a unique population of sensory neurons with immunomodulatory potential.

RESOURCE AVAILABILITY Lead Contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Ulrich H. von Andrian ( uva@hms.harvard.edu ).

Materials Availability

No new or unique reagents were generated in this study.

Data and Code Availability

The raw data and gene expression matrices for all scRNA-seq data have been deposited in the Gene Expression Omnibus ( https://www.ncbi.nlm.nih.gov/geo/ ) under the accession number GSE139600 for SMART-seq2 of LN-innervating or skin-innervating neurons, and accession number GSE139658 for Seq-Well of LN-resident cells. Raw data and metadata can also be downloaded and visualized on the Single Cell Portal: LN-resident cells, https://singlecell.broadinstitute.org/single_cell/study/SCP1186 . LN-innervating and skin-innervating neurons, https://singlecell.broadinstitute.org/single_cell/study/SCP1187 . Gene lists corresponding to differential expression tests in Figures 4E , 5A , 5B , S4F – S4H can be found in Table S1 . Cell type-identifying gene lists corresponding to Figures 6B , 6C , S5A – S5P can be found in Table S2 . All differentially expressed genes from optogenetic stimulation experiments presented in Figure 7C – 7H can be found in Table S3 . All R code for analysis available upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Mouse lines used in this study were all previously described and include Nav1.8 Cre (RRID:IMSR_EM:04582) ( Nassar et al., 2004 ), Rosa26 LSL-tdTomato (RRID:IMSR_JAX:007914), Bmx-CreER T2 (MGI:5513853) ( Ehling et al., 2013 ), Rosa26 LSL-DTA (RRID:IMSR_JAX:009669), Prox1-EGFP (MGI:4847348) ( Choi et al., 2011 ), Rosa26 LSL-ChR2-eYFP (RRID:IMSR_JAX:024109), Rosa26 LSL-eYFP (RRID:IMSR_JAX:007903) and ChAT BAC -eGFP (RRID:IMSR_JAX: 007902). All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols of Harvard Medical School. Unless indicated otherwise, adult mice of both sexes between 6–12 weeks of age were used for various experiments.

METHODS DETAILS Whole mount immunohistochemistry

Whole mount immunohistochemistry of LNs was performed using an iDISCO protocol with methanol pretreatment optimized for LNs ( Renier et al., 2014 ). Briefly, adult animals (6–12 weeks) were perfused with 25 mL of PBS (Hyclone) and 25 mL of 4% paraformaldehyde (PFA, Sigma) sequentially at room temperature (RT). Peripheral lymph nodes (PLNs), including popliteal and inguinal lymph nodes (popLNs and iLNs), were postfixed with 4% PFA for 1 hr at 4°C. For methanol pretreatment, fixed LNs were washed sequentially in 50% methanol (Fisher Scientific) (in PBS) for 1 hr, 100% methanol for 1 hr, 50% methanol for 1 hr, PBS for 1 hr twice, and PBS/0.2% Triton X-100 (VWR) for 1 hr twice at RT. LNs were then left in PBS/0.2% Triton X-100/20% DMSO (Sigma)/0.3 M glycine (BioRad) overnight at RT and blocked in PBS/0.2% Triton X-100/10% DMSO/6% donkey serum (Jackson Immunoresearch) or goat serum (Gibco)/anti-CD16/CD32 (Fc block) (Bio X cell) overnight at RT. LNs were subsequently washed in PBS/0.2% Tween-20 (Fisher Scientific)/10 μg/mL heparin (Sigma) (PTwH), for 1 hr twice at RT, before incubation with antibody mix in PTwH/5% DMSO/3% donkey or goat serum/Fc block 1:100 for 3 days at RT. LNs were extensively washed in PTwH for at least 6 times over the course of a day at RT. For unconjugated antibodies, LNs were further incubated with a secondary antibody mix including a panel of species-specific anti-IgG (H+L) Alexa Fluro 488, 546, 647 and 594-conjugated antibodies (Invitrogen or Jackson Immunoresearch) in PTwH/5% DMSO/3% donkey or goat serum/Fc block 1:100 for 3 more days at RT. LNs were washed in the same way as after primary antibody incubation for 1 day. Immunolabeled LNs following one round of antibody incubation for conjugated antibodies (or two for unconjugated antibodies) were then processed for clearing, which includes sequential incubation with 50% methanol for 1 hr, 100% methanol for 1 hr for three times and a mixture of 1-part benzyl alcohol (Sigma): 2-parts benzyl benzoate (Sigma) (BABB) overnight at RT. For tdTomato immunolabeling, goat anti-mCherry antibody (ACRIS) was preabsorbed against PLNs from tdTomato − animals overnight at RT prior to use. Whole mount immunohistochemistry of DRGs and the skin was performed as described previously ( Li et al., 2011 ). Briefly, DRGs inside vertebral column and the depilated hairy skin from PFA-perfused animals (6–12 weeks) were postfixed with 4% PFA for 1 hr or Zamboni fixative (Fisher Scientific) overnight, respectively at 4°C. Samples were washed every 30 min with PBS/0.3% Triton-100 (0.3% PBST) for 4–6 hr, then incubated with primary antibodies in antibody diluent (0.3% PBST/20% DMSO/5% donkey or goat serum) for 2–3 days at RT. Samples were then washed with 0.3% PBST every 30 min for 5–8 hr before incubation with secondary antibodies in antibody diluent for 2–3 days at RT. After extensive washes as described above, samples were dehydrated and cleared in 50% methanol for 1 hr, 100% methanol for 1 hr for three times and BABB overnight at RT. Cleared whole-mount tissues were imaged in BABB between two coverglasses using Olympus FV3000 confocal imaging system, except for those shown in Figures 7K and 7L , which were acquired on BioRad 2100MP system and those shown in Figures 3B , S3B and S3E , which were acquired on Zeiss Stereo Discovery V16. The antibodies used were: rabbit anti-CGRP (Immunostar, 24112, 1:500), chicken anti-GFP (Aves Labs, GFP-1020, 1:500), chicken anti-NF200 (Aves Labs, NFH, 1:500), rabbit anti-Tyrosine Hydroxylase (Millipore, AB152, 1:500), goat anti mCherry antibody (1:500, ACRIS AB0040-200), rabbit anti-βIII-Tubulin (Biolegend, 802001, 1:500), Alexa Fluor 647-conjugated rat anti-CD31 (Biolegend, 102416, 1:50), FITC–conjugated mouse anti-smooth muscle actin (aSMA) (Sigma, F3777-.2ML, 1:500), eFluor 660-conjugated mouse anti-smooth muscle actin (aSMA) (Thermo Fisher, 50-9760-82, 1:100), eFluor 660-conjugated rat anti-CD169 (Thermo Fisher, 50-5755-80, 1:50), Pacific Blue-conjugated rat anti-CD45 (Biolegend, 103126, 1:50), Alexa Fluor 488-conjugated rat anti-PNAd (Thermo Fisher, 53-6036-82, 1:50) Retrograde labeling of LN- and skin-innervating neurons To retrogradely label LN-innervating neurons, adult animals (6–12 weeks) were anesthetized by intraperitoneal injection of ketamine (Patterson Vet) (50 mg kg −1 ) and xylazine (Patterson Vet) (10 mg kg −1 ). The skin overlying the targeted iLN was shaved and depilated so that the LN underneath was visible percutaneously. A 5 mm incision was made directly on top of the iLN. The iLN was microdissected without perturbing afferent lymphatic vessels and surrounding blood vessels. 1 μl of Adeno-Associated Virus (AAV) (AAV2/1.CMV.HI.eGFP-Cre.WPRE.SV40, titer >=8E+12 vg/mL, Addgene) mixed with 0.5 μl of fast green (Sigma) was injected into the iLN of Rosa26 LSL-tdTomato/LSL-tdTomato animals using a pulled and trimmed glass pipette (FHC) which was connected to a 5 mL syringe through the aspiration assembly system (Sigma). The injection site was immediately rinsed with 2 mL of saline (Patterson Vet) to wash away any off-target virus before the incision was closed with sutures. Animals were sacrificed between 1 month and 6 months after injection for histology or scRNA-seq. To directly visualize the axonal projections of sensory neurons retrogradely labeled from the iLN, AAV carrying Cre-dependent tdTomato cassette (AAV2/1.CAG.Flex.tdTomato.WPRE.bGH, titer ≥10 13 vg/mL, Addgene) was injected into the iLN Nav1.8 Cre/+ animals as described above. For WGA-based retrograde labeling, 1 μl of WGA-AF488 (2 mg/mL in PBS, Invitrogen) was injected into the iLN of Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ animals as described before and the animals were processed for histology 4 days post injection. Retrograde labeling of skin-innervating neurons was described previously ( Kuehn et al., 2019 ). Briefly, following ketamine-xylazine mediated anesthesia, a single injection of 0.2 μl of various AAV2/1 viruses as described above and 0.1 μl of fast green was delivered using the injection device described above intradermally into the patch of depilated skin overlying the iLN of adult mice (6–12 weeks). Animals were sacrificed between 1 month and 6 months after injection for immunohistochemistry, RNAscope, or scRNA-seq.

Immunohistochemistry of tissue sections

Adult animals (6–12 weeks) were perfused with 25 mL of PBS and 25 mL of 4% PFA sequentially at RT. The intact vertebral column was postfixed overnight with 4% PFA at 4°C. DRGs were subsequently dissected and processed for cryosectioning. 14 μm serial cryosections were collected and processed for immunohistochemistry as described previously ( Li et al., 2011 ). In brief, sections were postfixed with 4% PFA for 10 min at RT. Following three washes with PBS, they were incubated with blocking buffer (PBS with 5% normal goat serum and 0.3% Triton-100) for 1 hr at RT. The sections were then incubated with Rabbit anti-TH (Millipore) in the same blocking buffer overnight at 4°C. The following day, sections were washed three times with wash buffer (PBS with 0.3% Triton-100) before incubation with goat Alexa Fluor 647-conjugated anti-rabbit (Invitrogen) for 1 hr at RT. Sections were then washed for three times with wash buffer before mounting in Fluoromount Aqueous Mounting Medium (Sigma). WGA-488 and tdTomato were visualized directly based on endogenous fluorescence. All the sections with tdTomato + cells were imaged at 20x using Olympus FV3000 confocal imaging system. Intravital two-photon microscopy Adult Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ animals (6–12 weeks) were given 1 μg of FITC-conjugated rat anti-CD169 antibody (BioRad) diluted in a total volume of 20 μl of PBS into the right footpad to label CD169 + subscapular macrophages inside the draining LN. Immediately after, the animals were prepared microsurgically for intravital two-photon microscopy as described before ( Mempel et al., 2004 ). Briefly, anesthesia during surgical preparation and imaging was achieved through the ketamine-xylazine method as described above. The right popLN was exposed and positioned with the cortex facing outwards with minimal perturbation to afferent lymphatic vessels and surrounding blood vessels, while the animal was immobilized onto a custom-built stage by its hip bone and the vertebral column. The imaging chamber was created around the exposed LN with high vacuum grease (VWR) on the side and a coverslip on top. A thermocouple (Omega) was placed next to the LN to monitor the local temperature, which was maintained between 36.5 and 37°C by a custom-built water bath heating system. Two-photon imaging was performed on a Bio-Rad Radiance 2100MP Confocal/Multiphoton microscopy system with two MaiTai Ti:sapphire lasers (Spectra-Physics) tuned to 800 nm and 900 nm for two photon excitation and second harmonic generation. Z-stacks of sensory innervation of the capsular/subcapsular space on the cortical side were acquired in 1 μm steps with a 20×, 0.95 numerical aperture objective (Olympus). Subcutaneous challenge by footpad injection of LPS or Pam3CSK4 LPS (Sigma) and Pam3CSK4 (Invivogen) were dissolved in water to a stock concentration of 5mg/ml and 2mg/ml respectively. 50μg of LPS or 40μg of Pam3CSK4 were injected into one footpad of adult Nav1.8 Cre/+ × Rosa26 LSL-tdTomato/LSL-tdTomato mice in a total volume of 20μl, while the contralateral footpad received equal volume of PBS. 3 days later, popLNs were harvested and processed for whole mount immunohistochemistry as described above. Manual cell sorting for scRNA-seq Adult mice with retrogradely-labeled LN- or skin-innervating neurons were sacrificed by CO 2 asphyxiation. T13 and L1 DRGs ipsilateral to the side of injection were quickly removed without nerves attached and checked for tdTomato labeling in cold HBSS (1X, no Ca 2+ or Mg 2+ ) (VWR) under Leica MZ10 F stereomicroscope with fluorescence. DRGs were immediately digested with 1 mL of papain solution (HBSS/10 mM HEPES (VWR)/500 μM EDTA (Westnet)/0.4 mg/mL L-Cysteine (Sigma)/1.5 mM CaCl 2 (Sigma)/20 unit/mL Papain (Worthington)) in a 37°C water bath for 10 min, with agitation every 2 min. DRGs were further digested with 1 mL of collagenase type II/dispase solution (HBSS/10 mM HEPES/4 mg/mL collagenase type II (Worthington)/5 mg/mL dispase (Thermo Fisher)) in a 37°C water bath for 30 min, with agitation every 10 min. Following centrifugation at 400 g for 4 min, digested DRGs were mechanically disrupted in 0.2 mL of complete L15 medium (L15 (Invitrogen)/10 mM HEPES/10% FBS (Germini)) by passing them first through a 1000 μL pipette tip up to 10 times, and then through a 200 μL pipette tip up to 5 times until the tissues were fully dissociated. To remove myelin/axonal debris, the cell suspension diluted in 1 mL of complete L15 medium was carefully layered on top of 5 mL of Percoll gradient (L15/10 mM HEPES/20% Percoll (GE Healthcare) and centrifuged at 400 g for 9 min. After removing the supernatant, cells were washed in 2 mL of L15/10 mM HEPES and centrifuged at 750 g for 3 min. Finally, cells were resuspended in 1 mL of cold sorting buffer (L15/10 mM HEPES/1 mg/mL BSA (VWR)/25 μg/mL DNase I (Roche)), and subjected to fluorescence-assisted single-cell picking as described previously ( Hempel et al., 2007 ). Briefly, the cell suspension diluted in 3 mL of sorting buffer was immediately transferred to a 35 mm petri dish (Scanning dish) with lane markings 6 mm apart and let sit on ice until most cells had settled to the bottom which normally takes 15–20 min. Rare fluorescent cells were readily identified under Leica MZ10 F stereomicroscope with fluorescence (transillumination off) by scanning the bottom of the dish lane by lane to maximize recovery and avoid rescanning. Zoom was set such that the field of view corresponded to the width of a single lane. To pick out fluorescent cells with minimal contamination from nonfluorescent cells, a pulled and trimmed micropipette (World Precision Instruments) was carefully lowered under transillumination into the sorting buffer until it was in the vicinity of the target cell. Simultaneous positive pressure was applied by mouth through the aspiration assembly system, as described above for retrograde labeling. Once the micropipette was in position, the target cell was gently aspirated into the micropipette through capillary action by transient release of positive pressure. The micropipette was quickly removed to prevent aspiration of unwanted cells or debris. The content of the micropipette, including the target cell, was expelled gently into a droplet of cold fresh sorting buffer on a different 35 mm petri dish (wash dish 1) under transillumination. Wash dish 1 was kept on ice while subsequent scans for fluorescent cells occurred. Once 16 or all the fluorescent cells, whichever comes first, were collected in wash dish 1, cells were washed two additional times by moving them one by one into a new droplet of sorting buffer on clean 35 mm petri dishes. Micropipettes were not reused for different cells to avoid cross contamination. After the final wash, each fluorescent cell was pipetted up and down the micropipette three times to remove unwanted contamination before being ejected into 10 μl of cold RLT (Qiagen) supplemented with 1% β-mercaptoethanol (Sigma) in a 96-well plate, and snap-frozen on dry ice and stored at −80°C. The entire manual sorting procedure was routinely completed in 1.5 hr. scRNA-seq of neurons using Smart-Seq2 Single-cell libraries were generated according to the SMART-seq2 protocol ( Picelli et al., 2014 ; Trombetta et al., 2014 ). Briefly, RNA from single-cell lysates was purified using AMPure RNA Clean Spri beads (Beckman Coulter) at a 2.2× volume ratio, and mixed with oligo-dT primer (SMART-seq2 3’ Oligo-dT Primer), dNTPs (NEB), and RNase inhibitor (Fisher Scientific) at 72°C for 3 minutes on a thermal cycler to anneal the 3’ primer to polyadenylated mRNA. Reverse transcription was carried out in a master mix of Maxima RNaseH-minus RT enzyme and buffer (Fisher Scientific), MgCl 2 (Sigma), Betaine (Sigma), RNase inhibitor, and a 5’ template switch oligonucleotide (SMART-seq2 5’ TSO) using the following protocol: 42°C for 90 minutes, followed by 10 cycles of 50°C for 2 minutes, 42°C for 2 minutes, and followed by inactivation at 70°C for 15 minutes. Whole transcriptome amplification was achieved by addition of KAPA HiFi HotStart ReadyMix (Kapa Biosystems) and IS PCR primer (ISPCR) to the reverse transcription product and amplification on a thermal cycler using the following protocol: 98°C for 3 minutes, followed by 21 cycles of 98°C for 15 seconds, 67°C for 20 seconds, 72°C for 6 minutes, followed by a final 5-minute extension at 72°C. Libraries were purified using AMPure XP SPRI beads at a volume ratio of 0.8x followed by 0.9x. Library size was assessed using a High-Sensitivity DNA chip (Agilent Bioanalyzer), confirming the expected size distribution of ~1000–2000 bp. Tagmentation reactions were carried out with the Nextera XT DNA Sample Preparation Kit (Illumina) using 250 pg of cDNA per single cell as input, with modified manufacturer’s instructions as described. Libraries were purified twice with AMPure XP SPRI beads at a volume ratio of 0.9x, size distribution assessed using a High Sensitivity DNA chip (Agilent Bioanalyzer) and Qubit High-Sensitivity DNA kit (Invitrogen). Libraries were pooled and sequenced using NextSeq500/550 High Output v2 kits (75 cycles, Illumina) using 30–30 paired end sequencing with 8-mer dual indexing. RNAscope The RNAscope Fluorescent Multiplex Assay (ACD Biosystems) was performed according to RNAscope Multiplex Fluorescent Reagent Kit v2 user manual for fresh-frozen tissue samples. Briefly, 14 μm fresh frozen sections from T13 and L1 DRGs with each side containing retrogradely-labeled tdTomato + LN- or skin- innervating neurons from the same animal were hybridized with RNAscope probes for Ptgir (487851), tdTomato (317041-C2), and Prokr2 (498431-C3) simultaneously. The probes were amplified and detected with TSA plus fluorescein, cyanine 3 and cyanine 5 (Perkin Elmer). The ACD 3-plex negative control probe was run in parallel on separate sections in each experiment to assess the background level and set the acquisition parameter. All sections with tdTomato + cells were imaged at 20x using an Olympus FV3000 confocal imaging system. The frequency of Ptgir + or Prokr2 + DRG neurons among the tdTomato + LN- or skin-innervating population was determined by considering all the tdTomato + cells that were recovered and uniquely-defined from a single animal. Tamoxifen treatment Tamoxifen (Sigma) was dissolved in corn oil (Sigma) at a concentration of 20 mg/mL by shaking overnight at 37°C, and stored at 4°C for the duration of the injections. For labeling arterial vessels with Bmx-CreER T2 , 0.5 mg of tamoxifen was delivered intraperitoneally to label arterial ECs into Bmx-CreER T2 ; Rosa26 eYFP/+ animals between 4–6 weeks of age daily for three consecutive days. Animals were analyzed between 1–3 weeks later.

6-OHDA treatment

For sympathetic denervation, the stock solution of 6-hydroxydopamine (6-OHDA) (Sigma) was prepared in water at 42 mg/mL and stored at −20°C. Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ animals from the same litter between the ages of 6–12 weeks were injected intraperitoneally with 6-OHDA (100 mg kg −1 ) or an equal volume of saline daily for 5 consecutive days. Animals were analyzed the following day.

Optogenetic stimulation of iLN-innervating sensory neurons

Age-matched adult Nav1.8 Cre/+ ; Rosa26 LSL-ChR2-eYFP/+ (ChR2+) or Nav1.8 Cre/+ ; Rosa26 LSL-eYFP/+ (ChR2−) animals (6–12 weeks) were deeply anesthetized (isoflurane, 1.5%–2%, Patterson Vet) maintained at normal body temperature with a water bath heating system (Baxter) during surgical preparation and photostimulation. The animals were surgically prepared for intravital optogenetic stimulation using a method that was adapted from a previously-described protocol for intravital microcopy of iLNs ( von Andrian, 1996 ). Briefly, the skin with the left iLN was flipped inside out following a small incision immediately left to the midline and glued onto a metal block to keep the medulla side of LN exposed. Care was taken not to overstretch the skin flap and damage lymphatic and blood vessels. The site of illumination, the branch point of the antero-posterior-running segment of the y-shaped superficial epigastric artery from where LN feeding arterioles emerged was located and exposed with microdissection without compromising the blood vessel integrity while the tissue was kept moist with normal saline. The stimulation chamber was then built around the iLN with vacuum grease on the side to keep solution from leaking, as well as a metal hairpin shaped tubing with hot water flowing inside on top of vacuum grease to maintain the tissue between 36.5 and 37°C. A thermocouple was placed next to the branch point to monitor the temperature at the tissue. An optic fiber (200 μm core, Thorlabs) coupled to a DPSS laser light source (473 nm, Shanghai Laser & Optics Century) was positioned for focal illumination directly on top of the branch point. The stimulation chamber was subsequently filled to the metal tubing with GenTeal Tears Lubricant Eye Gel (Alcon) to keep the tissue from drying out during stimulation. Pulsed light stimulation (5 ms pulses, 125 mW/mm 2 intensity, 20 Hz) was delivered to the targeted region for 3 hr under the control of a shutter system (Uniblitz). iLNs from both sides were immediately removed after light stimulation and kept in ice cold LN media (HBSS (Corning)/2% FBS/10 mM HEPES/2 mM CaCl 2 ) until subsequent processing.

LN Dissociation and Single Cell Isolation

LNs were kept on ice until processing, < 60 minutes between animal sacrifice and tissue digestion. To minimize biases introduced during tissue dissociation, a gentle dissociation protocol optimized for reliable isolation of both stromal and hemopoietic LN cells was used to generate a suspension of single cells from both the non-immune and the immune compartments ( Fletcher et al., 2011 ). LN media was aspirated, and each LN was placed in 1 mL of pre-warmed digestion media (0.8 mg/mL dispase, 0.2 mg/mL collagenase P (Roche), 50 μg/mL of DNase I in LN media). Using a pair of needle-nose forceps, the capsule of each LN was gently pierced, and the LN in digestion media were placed in a 37°C water bath for 20 minutes with no agitation. Next, LNs were gently agitated without touching the tissue, pelleted by gravity, and the 1 mL of digestion media supernatant was removed and placed in a collection tube on ice containing 10 mL of quenching buffer (PBS/5 mM EDTA/5% FBS). A fresh 1 mL of digestion buffer was added to each LN, and the LNs were placed back in the 37°C water bath for an additional 5 minutes. The LN was gently agitated and triturated using a 1000 μL pipette tip, solid capsular and stromal matter was allowed to settle to the bottom of the tube without centrifuging, and the supernatant digestion media was added to the same collection tube containing quenching buffer. 5-minute incubation periods in fresh digestion buffer and trituration with a 1000 μL pipette tip continued until LNs were completely digested, typically requiring 3–4 additional digestion steps. The cellular suspension in quenching buffer was filtered through a 100 μm filter, and washed with an additional 15 mL of quenching buffer. Single-cell suspensions were centrifuged at 300g for 3 minutes at 4°C, and counted using a hemocytometer and light microscope. We recovered an average of 4.00 +/− 0.53 million cells per LN, and observed no differences in cellularity by treatment group or animal genotype. We saved an aliquot of 60,000 cells from each sample in quenching media on ice as the unenriched sample, and centrifuged the remaining cells at 300 g for 3 minutes at 4°C. Next, using the Miltenyi CD3ε microbead kit and CD19 mouse microbead kit, all remaining LN cells were stained according to manufacturer instructions with the following modifications. First, single cells were stained with CD3ε biotin for 10 minutes on ice, washed once with MACS buffer (PBS/0.5% BSA (Sigma)/2 mM EDTA) and stained simultaneously with CD19 microbeads and biotin microbeads. Cells were isolated using LD columns (Miltenyi) according to manufacturer specifications and the flow-through was collected as the non-T and non-B enriched sample. Single cells from both enriched and unenriched samples were pelleted by centrifugation at 300g for 3 minutes at 4°C, and counted using a hemocytometer with trypan blue staining to estimate cell viability. Across 14 LNs, we recovered an average of 270,000 +/− 31,000 (mean +/− SEM) cells per lymph node following CD3ε and CD19 depletion with > 90% viability. For LN cellularity analysis, single-cell suspensions of the two iLNs from the same ChR2+ or ChR2− mouse (6–12 weeks) were prepared as above. The cells were then filtered through steel mesh and resuspended at the appropriate cell density in FACS buffer before being acquired on a BD Accuri ™ C6 Plus flow cytometer (BD Biosciences). LN scRNA-seq using Seq-Well Single cells from each lymph node prior to and post CD3ε and CD19 depletion were kept separate and diluted to 15,000 cells in 200 μL complete media (RPMI 1640/10% FBS). Seq-Well was performed as described with changes noted below ( Aicher et al., 2019 ; Gierahn et al., 2017 ). Briefly, a pre-functionalized PDMS array containing ~86,000 nanowells was loaded with uniquely-barcoded mRNA capture beads (ChemGenes) ( Macosko et al., 2015 ) and suspended in complete media for at least 20 minutes. 15,000 cells were deposited onto the top of each PDMS array and let settle by gravity into distinct wells. The array was gently washed with PBS, and sealed using a functionalized polycarbonate membrane with a pore size of 0.01 μm, which allows exchange of buffers without permitting mixing of cell materials between different wells. Seq-Well arrays were sealed in a dry 37°C oven for 40 minutes, and submerged in a lysis buffer containing 5 M guanidium thiocyanate (Sigma), 1 mM EDTA, 1% beta-mercaptoethanol and 0.05% sarkosyl (Sigma) for 20 minutes at room temperature. Arrays were transferred to hybridization buffer containing 2 M NaCl (Fisher Scientific) with 8% (v/v) polyethylene glycol (PEG, Sigma) and agitated for 40 minutes at room temperature, mRNA capture beads with mRNA hybridized were collected from each Seq-Well array, and beads were resuspended in a master mix for reverse transcription containing Maxima H Minus Reverse Transcriptase and buffer, dNTPs, RNase inhibitor, a 5’ template switch oligonucleotide (Seq-Well 5’ TSO), and PEG for 30 minutes at room temperature, and overnight at 52°C with end-over-end rotation. Exonuclease digestion was carried out as described previously: beads were washed with TE with 0.01% tween-20 (Fisher Scientific) and TE with 0.5% SDS (Sigma), denatured while rotating for 5 minutes in 0.2 mM NaOH, and resuspended in ExoI (NEB) for 1 hour at 37°C with end-over-end rotation ( Hughes et al., 2019 ). Next, beads were washed with TE + 0.01% tween-20, and second strand synthesis was carried out by resuspending beads in a master mix containing Klenow Fragment (NEB), dNTPs, PEG, and the dN-SMRT oligonucleotide (Seq-Well Second Strand Primer) to enable random priming off of the beads. PCR was carried out as described using 2X KAPA HiFi Hotstart Readymix and ISPCR primer (Seq-Well ISPCR), and placed on a thermal cycler using the following protocol: 95°C for 3 minutes, followed by 4 cycles of 98°C for 20 seconds, 65°C for 45 seconds, 72°C for 3 minutes, followed by 12 cycles of 98°C for 20 seconds, 67°C for 20 seconds, 72°C for 3 minutes, followed by a final 5-minute extension at 72°C. Post-whole transcriptome amplification proceeded as described above for SMART-seq2 libraries, with the following exceptions: AMPure XP SPRI bead cleanup occurred first at a 0.6 × volume ratio, followed by 0.8x. Library size was analyzed using an Agilent Tapestation hsD5000 kit, confirming the expected peak at ~1000 bp, and absence of smaller peaks corresponding to primer. Libraries were quantified using Qubit High-Sensitivity DNA kit and prepared for Illumina sequencing using Nextera XT DNA Sample Preparation kit using 900 pg of cDNA library as input to tagmentation reactions. Amplified final libraries were purified twice with AMPure XP SPRI beads as before, with a volume ratio of 0.6x followed by 0.8x. Libraries from 3 Seq-Well arrays were pooled and sequenced together using a NextSeq 500/550 High Output v2 kit (75 cycles) using a paired end read structure with custom read 1 primer (Seq-Well CR1P): read 1: 20 bases, read 2: 50 bases, read 1 index: 8 bases.

QUANTIFICATION AND STATISTICAL ANALYSIS Image analysis

All image analyses were performed in Imaris 9.2.1 or 7.4.2 as detailed below. To better visualize neuronal architecture in or/and around LNs, for all LN images except for Figures 7K , 7L , S1B , S2D , S2E , S3A , S3D and S7A , an isosurface for the LN was generated by manually drawing LN contours on 2D slices every fifth slice and was used to mask the original images so that only what was inside the LN mask was shown. Depending on the purpose of the experiment, LN isosurfaces were defined with varying degrees of stringency: based on the outermost layer of LECs in Figures 2A , 2C and S2A , on collagen type I staining in Figure 2F , on SMA staining in Figures 2D , 2E and S2B , or on GFP background staining in Figure S1A or on tdTomato background staining everywhere else. In Figures 2A , 2C , 2D , 2E , 2F , S2A S2B and S3H , additional masking of the channel(s) where nerves were stained was performed with isosurfaces generated for neuronal signal within LNs based on morphology, i.e. fiber-like structures that can be traced through multiple slices, to highlight neuronal structures. To better visualize fibers in the capsular/subcapsular space of LNs as shown in Figure 2F , intranodal sensory fibers and total sensory fibers within and below the capsule were isolated by masking the original channel with LN isosurfaces defined based on GFP (LECs) and collagen type I staining, respectively. The resulting channel after subtracting the former channel from the latter one corresponds to the capsular/subcapsular plexus. Original, processed images and rendered isosurfaces were viewed as 3D reconstructions in surpass view with orthogonal camera setting unless indicated otherwise. For quantification of innervation density of LNs as in Figures 1C – 1F , relevant channels were first masked with the LN isosurface as described above. Isosurfaces for sensory and sympathetic fibers within the masked channels, i.e., inside the LN, were then generated by automatic creation based on features that distinguish neuronal signal from everything else, e.g., intensity, sphericity, followed by manual editing. Sensory or sympathetic fiber density for a given LN was defined as the ratio of the volume of isosurfaces for sensory or sympathetic fibers within the LN to that for the LN. Similar steps were taken to generate isosurfaces for sensory fibers inside the LNs for 3D rendering and quantification in Figures 5G – 5N . For quantification of penetration depth of intranodal sensory fibers, the outermost layer of LECs, which demarcates the LN boundary, was used to precisely segment LNs into isosurfaces. Isosurfaces for intranodal sensory fibers, sensory fibers within the relevant channel after applying the LN isosurface as a mask, were generated as described above. Using the distance transformation function, the closest distance from any given voxel within the LN isosurface to the surface of the LN in μm was computed and converted into an intensity value for that given voxel in a separate channel. To determine penetration depth of intranodal sensory fibers, the distance transformation channel was masked against isosurfaces for intranodal sensory fibers to generate a new channel where the penetration depth at any given voxel within the intranodal sensory fibers was encoded as the intensity value for that specific voxel with the maximum intensity value representing the maximum penetration depth for a given LN. Such a channel, when displayed in surpass view as in Figure 2A , allowed direct visualization of the spatial relationship between intranodal sensory fibers and the nearest LN surface. Additionally, the penetration depth of intranodal sensory fibers was described in Figure 2B in the form of the percentage of total intranodal fibers found within LN spaces with increasing distance away from the LN surface. For that analysis, the original distance transformation channel, as described above, was used to create a series of isosurfaces of decreasing sizes which represent increasingly-deep LN spaces with its closest distance to the LN surface increasing from 0 to 100 μm with 10 μm intervals. For example, 10 was set as the intensity threshold cutoff during automatic creation so that all voxels with intensity value larger than and equal to 10 were selected in one single surface which corresponds to the LN space 10 μm and more below the LN surface. To calculate the percentage of total intranodal sensory fibers in any of those LN spaces, the isosurface for total sensory fibers and that for a said LN space, e.g., 10 μm and more below the surface, as described above, were each used to generate their corresponding binary channels, where all voxels outside of a surface were set as 0, while those inside were set at 100. Colocalization analysis was then performed on those two binary channels, and the percentage of non 0 voxels in the binary nerve channel that were colocalized corresponded to the percentage of total intranodal sensory fibers found 10 μm and more below the LN surface. This process was repeated for increasingly smaller LN spaces with 10 μm interval. The percentage of total intranodal fibers in LN spaces in the form of 10um bins from 0 to 100 μm, e.g., 10–20 μm, was further derived from serial subtraction of the percentage in the LN space that is 10um deeper, e.g., (20 μm and more below the surface), from that in the current LN space, e.g., (10 μm and more below the surface), as shown in Figure 2B . Neuron scRNA-seq data preprocessing Single cells were sequenced to a depth of 1.6 +/− 0.1 million (mean +/− SEM) reads per cell. Pooled libraries were demultiplexed using bcl2fastq (v2.17.1.14) with default settings, and aligned using STAR ( Dobin et al., 2013 ) to the mouse UCSC genome reference (version mm10), and a gene expression matrix was generated using RSEM (v1.2.3) in paired-end mode. Single-cell libraries with fewer than 3,000 unique genes and fewer than 17% of reads mapping to transcriptomic regions were excluded from subsequent analysis, resulting in a final dataset of 52 LN-innervating neurons collected from 8 mice, and 31 skin-innervating neurons collected from 4 mice. Among cells retained for analysis, the number of unique genes captured was 9,843 +/− 229 (mean +/− SEM) among LN-innervating neurons and 9,653 +/− 302 among skin-innervating neurons. Libraries from LN-innervating neurons contained 50.45 +/− 2.3% transcriptome-aligning fragments, libraries from skin-innervating neurons contained 58.33 +/− 2.9 %. Among all alignment and library quality metrics assessed, we found no significant differences between LN-innervating and skin-innervating neurons. All analysis of gene expression was completed using the normalized RSEM output as transcripts per million (TPM).

Neuron scRNA-seq differential gene expression

All analysis of scRNA-seq data was carried out using the R language for Statistical Computing. Single-cell libraries were first assessed for expression of canonical neuronal markers and known lineage-defining genes from accompanying imaging data, such as Nav1.8 ( Scn10a ) and tyrosine hydroxylase ( Th ). The full list of markers is supplied in Table S1 . To directly assess differences in gene expression between LN-innervating and skin-innervating neurons, we used the R package Single Cell Differential Expression (SCDE, version 1.99.1) with default input parameters ( Fan et al., 2016 ; Kharchenko et al., 2014 ). A cutoff of Holm corrected Z score >1.96 or < −1.96 (corresponding to a corrected p-value < 0.05) was used to identify significantly DE genes for subsequent analysis. Heatmaps were created using the R package gplots (version 3.0.1). DAVID was used for analysis of over-represented gene ontologies over significantly DE genes ( Huang da et al., 2009a , b ). Analysis of neuron scRNA-seq in the context of previously published datasets As our target-specific single cells do not represent the full diversity of neurons contained in the DRG, we utilized the scRNA-seq atlas published by Usoskin, Furlan et al. Nature Neuroscience 2015 (subsequently referred to as the “Usoskin, Furlan Atlas” to classify our cells ( Usoskin et al., 2015 ). Using the raw data and accompanying metadata hosted at http://linnarssonlab.org/drg/ , we first identified the intersection of expressed genes from the Usoskin, Furlan Atlas and LN-innervating and skin-innervating single cells, and eliminated cells identified as non-neuronal (“NoN” and “NoN outlier”) from the Usoskin, Furlan Atlas, resulting in a dataset of 148 neurofilamentous (NF), 81 peptidergic (PEP), 251 tyrosine hydroxylase (TH), 169 non-peptidergic (NP), and 39 “Central, unsolved” cells. To mimic the dimensionality reduction methods the previous authors used to identify major neuronal cell types, we transformed the data as log 2 (1+TPM), and calculated the gene variance across all cells. We cut to genes with a variance log 2 (1+TPM) > 0.5, resulting in 11,778 genes. Next, we performed principal component analysis over the log 2 -transformed, mean-centered data, and found that PC 2 and PC 4 reflected major axes of variability between TH, PEP, NF, and NP cell types – identified by the authors of the previous study as “Level.1” cell type subsets ( Figure S4A ). To identify how LN-innervating and skin-innervating cells related to major DRG cell types in a reduced dimensional space, we projected our target-specific data into PC 2 and PC 4 of the Usoskin, Furlan Atlas. This was completed by first calculating the principal components: (Equation 1) X − c m = U S V T where X is the log 2 (1+TPM) data matrix of M genes by N cells from the Usoskin, Furlan Atlas. Equation 1 calculates the singular value decomposition of this matrix after subtracting the average of each row (gene) of X , denoted c m , from X . U represents a matrix of M orthonormal vectors corresponding to M genes and V represents a matrix of N orthonormal vectors corresponding to N cells. To apply this same dimensionality reduction transformation to our new dataset of LN-innervating and skin-innervating single cells, Y , we use Equation 2 : (Equation 2) P C i = ∑ m = 1 M Y − c m u i Y represents the log 2 (1+TPM) transformed matrix of our innervation-target-specific data, c m refers to the same vector of row (gene) averages calculated from X . The centered Y matrix is multiplied as a dot product with the i th principal component gene eigenvector, or the i th column vector of U , denoted u i . By taking the sum over all transformed rows for each column (cell), we project the LN-innervating and skin-innervating data ( Y ) into the principal component space calculated for the Usoskin, Furlan Atlas ( X ), denoted PC i . This data is visualized by plotting the PC 2 and PC 4 vectors from the Usoskin, Furlan Atlas (transparent circles, Figure S4A ), with the PC 2 and PC 4 vectors from the transformed LN-innervating and skin-innervating cells (filled squares). The Euclidean distance between each innervation-target-specific single cell and all cells within the Usoskin, Furlan Atlas was calculated over PC2 and PC4 ( Figure S4B ). The range of cell-to-cell Euclidean distances between like-cells (e.g. PEP-to-PEP) within the Usoskin, Furlan Atlas is represented by a dashed line corresponding to the 99%ile. We similarly analyzed our LN-innervating and skin-innervating single cell transcriptomes in the context of an expanded dataset of somatosensory neurons published by Sharma et al. Nature 2020 (subsequently referred to as the “Sharma Atlas”) ( Sharma et al. 2020 ). Here, authors generated data from 10,922 cells from cervical, thoracic, and lumbar dorsal root ganglia, uncovering at least 15 transcriptionally distinct cell types/states. Unlike the Usoskin, Furlan Atlas, the cell types in the Sharma Atlas could not be sufficiently represented in two-dimensional principal components space. When we used the JackStraw method from the Seurat package to identify significant principal components, resulting in 32 significant principal components. As with the Usoskin, Furlan Atlas, we next projected the LN-innervating and skin-innervating single cell transcriptomes into the 32-dimensional principal components space defined by the Sharma Atlas according to Equation 2 , and created a joint principal component matrix of cells × PC 1–32 including cells from the Sharma Atlas, LN-innervating neurons, and skin-innervating neurons. Finally, we constructed a shared nearest-neighbor graph over this joint matrix and visualized single cells in a two-dimensional uniform manifold approximation and projection (UMAP) plot ( Figure 4A , 4B ). To establish cell similarity, we used label propagation methods over a shared nearest-neighbor graph built from the 32-dimensional projected space, similar to other published methods for single-cell transcriptomic integration ( Barkas et al. 2019 , Stuart et al. 2019 ). Here, “unknown” cells (from the LN-innervating or skin-innervating datasets) were each queried for their top 100 nearest neighbors, and cellular identity was predicted based on a consensus identity of the cell’s neighbors. Results from label propagation are presented in Figure 4C . To analyze the expression similarity between each single cell from our target-specific dataset and the data in the Usoskin, Furlan Atlas and Sharma Atlas in a more directed, supervised manner, we assessed how each single cell correlated with each neuron subtype defined by each Atlas. For the Usoskin, Furlan Atlas, we elected to use the more detailed neuronal subtypes, termed “Level.3”, which breaks some of the major neuron subtypes, NP, PEP, and NF, into subtypes based on intra-population diversity. For both Atlases, we calculated the average gene expression for each neuron subtype (e.g. NP1) over the log 2 (1+TPM) transformed single-cell data, generating pseudo-population averages for each neuron subtype. Next, we only considered genes in our pseudo-population averages that were designated as “subtype-defining” by the Usoskin, Furlan et al. analysis, corresponding to the top 50 genes upregulated within each cell type when compared to all other cell types in their dataset, yielding 379 unique genes. For the Sharma Atlas, we restricted genes to the 700 genes defined by the authors as “subtype specific” and provided as supplementary data 1 in their manuscript. We similarly restricted our LN-innervating and skin-innervating single-cell libraries to only these 379 unique genes (for comparison to the Usoskin, Furlan Atlas) or 700 unique genes (for comparison to the Sharma Atlas), and calculated the Spearman correlation between each target-specific single cell (following log 2 (1+TPM) transformation) and the pseudo-population averages. The results from this analysis are presented in Figure 4D for comparison to the Sharma Atlas, and Figure S4C for comparison to the Usoskin, Furlan Atlas. We clustered LN-innervating and skin-innervating single cells by their correlation with the Usoskin, Furlan Atlas pseudo-population using complete linkage clustering, and using a cut height of 0.8 retained 4 distinct Neuron Types: Neuron type 1 (LN-innervating cells: 25, skin-innervating cells: 9), Neuron type 2 (LN-innervating cells: 1, skin-innervating cells: 14), Neuron type 3 (LN-innervating cells: 23, skin-innervating cells: 5), and Neuron type 4 (LN-innervating cells: 3, skin-innervating cells: 3) ( Figure S4C ). We confirmed these clusters were largely stable and in agreement with the partitions created in comparison with the Sharma Atlas ( Figure 4D ). Neuron Types 1 and 3, which share the strongest similarity with CGRP-η, -ξ, -α, -β -γ, and -ε subtypes defined by Sharma et al. (and highest similarity with PEP1 and PEP2 subtypes defined by Usoskin, Furlan et al.), were enriched in the LN-innervating population relative to the skin-innervating population (LN-innervating: 48% Neuron Type 1, 44% Neuron Type 3; skin-innervating: 29% Neuron Type 1, 16% Neuron Type 3). Conversely, Neuron Types 2 and 4 were underrepresented in the LN-innervating population compared to the skin-innervating population (LN-innervating: 2% Neuron Type 2, 6% Neuron Type 4; skin-innervating: 45% Neuron Type 2, 10% Neuron Type 4). To assess the gene expression phenotype of each Neuron Type, we used SCDE to identify DE genes between cells of each Neuron Type compared to all cells of the 3 remaining Neuron Types. SCDE was run as described above with default input parameters, genes with a Holm-corrected p-value < 0.01 were considered significant and presented in Figure 4F and Table S1 . LN Seq-Well data preprocessing Reads were aligned and processed according to the Drop-Seq Computational Protocol v2.0 ( https://github.com/broadinstitute/Drop-seq ). Briefly, reads were first demultiplexed according to index read 1 using bcl2fastq (v2.17.1.14) with default settings. Read 1 was split into the first 12 base pairs corresponding to the cell barcode (CB), and the 13–20 th base pairs, which encode the unique molecular identifier (UMI). CBs, UMIs, and read 2 sequences with low base quality were discarded, as were any that contained non-random sequences (e.g. primer sequences, poly-A tails). Following CB and UMI tagging, read 2 was aligned to the mouse genome (version mm10) using STAR v2.5.2b with default parameters including “--limitOutSJcollapsed 1000000 --twopassMode Basic”. STAR alignments were merged to recover cell and molecular barcodes, and any sequences within hamming edit distance 1 were merged, as these likely originated from the same original sequence. Additional methods to correct for bead synthesis errors in the CB or UMI are detailed in the Drop-Seq Computational Protocol v2.0 (“DetectBeadSynthesisErrors” function). Digital gene expression matrices for each array were retained following quality filtering and UMI-correction, and further processed using the R language for Statistical Computing. Cells with fewer than 300 unique genes were removed from analysis. Dimensionality reduction, clustering, visualization, and cell type identification of LN Seq-Well data We restricted our primary analysis of LN-resident cell types to only arrays corresponding to steady state iLN without surgical manipulation or optogenetic stimulation. A total of 9,662 cells were retained with 25,929 unique genes expressed across 7 mice with 1 LN per mouse. For 2 mice, we sequenced arrays corresponding to all LN cells prior to CD3ε/CD19 depletion as well as CD3 ε /CD19 depleted cells on a separate array. The average cell recovery per array was 1,074 +/− 141 (mean +/− SEM) cells, with an average gene count of 1,581 +/− 11 genes and average UMI per cell of 4,251 +/− 48 UMI (mean +/− SEM). Data was normalized and scaled using the Seurat R package ( https://github.com/satija.lab/seurat ) ( Butler et al., 2018 ): transforming the data to log e (UMI+1) and applying a scale factor of 10,000. We confirmed equivalent depth and cell quality across each of our arrays and the absence of major batch effects introduced by sequencing work-up day or other technical factors, and thus did not regress any batch-related covariates out of our data, including individual cell quality or mitochondrial percent. To identify major axes of variation within our data, we first subsetted our data to only highly-variable genes across all cells – all genes with dispersion (calculated as the variance to mean ratio) > 1.1 were kept, resulting in 2,348 variable genes. Principal component analysis was applied to the cells cut to variable genes for the top 100 principal components. Using the JackStraw function within Seurat, we identified the top significant PCs, and compared these significant PCs to the variance explained by each dimension, ultimately choosing 41 PCs for subsequent clustering and further dimensionality reduction ( Shekhar et al., 2016 ). Critically, we completed all of the following analysis over a range of variable gene cutoffs and principal components to ensure that our cell identification results were robust to parameter choice. For 2D visualization, we used the Barnes-Hut implementation of t-distributed stochastic neighbor embedding (t-SNE) with “perplexity” set to 40. This tSNE projection of the steady state LN atlas is represented in Figure 6B , S5A . To identify clusters of transcriptionally-similar cells, we employed unsupervised clustering with the Louvain algorithm with the Jaccard correction ( Blondel et al., 2008 ). Briefly, this method involves constructing a k-nearest neighbor graph over the Euclidean distance between cells in the 41-PC reduced space, followed by a shared nearest neighbor (SNN)-based clustering and modularity optimization ( Waltman and Van Eck, 2013 ). We implemented this using the FindClusters tool within the Seurat R package with default parameters and k.param set to 20 and resolution set to 0.4. Here, we intentionally underclustered our data to avoid erroneously splitting cells with shared cell type functions, as the variable genes calculated for this dimensionally-reduced space likely did not fully reflect more nuanced cell type differences (e.g. variable behavior between Neutrophil subtypes). The “Parent Cluster” results from first-pass cell type clustering are represented in the tSNE plot and clusters identified in Figure S5A . We used the Seurat function FindAllMarkers to identity differentially-expressed genes upregulated within each cluster compared to all other cells in the dataset, and tested differential expression using the likelihood-ratio test for single-cell gene expression (by setting test.use to “bimod”) ( McDavid et al., 2013 ). The top 100 differentially-expressed genes for each cluster were analyzed, as ranked by the average fold change and restricted to only those with FDR-corrected p-values < 0.05. Next, to assess if any cell subtypes existed within each cluster, we restricted our data to only cells within a single “Parent Cluster”, and recalculated the variable genes over these cells. The above analysis, from calculation of variable genes to tSNE visualization and cluster identification, was repeated for each cluster listed in Figure S5A . Cell types for which we could identify sub-clusters with significant differentially-expressed genes are marked with asterisks next to their names in Figure S5A , and the sub-cluster tSNE projections and top differentially expressed genes are represented in Figures S5B – S5O . For the T cell parent cluster, we required two iterative sub-clustering steps to fully enumerate all constituent cell types: the first clustering step differentiated regulatory T cells (Tregs) from the remaining T cells ( Figure S5B , S5C ), and subsequent clustering on the non-Treg T cells uncovered CD4 T cells vs. CD8 T cells. All differentially-expressed genes within each sub-cluster can be found in Table S2 . After exhaustive assessment for cell subclusters within each cell type, we identified 24 unique cell types within our steady state dataset ( Figure 6B ). We calculated the differentially expressed genes between each cell type and all other cells using a likelihood ratio test (using the FindAllMarkers function with test.use set to “bimod”), the results of this analysis are presented in Figures 6C , S5P , and Table S2 . By identifying canonical marker genes within these DE gene lists from the literature and using resources such as ImmGen ( Heng et al., 2008 ), we attributed cell identities to each cell type within our dataset, as named in Figures 6B , 6C , and S5P . We recovered two populations of blood endothelial cells, one of which likely corresponds to HEV endothelial cells (HEC), based on uniform and high expression of well-established HEC markers, e.g., Ccl21a , Fut7 , Chst4 ( Figures S5H , S5I ) ( Homeister et al., 2001 ; Kawashima et al., 2005 ; Stein et al., 2000 ; Uchimura et al., 2005 ). The second blood endothelial cell cluster, termed BEC, includes a heterogeneous set of non-HEV endothelial cells, including arterial, capillary, and non-HEV venular cells ( Figures S5H , S5I ) ( Vanlandewijck et al., 2018 ). Additionally, we identified two distinct populations of Prox1 + lymphatic endothelial cells (LECs), with LEC 1 likely representing a mixture of LECs lining the floor of SCS and medullary sinuses based on Madcam1, Msr1, Bmp2, Vcam1, and CD274 (PD-L1) expression profiles ( Cohen et al., 2014 ; Cordeiro et al., 2016 ; Takeda et al., 2019 ), and LEC 2 defined by unique expression of multiple extracellular matrix or structural proteins, including Fbln2, Aqp1, Fbln5, Tnc and Reln ( Figures S5N , S5O ). Interestingly, cells within the LEC 2 cluster shared molecular signatures with LECs positioned at the ceiling of SCS, within the cortical sinuses, and lining lymphatic vessels, including selective expression of Emcn , Ackr4 , Ackr2 , Klf2 , Fabp4 and Cav1 ( Iftakhar et al., 2016 ; Takeda et al., 2019 ; Ulvmar et al., 2014 ). We also identified a subtype of dendritic-cell-like cells (Aire + APC), which likely represent the Aire-expressing ILC3-like cells that were described recently ( Yamano et al., 2019 ). Similarly, we detected two subtypes of neutrophils, Neutrophils 1 and 2, which potentially reflect distinct maturation states similar to what was recently described for bone marrow neutrophils ( Evrard et al., 2018 ) ( Figure S5L and S5M ): Neutrophils 1 expressed components of neutrophil granules and effector molecules at high levels, including Elane, Prtn3, Ctsg, Ngp, Ltf, Camp, and Mpo , whereas Neutrophils 2 expressed little or no effector molecules, but elevated levels of pro-inflammatory genes, including leukocyte traffic molecules, chemokines, cytokines, and cytokine receptors such as Sell, Ccl4, Cxcr2, Cxcl2, Ccl6, Il1b, and Csf3r .

Analysis of cellular receptor-ligand pairs

We reasoned that cells or cell types within the LN that interact with innervating neurons would likely express proteins that enable such contact or communication. As we generated unbiased single-cell transcriptomic data from LN-innervating neurons and the potential targeted cell types, we incorporated databases of ligand and receptor pairs to understand if any of the LN-resident cell types expressed a high abundance of cognate molecules, and would thus be poised to interact with innervating neurons. A general schematic of this method is provided in Figure S6A . We used the database of receptor-ligand interactions curated by Ramilowski et al ( Ramilowski et al., 2015 ), which consists of 2,422 total interactions over 708 unique genes (originally provided as human genes, and converted to mouse orthologs using the HUGO database). First, data from LN-innervating neurons was limited to only genes with non-negligible expression, using a cutoff of average log 2 (1+TPM) > 3, yielding 6,666 total genes for subsequent analysis. The intersection of genes within the Ramilowski interaction database and those expressed at non-negligible levels among LN-innervating neurons yielded 184 total genes. After limiting to only interactions with at least one participating gene expressed in the LN-innervating neurons, the interaction database was restricted to 750 total receptor-ligand pairs, and 471 unique potential cognates. We next assessed the expression of these 471 cognate genes within the LN-resident cell atlas. First, we summarized the expression of individual cells within the LN-resident atlas by taking the pseudo-population average of each cell type (over non-log single-cell data). We limited the LN-resident atlas data to only genes with non-negligible expression across all cell type pseudo-populations, cutting to genes with an average UMI expression > 1, yielding 256 total potential cognates (from the previous 471). Next, we developed a summary statistic to reflect the abundance of neuron cognates expressed within LN-resident cell types. First, we scaled our data by subtracting the mean and dividing by the standard deviation for each individual gene – this enabled us to assess the contribution of all genes equally such that signal was not dominated by genes with high total expression ( Figure 6E ). Finally, we calculated the “Interaction Potential” (IP) as the mean of these scaled values for each cell type: cell types that expressed relatively higher abundances of all candidate neuron-cognates received a higher IP score. Our null model states that the interaction potentials we calculated are no more extreme than the IP we would have recovered by chance. To test our experimentally-derived IP, we generated a null distribution by shuffling the cell type labels over all single cells within the LN-resident cell atlas, and repeated the “cell type” averaging, scaling, and IP calculation for 1,000 permutations. By comparing our true IP scores to the null distribution, we were able to identify certain cell types with significantly higher IP than observed by chance, and could attribute a P-value to each cell type ( Figure 6F , 99% confidence interval denoted by dashed line). IP scores were re-scaled such that the lower bound of the 99% confidence interval was equal to 0 for clarity. The results of this approach are presented in Figures 6D – 6F , S6A , S6B . Crucially, we were concerned that the method of calculation of the IP, the summary statistics applied, the choice of raw vs. scaled data, or confounding factors that differentiate cell types, including average genes/cell and number of cells per cell type, would influence our ranked list of top interacting cell types and bias our results. For example, we wondered whether differences in quality metrics or other technical factors between cell types might result in higher or lower IP rankings – for instance, a cell type with significantly higher RNA recovery per cell than another cell type would appear to have a higher interaction potential. We found no correlation between the IP (as reported in Figure 6F ) and the median UMI per cell for each cell type ( Figure S6B , p = 0.32). To address bias introduced by our choice in summary statistic or data normalization, we repeated the above pipeline without gene-wise scaling across cell types ( Figure S6C ), or by calculating the percent of cells with non-zero expression of a given gene, in the place of calculating of average expression per cell type ( Figure S6D ). In both of these cases, we observed that non-endothelial stroma, LEC 1, LEC 2, BEC, and HEC remained the top-scoring cell types for Interaction Potential (significance calculated by permutation test as described above). Finally, we reasoned that variations in the number of cells per cell type might limit our ability to compare between different cell types. We iteratively down-sampled our single-cell data to analyze interaction potentials (using the method in Figures 6D – 6F ) for only 25 total cells per cell type – the histograms of these calculations after 1,000 iterations are plotted in Figure S6E . Critically, non-endothelial stroma, LEC 1, LEC 2, BEC, and HEC cell types remained top-ranking in Interaction Potential after controlling for cell abundance per cell type. Finally, we derived an alternative statistical testing strategy to assess the overrepresentation of neuron-interaction cognates among expressed genes between different cell types. Here, we binarized our data to classify genes as “expressed” or “not expressed” within a cell type, using an average gene expression cutoff of 1. We considered the list of 256 potential neuronal cognate genes, and used a Fisher’s Exact Test to assess whether the cognate gene list was overrepresented among expressed genes for a given cell type (mimicking the field-standard for gene ontology enrichment analysis) ( Huang da et al., 2009a , b ), and a Holm correction to adjust for multiple tests. In close agreement with the results from our interaction potential statistic above, we found significant overrepresentation of potential neuronal cognate genes in the following cell types (listed in decreasing statistical significance): non-endothelial stroma (p = 1.6 × 10 −28 ), BEC (p = 2.5 × 10 −22 ), LEC 1 (p = 4.5 × 10 −22 ), HEC (p = 8.3 × 10 −21 ), LEC 2 (p = 9.6 × 10 −20 ), Macrophages (p = 8.7 × 10 −9 ), Mast Cells (p = 6.5 × 10 −8 ), Neutrophils 2 (p = 5.2 × 10 −6 ), Neutrophils 1 (p = 1.8 × 10 −4 ), pDC (p = 1.7 × 10 −3 ), Aire + APCs (p = 3.4 × 10 −3 ), and cDC2 (8.9 × 10 −3 ). All other cell types were non-significant by a Holm-adjusted p-value cutoff of 0.01. Critically, this ranking was not sensitive to the choice of binarization cutoff, tested over a range of 0.5 – 10 UMI, data not shown).

Differential gene expression following optogenetic stimulation

Cells were partitioned into the cell types annotated in Figure 7B . Using the Seurat function DiffExpTest, which employs a likelihood ratio test to identify differentially expressed genes, we analyzed cells for each cell type from ChR2+Light+ LN vs. ChR2+Light- LN. Similarly, we identified differentially expressed genes by cell type between ChR2-Light+ LN vs. Chr2-Light- LN. We reasoned that the DE genes in ChR2+ mice represented both the effects of neuronal stimulation, as well as changes induced by surgery and/or phototoxicity, while the DE genes in the ChR2− mice only correspond to changes due to surgery and/or phototoxicity. For each cell type, we identified genes DE in ChR2+ animals by a Holm-adjusted p-value cutoff of 0.05, and eliminated genes from these lists that were also DE (using the same cutoff) in ChR2− LN. We calculated the effect size using Cohen’s d, and restricted our gene lists to only those genes with a non-negligible effect size, using a cutoff of 0.2 (analysis the effect of various effect-size cutoffs in Figure S7F ). The results of these analyses for each cell type can be found in Table S3 . In Figure 7H , we further restricted our DE gene lists for heatmap visualization, and in Figure 7J for gene ontology analysis (using DAVID, as described above) by only including genes that were also DE between LNs harvested from the same mouse in at least 2 of 4 ChR2+ mice.

Statistical testing

Using prism software, we performed unpaired two-tailed Student’s t-tests for Figures 1E , 1H and Figure 5F , Welch’s t-tests for Figures 5G , 5H , 5K , 5L , S1G , S3C , S4N , S4O , 2-way ANOVA with Sidak’s multiple comparisons test with for Figure S7E . All other statistical tests corresponding to differential gene expression or assessment of interaction potential are described above and completed using R language for Statistical Computing. Tests of correlation and correlation significance are annotated by the correlation model used (Pearson vs. Spearman) were completed using R language for Statistical Computing. Parameters such as sample size, number of replicates, number of independent experiments, measures of center, dispersion, and precision (mean +/− SEM) and statistical significances are reported in Figures and Figure Legends. A p-value less than 0.05 was considered significant unless otherwise reported; a more stringent cutoff of 0.01 was used in some instances, and annotated as such. Where appropriate, a Holm correction was used to account for multiple tests, as noted in the figure legends or STAR Methods .

Materials Availability

No new or unique reagents were generated in this study.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Mouse lines used in this study were all previously described and include Nav1.8 Cre (RRID:IMSR_EM:04582) ( Nassar et al., 2004 ), Rosa26 LSL-tdTomato (RRID:IMSR_JAX:007914), Bmx-CreER T2 (MGI:5513853) ( Ehling et al., 2013 ), Rosa26 LSL-DTA (RRID:IMSR_JAX:009669), Prox1-EGFP (MGI:4847348) ( Choi et al., 2011 ), Rosa26 LSL-ChR2-eYFP (RRID:IMSR_JAX:024109), Rosa26 LSL-eYFP (RRID:IMSR_JAX:007903) and ChAT BAC -eGFP (RRID:IMSR_JAX: 007902). All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols of Harvard Medical School. Unless indicated otherwise, adult mice of both sexes between 6–12 weeks of age were used for various experiments.

METHODS DETAILS Whole mount immunohistochemistry

Whole mount immunohistochemistry of LNs was performed using an iDISCO protocol with methanol pretreatment optimized for LNs ( Renier et al., 2014 ). Briefly, adult animals (6–12 weeks) were perfused with 25 mL of PBS (Hyclone) and 25 mL of 4% paraformaldehyde (PFA, Sigma) sequentially at room temperature (RT). Peripheral lymph nodes (PLNs), including popliteal and inguinal lymph nodes (popLNs and iLNs), were postfixed with 4% PFA for 1 hr at 4°C. For methanol pretreatment, fixed LNs were washed sequentially in 50% methanol (Fisher Scientific) (in PBS) for 1 hr, 100% methanol for 1 hr, 50% methanol for 1 hr, PBS for 1 hr twice, and PBS/0.2% Triton X-100 (VWR) for 1 hr twice at RT. LNs were then left in PBS/0.2% Triton X-100/20% DMSO (Sigma)/0.3 M glycine (BioRad) overnight at RT and blocked in PBS/0.2% Triton X-100/10% DMSO/6% donkey serum (Jackson Immunoresearch) or goat serum (Gibco)/anti-CD16/CD32 (Fc block) (Bio X cell) overnight at RT. LNs were subsequently washed in PBS/0.2% Tween-20 (Fisher Scientific)/10 μg/mL heparin (Sigma) (PTwH), for 1 hr twice at RT, before incubation with antibody mix in PTwH/5% DMSO/3% donkey or goat serum/Fc block 1:100 for 3 days at RT. LNs were extensively washed in PTwH for at least 6 times over the course of a day at RT. For unconjugated antibodies, LNs were further incubated with a secondary antibody mix including a panel of species-specific anti-IgG (H+L) Alexa Fluro 488, 546, 647 and 594-conjugated antibodies (Invitrogen or Jackson Immunoresearch) in PTwH/5% DMSO/3% donkey or goat serum/Fc block 1:100 for 3 more days at RT. LNs were washed in the same way as after primary antibody incubation for 1 day. Immunolabeled LNs following one round of antibody incubation for conjugated antibodies (or two for unconjugated antibodies) were then processed for clearing, which includes sequential incubation with 50% methanol for 1 hr, 100% methanol for 1 hr for three times and a mixture of 1-part benzyl alcohol (Sigma): 2-parts benzyl benzoate (Sigma) (BABB) overnight at RT. For tdTomato immunolabeling, goat anti-mCherry antibody (ACRIS) was preabsorbed against PLNs from tdTomato − animals overnight at RT prior to use. Whole mount immunohistochemistry of DRGs and the skin was performed as described previously ( Li et al., 2011 ). Briefly, DRGs inside vertebral column and the depilated hairy skin from PFA-perfused animals (6–12 weeks) were postfixed with 4% PFA for 1 hr or Zamboni fixative (Fisher Scientific) overnight, respectively at 4°C. Samples were washed every 30 min with PBS/0.3% Triton-100 (0.3% PBST) for 4–6 hr, then incubated with primary antibodies in antibody diluent (0.3% PBST/20% DMSO/5% donkey or goat serum) for 2–3 days at RT. Samples were then washed with 0.3% PBST every 30 min for 5–8 hr before incubation with secondary antibodies in antibody diluent for 2–3 days at RT. After extensive washes as described above, samples were dehydrated and cleared in 50% methanol for 1 hr, 100% methanol for 1 hr for three times and BABB overnight at RT. Cleared whole-mount tissues were imaged in BABB between two coverglasses using Olympus FV3000 confocal imaging system, except for those shown in Figures 7K and 7L , which were acquired on BioRad 2100MP system and those shown in Figures 3B , S3B and S3E , which were acquired on Zeiss Stereo Discovery V16. The antibodies used were: rabbit anti-CGRP (Immunostar, 24112, 1:500), chicken anti-GFP (Aves Labs, GFP-1020, 1:500), chicken anti-NF200 (Aves Labs, NFH, 1:500), rabbit anti-Tyrosine Hydroxylase (Millipore, AB152, 1:500), goat anti mCherry antibody (1:500, ACRIS AB0040-200), rabbit anti-βIII-Tubulin (Biolegend, 802001, 1:500), Alexa Fluor 647-conjugated rat anti-CD31 (Biolegend, 102416, 1:50), FITC–conjugated mouse anti-smooth muscle actin (aSMA) (Sigma, F3777-.2ML, 1:500), eFluor 660-conjugated mouse anti-smooth muscle actin (aSMA) (Thermo Fisher, 50-9760-82, 1:100), eFluor 660-conjugated rat anti-CD169 (Thermo Fisher, 50-5755-80, 1:50), Pacific Blue-conjugated rat anti-CD45 (Biolegend, 103126, 1:50), Alexa Fluor 488-conjugated rat anti-PNAd (Thermo Fisher, 53-6036-82, 1:50) Retrograde labeling of LN- and skin-innervating neurons To retrogradely label LN-innervating neurons, adult animals (6–12 weeks) were anesthetized by intraperitoneal injection of ketamine (Patterson Vet) (50 mg kg −1 ) and xylazine (Patterson Vet) (10 mg kg −1 ). The skin overlying the targeted iLN was shaved and depilated so that the LN underneath was visible percutaneously. A 5 mm incision was made directly on top of the iLN. The iLN was microdissected without perturbing afferent lymphatic vessels and surrounding blood vessels. 1 μl of Adeno-Associated Virus (AAV) (AAV2/1.CMV.HI.eGFP-Cre.WPRE.SV40, titer >=8E+12 vg/mL, Addgene) mixed with 0.5 μl of fast green (Sigma) was injected into the iLN of Rosa26 LSL-tdTomato/LSL-tdTomato animals using a pulled and trimmed glass pipette (FHC) which was connected to a 5 mL syringe through the aspiration assembly system (Sigma). The injection site was immediately rinsed with 2 mL of saline (Patterson Vet) to wash away any off-target virus before the incision was closed with sutures. Animals were sacrificed between 1 month and 6 months after injection for histology or scRNA-seq. To directly visualize the axonal projections of sensory neurons retrogradely labeled from the iLN, AAV carrying Cre-dependent tdTomato cassette (AAV2/1.CAG.Flex.tdTomato.WPRE.bGH, titer ≥10 13 vg/mL, Addgene) was injected into the iLN Nav1.8 Cre/+ animals as described above. For WGA-based retrograde labeling, 1 μl of WGA-AF488 (2 mg/mL in PBS, Invitrogen) was injected into the iLN of Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ animals as described before and the animals were processed for histology 4 days post injection. Retrograde labeling of skin-innervating neurons was described previously ( Kuehn et al., 2019 ). Briefly, following ketamine-xylazine mediated anesthesia, a single injection of 0.2 μl of various AAV2/1 viruses as described above and 0.1 μl of fast green was delivered using the injection device described above intradermally into the patch of depilated skin overlying the iLN of adult mice (6–12 weeks). Animals were sacrificed between 1 month and 6 months after injection for immunohistochemistry, RNAscope, or scRNA-seq.

Immunohistochemistry of tissue sections

Adult animals (6–12 weeks) were perfused with 25 mL of PBS and 25 mL of 4% PFA sequentially at RT. The intact vertebral column was postfixed overnight with 4% PFA at 4°C. DRGs were subsequently dissected and processed for cryosectioning. 14 μm serial cryosections were collected and processed for immunohistochemistry as described previously ( Li et al., 2011 ). In brief, sections were postfixed with 4% PFA for 10 min at RT. Following three washes with PBS, they were incubated with blocking buffer (PBS with 5% normal goat serum and 0.3% Triton-100) for 1 hr at RT. The sections were then incubated with Rabbit anti-TH (Millipore) in the same blocking buffer overnight at 4°C. The following day, sections were washed three times with wash buffer (PBS with 0.3% Triton-100) before incubation with goat Alexa Fluor 647-conjugated anti-rabbit (Invitrogen) for 1 hr at RT. Sections were then washed for three times with wash buffer before mounting in Fluoromount Aqueous Mounting Medium (Sigma). WGA-488 and tdTomato were visualized directly based on endogenous fluorescence. All the sections with tdTomato + cells were imaged at 20x using Olympus FV3000 confocal imaging system. Intravital two-photon microscopy Adult Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ animals (6–12 weeks) were given 1 μg of FITC-conjugated rat anti-CD169 antibody (BioRad) diluted in a total volume of 20 μl of PBS into the right footpad to label CD169 + subscapular macrophages inside the draining LN. Immediately after, the animals were prepared microsurgically for intravital two-photon microscopy as described before ( Mempel et al., 2004 ). Briefly, anesthesia during surgical preparation and imaging was achieved through the ketamine-xylazine method as described above. The right popLN was exposed and positioned with the cortex facing outwards with minimal perturbation to afferent lymphatic vessels and surrounding blood vessels, while the animal was immobilized onto a custom-built stage by its hip bone and the vertebral column. The imaging chamber was created around the exposed LN with high vacuum grease (VWR) on the side and a coverslip on top. A thermocouple (Omega) was placed next to the LN to monitor the local temperature, which was maintained between 36.5 and 37°C by a custom-built water bath heating system. Two-photon imaging was performed on a Bio-Rad Radiance 2100MP Confocal/Multiphoton microscopy system with two MaiTai Ti:sapphire lasers (Spectra-Physics) tuned to 800 nm and 900 nm for two photon excitation and second harmonic generation. Z-stacks of sensory innervation of the capsular/subcapsular space on the cortical side were acquired in 1 μm steps with a 20×, 0.95 numerical aperture objective (Olympus). Subcutaneous challenge by footpad injection of LPS or Pam3CSK4 LPS (Sigma) and Pam3CSK4 (Invivogen) were dissolved in water to a stock concentration of 5mg/ml and 2mg/ml respectively. 50μg of LPS or 40μg of Pam3CSK4 were injected into one footpad of adult Nav1.8 Cre/+ × Rosa26 LSL-tdTomato/LSL-tdTomato mice in a total volume of 20μl, while the contralateral footpad received equal volume of PBS. 3 days later, popLNs were harvested and processed for whole mount immunohistochemistry as described above. Manual cell sorting for scRNA-seq Adult mice with retrogradely-labeled LN- or skin-innervating neurons were sacrificed by CO 2 asphyxiation. T13 and L1 DRGs ipsilateral to the side of injection were quickly removed without nerves attached and checked for tdTomato labeling in cold HBSS (1X, no Ca 2+ or Mg 2+ ) (VWR) under Leica MZ10 F stereomicroscope with fluorescence. DRGs were immediately digested with 1 mL of papain solution (HBSS/10 mM HEPES (VWR)/500 μM EDTA (Westnet)/0.4 mg/mL L-Cysteine (Sigma)/1.5 mM CaCl 2 (Sigma)/20 unit/mL Papain (Worthington)) in a 37°C water bath for 10 min, with agitation every 2 min. DRGs were further digested with 1 mL of collagenase type II/dispase solution (HBSS/10 mM HEPES/4 mg/mL collagenase type II (Worthington)/5 mg/mL dispase (Thermo Fisher)) in a 37°C water bath for 30 min, with agitation every 10 min. Following centrifugation at 400 g for 4 min, digested DRGs were mechanically disrupted in 0.2 mL of complete L15 medium (L15 (Invitrogen)/10 mM HEPES/10% FBS (Germini)) by passing them first through a 1000 μL pipette tip up to 10 times, and then through a 200 μL pipette tip up to 5 times until the tissues were fully dissociated. To remove myelin/axonal debris, the cell suspension diluted in 1 mL of complete L15 medium was carefully layered on top of 5 mL of Percoll gradient (L15/10 mM HEPES/20% Percoll (GE Healthcare) and centrifuged at 400 g for 9 min. After removing the supernatant, cells were washed in 2 mL of L15/10 mM HEPES and centrifuged at 750 g for 3 min. Finally, cells were resuspended in 1 mL of cold sorting buffer (L15/10 mM HEPES/1 mg/mL BSA (VWR)/25 μg/mL DNase I (Roche)), and subjected to fluorescence-assisted single-cell picking as described previously ( Hempel et al., 2007 ). Briefly, the cell suspension diluted in 3 mL of sorting buffer was immediately transferred to a 35 mm petri dish (Scanning dish) with lane markings 6 mm apart and let sit on ice until most cells had settled to the bottom which normally takes 15–20 min. Rare fluorescent cells were readily identified under Leica MZ10 F stereomicroscope with fluorescence (transillumination off) by scanning the bottom of the dish lane by lane to maximize recovery and avoid rescanning. Zoom was set such that the field of view corresponded to the width of a single lane. To pick out fluorescent cells with minimal contamination from nonfluorescent cells, a pulled and trimmed micropipette (World Precision Instruments) was carefully lowered under transillumination into the sorting buffer until it was in the vicinity of the target cell. Simultaneous positive pressure was applied by mouth through the aspiration assembly system, as described above for retrograde labeling. Once the micropipette was in position, the target cell was gently aspirated into the micropipette through capillary action by transient release of positive pressure. The micropipette was quickly removed to prevent aspiration of unwanted cells or debris. The content of the micropipette, including the target cell, was expelled gently into a droplet of cold fresh sorting buffer on a different 35 mm petri dish (wash dish 1) under transillumination. Wash dish 1 was kept on ice while subsequent scans for fluorescent cells occurred. Once 16 or all the fluorescent cells, whichever comes first, were collected in wash dish 1, cells were washed two additional times by moving them one by one into a new droplet of sorting buffer on clean 35 mm petri dishes. Micropipettes were not reused for different cells to avoid cross contamination. After the final wash, each fluorescent cell was pipetted up and down the micropipette three times to remove unwanted contamination before being ejected into 10 μl of cold RLT (Qiagen) supplemented with 1% β-mercaptoethanol (Sigma) in a 96-well plate, and snap-frozen on dry ice and stored at −80°C. The entire manual sorting procedure was routinely completed in 1.5 hr. scRNA-seq of neurons using Smart-Seq2 Single-cell libraries were generated according to the SMART-seq2 protocol ( Picelli et al., 2014 ; Trombetta et al., 2014 ). Briefly, RNA from single-cell lysates was purified using AMPure RNA Clean Spri beads (Beckman Coulter) at a 2.2× volume ratio, and mixed with oligo-dT primer (SMART-seq2 3’ Oligo-dT Primer), dNTPs (NEB), and RNase inhibitor (Fisher Scientific) at 72°C for 3 minutes on a thermal cycler to anneal the 3’ primer to polyadenylated mRNA. Reverse transcription was carried out in a master mix of Maxima RNaseH-minus RT enzyme and buffer (Fisher Scientific), MgCl 2 (Sigma), Betaine (Sigma), RNase inhibitor, and a 5’ template switch oligonucleotide (SMART-seq2 5’ TSO) using the following protocol: 42°C for 90 minutes, followed by 10 cycles of 50°C for 2 minutes, 42°C for 2 minutes, and followed by inactivation at 70°C for 15 minutes. Whole transcriptome amplification was achieved by addition of KAPA HiFi HotStart ReadyMix (Kapa Biosystems) and IS PCR primer (ISPCR) to the reverse transcription product and amplification on a thermal cycler using the following protocol: 98°C for 3 minutes, followed by 21 cycles of 98°C for 15 seconds, 67°C for 20 seconds, 72°C for 6 minutes, followed by a final 5-minute extension at 72°C. Libraries were purified using AMPure XP SPRI beads at a volume ratio of 0.8x followed by 0.9x. Library size was assessed using a High-Sensitivity DNA chip (Agilent Bioanalyzer), confirming the expected size distribution of ~1000–2000 bp. Tagmentation reactions were carried out with the Nextera XT DNA Sample Preparation Kit (Illumina) using 250 pg of cDNA per single cell as input, with modified manufacturer’s instructions as described. Libraries were purified twice with AMPure XP SPRI beads at a volume ratio of 0.9x, size distribution assessed using a High Sensitivity DNA chip (Agilent Bioanalyzer) and Qubit High-Sensitivity DNA kit (Invitrogen). Libraries were pooled and sequenced using NextSeq500/550 High Output v2 kits (75 cycles, Illumina) using 30–30 paired end sequencing with 8-mer dual indexing. RNAscope The RNAscope Fluorescent Multiplex Assay (ACD Biosystems) was performed according to RNAscope Multiplex Fluorescent Reagent Kit v2 user manual for fresh-frozen tissue samples. Briefly, 14 μm fresh frozen sections from T13 and L1 DRGs with each side containing retrogradely-labeled tdTomato + LN- or skin- innervating neurons from the same animal were hybridized with RNAscope probes for Ptgir (487851), tdTomato (317041-C2), and Prokr2 (498431-C3) simultaneously. The probes were amplified and detected with TSA plus fluorescein, cyanine 3 and cyanine 5 (Perkin Elmer). The ACD 3-plex negative control probe was run in parallel on separate sections in each experiment to assess the background level and set the acquisition parameter. All sections with tdTomato + cells were imaged at 20x using an Olympus FV3000 confocal imaging system. The frequency of Ptgir + or Prokr2 + DRG neurons among the tdTomato + LN- or skin-innervating population was determined by considering all the tdTomato + cells that were recovered and uniquely-defined from a single animal. Tamoxifen treatment Tamoxifen (Sigma) was dissolved in corn oil (Sigma) at a concentration of 20 mg/mL by shaking overnight at 37°C, and stored at 4°C for the duration of the injections. For labeling arterial vessels with Bmx-CreER T2 , 0.5 mg of tamoxifen was delivered intraperitoneally to label arterial ECs into Bmx-CreER T2 ; Rosa26 eYFP/+ animals between 4–6 weeks of age daily for three consecutive days. Animals were analyzed between 1–3 weeks later.

6-OHDA treatment

For sympathetic denervation, the stock solution of 6-hydroxydopamine (6-OHDA) (Sigma) was prepared in water at 42 mg/mL and stored at −20°C. Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ animals from the same litter between the ages of 6–12 weeks were injected intraperitoneally with 6-OHDA (100 mg kg −1 ) or an equal volume of saline daily for 5 consecutive days. Animals were analyzed the following day.

Optogenetic stimulation of iLN-innervating sensory neurons

Age-matched adult Nav1.8 Cre/+ ; Rosa26 LSL-ChR2-eYFP/+ (ChR2+) or Nav1.8 Cre/+ ; Rosa26 LSL-eYFP/+ (ChR2−) animals (6–12 weeks) were deeply anesthetized (isoflurane, 1.5%–2%, Patterson Vet) maintained at normal body temperature with a water bath heating system (Baxter) during surgical preparation and photostimulation. The animals were surgically prepared for intravital optogenetic stimulation using a method that was adapted from a previously-described protocol for intravital microcopy of iLNs ( von Andrian, 1996 ). Briefly, the skin with the left iLN was flipped inside out following a small incision immediately left to the midline and glued onto a metal block to keep the medulla side of LN exposed. Care was taken not to overstretch the skin flap and damage lymphatic and blood vessels. The site of illumination, the branch point of the antero-posterior-running segment of the y-shaped superficial epigastric artery from where LN feeding arterioles emerged was located and exposed with microdissection without compromising the blood vessel integrity while the tissue was kept moist with normal saline. The stimulation chamber was then built around the iLN with vacuum grease on the side to keep solution from leaking, as well as a metal hairpin shaped tubing with hot water flowing inside on top of vacuum grease to maintain the tissue between 36.5 and 37°C. A thermocouple was placed next to the branch point to monitor the temperature at the tissue. An optic fiber (200 μm core, Thorlabs) coupled to a DPSS laser light source (473 nm, Shanghai Laser & Optics Century) was positioned for focal illumination directly on top of the branch point. The stimulation chamber was subsequently filled to the metal tubing with GenTeal Tears Lubricant Eye Gel (Alcon) to keep the tissue from drying out during stimulation. Pulsed light stimulation (5 ms pulses, 125 mW/mm 2 intensity, 20 Hz) was delivered to the targeted region for 3 hr under the control of a shutter system (Uniblitz). iLNs from both sides were immediately removed after light stimulation and kept in ice cold LN media (HBSS (Corning)/2% FBS/10 mM HEPES/2 mM CaCl 2 ) until subsequent processing.

LN Dissociation and Single Cell Isolation

LNs were kept on ice until processing, < 60 minutes between animal sacrifice and tissue digestion. To minimize biases introduced during tissue dissociation, a gentle dissociation protocol optimized for reliable isolation of both stromal and hemopoietic LN cells was used to generate a suspension of single cells from both the non-immune and the immune compartments ( Fletcher et al., 2011 ). LN media was aspirated, and each LN was placed in 1 mL of pre-warmed digestion media (0.8 mg/mL dispase, 0.2 mg/mL collagenase P (Roche), 50 μg/mL of DNase I in LN media). Using a pair of needle-nose forceps, the capsule of each LN was gently pierced, and the LN in digestion media were placed in a 37°C water bath for 20 minutes with no agitation. Next, LNs were gently agitated without touching the tissue, pelleted by gravity, and the 1 mL of digestion media supernatant was removed and placed in a collection tube on ice containing 10 mL of quenching buffer (PBS/5 mM EDTA/5% FBS). A fresh 1 mL of digestion buffer was added to each LN, and the LNs were placed back in the 37°C water bath for an additional 5 minutes. The LN was gently agitated and triturated using a 1000 μL pipette tip, solid capsular and stromal matter was allowed to settle to the bottom of the tube without centrifuging, and the supernatant digestion media was added to the same collection tube containing quenching buffer. 5-minute incubation periods in fresh digestion buffer and trituration with a 1000 μL pipette tip continued until LNs were completely digested, typically requiring 3–4 additional digestion steps. The cellular suspension in quenching buffer was filtered through a 100 μm filter, and washed with an additional 15 mL of quenching buffer. Single-cell suspensions were centrifuged at 300g for 3 minutes at 4°C, and counted using a hemocytometer and light microscope. We recovered an average of 4.00 +/− 0.53 million cells per LN, and observed no differences in cellularity by treatment group or animal genotype. We saved an aliquot of 60,000 cells from each sample in quenching media on ice as the unenriched sample, and centrifuged the remaining cells at 300 g for 3 minutes at 4°C. Next, using the Miltenyi CD3ε microbead kit and CD19 mouse microbead kit, all remaining LN cells were stained according to manufacturer instructions with the following modifications. First, single cells were stained with CD3ε biotin for 10 minutes on ice, washed once with MACS buffer (PBS/0.5% BSA (Sigma)/2 mM EDTA) and stained simultaneously with CD19 microbeads and biotin microbeads. Cells were isolated using LD columns (Miltenyi) according to manufacturer specifications and the flow-through was collected as the non-T and non-B enriched sample. Single cells from both enriched and unenriched samples were pelleted by centrifugation at 300g for 3 minutes at 4°C, and counted using a hemocytometer with trypan blue staining to estimate cell viability. Across 14 LNs, we recovered an average of 270,000 +/− 31,000 (mean +/− SEM) cells per lymph node following CD3ε and CD19 depletion with > 90% viability. For LN cellularity analysis, single-cell suspensions of the two iLNs from the same ChR2+ or ChR2− mouse (6–12 weeks) were prepared as above. The cells were then filtered through steel mesh and resuspended at the appropriate cell density in FACS buffer before being acquired on a BD Accuri ™ C6 Plus flow cytometer (BD Biosciences). LN scRNA-seq using Seq-Well Single cells from each lymph node prior to and post CD3ε and CD19 depletion were kept separate and diluted to 15,000 cells in 200 μL complete media (RPMI 1640/10% FBS). Seq-Well was performed as described with changes noted below ( Aicher et al., 2019 ; Gierahn et al., 2017 ). Briefly, a pre-functionalized PDMS array containing ~86,000 nanowells was loaded with uniquely-barcoded mRNA capture beads (ChemGenes) ( Macosko et al., 2015 ) and suspended in complete media for at least 20 minutes. 15,000 cells were deposited onto the top of each PDMS array and let settle by gravity into distinct wells. The array was gently washed with PBS, and sealed using a functionalized polycarbonate membrane with a pore size of 0.01 μm, which allows exchange of buffers without permitting mixing of cell materials between different wells. Seq-Well arrays were sealed in a dry 37°C oven for 40 minutes, and submerged in a lysis buffer containing 5 M guanidium thiocyanate (Sigma), 1 mM EDTA, 1% beta-mercaptoethanol and 0.05% sarkosyl (Sigma) for 20 minutes at room temperature. Arrays were transferred to hybridization buffer containing 2 M NaCl (Fisher Scientific) with 8% (v/v) polyethylene glycol (PEG, Sigma) and agitated for 40 minutes at room temperature, mRNA capture beads with mRNA hybridized were collected from each Seq-Well array, and beads were resuspended in a master mix for reverse transcription containing Maxima H Minus Reverse Transcriptase and buffer, dNTPs, RNase inhibitor, a 5’ template switch oligonucleotide (Seq-Well 5’ TSO), and PEG for 30 minutes at room temperature, and overnight at 52°C with end-over-end rotation. Exonuclease digestion was carried out as described previously: beads were washed with TE with 0.01% tween-20 (Fisher Scientific) and TE with 0.5% SDS (Sigma), denatured while rotating for 5 minutes in 0.2 mM NaOH, and resuspended in ExoI (NEB) for 1 hour at 37°C with end-over-end rotation ( Hughes et al., 2019 ). Next, beads were washed with TE + 0.01% tween-20, and second strand synthesis was carried out by resuspending beads in a master mix containing Klenow Fragment (NEB), dNTPs, PEG, and the dN-SMRT oligonucleotide (Seq-Well Second Strand Primer) to enable random priming off of the beads. PCR was carried out as described using 2X KAPA HiFi Hotstart Readymix and ISPCR primer (Seq-Well ISPCR), and placed on a thermal cycler using the following protocol: 95°C for 3 minutes, followed by 4 cycles of 98°C for 20 seconds, 65°C for 45 seconds, 72°C for 3 minutes, followed by 12 cycles of 98°C for 20 seconds, 67°C for 20 seconds, 72°C for 3 minutes, followed by a final 5-minute extension at 72°C. Post-whole transcriptome amplification proceeded as described above for SMART-seq2 libraries, with the following exceptions: AMPure XP SPRI bead cleanup occurred first at a 0.6 × volume ratio, followed by 0.8x. Library size was analyzed using an Agilent Tapestation hsD5000 kit, confirming the expected peak at ~1000 bp, and absence of smaller peaks corresponding to primer. Libraries were quantified using Qubit High-Sensitivity DNA kit and prepared for Illumina sequencing using Nextera XT DNA Sample Preparation kit using 900 pg of cDNA library as input to tagmentation reactions. Amplified final libraries were purified twice with AMPure XP SPRI beads as before, with a volume ratio of 0.6x followed by 0.8x. Libraries from 3 Seq-Well arrays were pooled and sequenced together using a NextSeq 500/550 High Output v2 kit (75 cycles) using a paired end read structure with custom read 1 primer (Seq-Well CR1P): read 1: 20 bases, read 2: 50 bases, read 1 index: 8 bases.

Supplementary Material Table S1 Table S1. Differentially expressed genes by Neuron Type and between LN-innervating and skin-innervating neurons, gene ontology lists. Related to Figures 4 and 5 Table S2 Table S2 . Steady state LN cell atlas: Cell type cluster defining genes and subcluster defining genes. Related to Figure 6 Table S3 Table S3 . Differentially expressed genes following neuronal stimulation by cell type and genotype. Related to Figure 7 Movie 3 Movie 3. A. 0’00” – 0’21”: 3D animation of whole-mount popLN from Nav1.8 Cre ; Rosa26 LSL-tdTomato/LSL-tdTomato animals, stained for tdTomato (green), SMA (magenta), CD31 (cyan). B. 0’22” – 1’19” : 3D animation and fly through view of whole-mount iLN from Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/LSL-tdTomato animals, stained for tdTomato (red), PNAd (green) and SMA (blue), demonstrating a general lack of contact between tdTomato + sensory fibers and PNAd + high endothelial venules (HEVs). C. 1’20” – 1’57” : 3D animation and fly through view of whole-mount popLN from Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ ; Prox1-EGFP animals, stained for tdTomato, GFP, and collagen type 1 (blue). Capsular/subcapsular plexus of sensory nerves (white); parenchymal sensory fibers (red). Related to Figures 2E , S2B , and 2F . Movie 2 Movie 2. A. 0’00” – 0’17”: 3D animation of tdTomato + sensory fibers within popLN of Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ ; Prox1-EGFP animals color-coded by penetration depth (based on the outermost layer of GFP + LECs). B. 0’18” – 0’39” : Slice view of whole-mount popLN from Nav1.8 Cre/+ ; Rosa26 LSL-tdTomato/+ ; Prox1-EGFP animals, stained for tdTomato (white), LYVE-1 (red) and CD45 (blue). Related to Figures 2A and 2C . Fig S4 Figure S4. Contextualizing LN- and skin-innervating sensory and evaluation of individual gene and gene pathway expression Related to Figures 4 and 5 A. Principal Components 2 vs. 4 of Usoskin, Furlan Atlas ( Usoskin et al., 2015 ), represented by transparent circles, colored by previously-defined cell types: non-peptidergic nociceptors (NP, light blue); peptidergic nociceptors (PEP, orange); neurofilament containing (NF, red); and tyrosine hydroxylase containing (TH, purple). LN-innervating (yellow squares) and skin-innervating (blue squares) neurons are projected onto the PC space (see STAR Methods ). B . Euclidean distance between each LN-innervating neuron (left) or skin-innervating neuron (right) and neurons in the Usoskin, Furlan Atlas, separated by cell type. Dashed lines represent the 99% confidence interval for distance between single cells categorized as the same cell type within the Atlas. Box represents 25–75 quantiles, error bars span min-max range. C. Spearman correlation between the scRNA-seq profiles of LN- or skin-innervating neurons and neuronal subsets from the Usoskin, Furlan Atlas. Hierarchical clustering divides LN- and skin-innervating neurons into 4 major subtypes: Neuron Type 1 (PEP1-like, black); Neuron Type 2 (NP-like, tan); Neuron Type 3 (mixed PEP2/NF123, turquoise); and Neuron Type 4 (mixed PEP2/NF12345, dark blue). D. and E. Expression of Prokr2 ( D ) and Ptgir ( E ) by innervation target and Neuron Type. Significance assessed by Mann-Whitney-Wilcoxon test * p