Peripheral Nerve Single-Cell Analysis Identifies Mesenchymal Ligands that Promote Axonal Growth

Abstract Peripheral nerves provide a supportive growth environment for developing and regenerating axons and are essential for maintenance and repair of many non-neural tissues. This capacity has largely been ascribed to paracrine factors secreted by nerve-resident Schwann cells. Here, we used single-cell transcriptional profiling to identify ligands made by different injured rodent nerve cell types and have combined this with cell-surface mass spectrometry to computationally model potential paracrine interactions with peripheral neurons. These analyses show that peripheral nerves make many ligands predicted to act on peripheral and CNS neurons, including known and previously uncharacterized ligands. While Schwann cells are an important ligand source within injured nerves, more than half of the predicted ligands are made by nerve-resident mesenchymal cells, including the endoneurial cells most closely associated with peripheral axons. At least three of these mesenchymal ligands, ANGPT1, CCL11, and VEGFC, promote growth when locally applied on sympathetic axons. These data therefore identify an unexpected paracrine role for nerve mesenchymal cells and suggest that multiple cell types contribute to creating a highly pro-growth environment for peripheral axons.


Introduction
Following injury, mammalian peripheral neurons can regenerate and reinnervate their target tissues. Their ability to do so is thought to be a consequence of a peripheral nerve environment that is highly supportive of axonal growth. Support for this idea comes from classic studies with CNS neurons, which normally fail to regenerate following brain or spinal cord injury but will regrow their axons when peripheral nerve segments are transplanted into the damaged region (David and Aguayo, 1981; for review, see Benowitz et al., 2017). Intriguingly, peripheral nerves are also important for maintenance, repair and regeneration of the non-neural tissues that they innervate. For example, normal peripheral innervation is essential for mammalian hair follicle and hematopoietic stem cells (Brownell et al., 2011;Yamazaki et al., 2011), for cardiac and dermal repair (Mahmoud et al., 2015;Johnston et al., 2013Johnston et al., , 2016 and for amphibian limb (for review, see Kumar and Brockes, 2012) and murine digit tip regeneration (Johnston et al., 2016).
The supportive peripheral nerve environment has largely been ascribed to growth factors made by nerve cells (for review, see Terenghi, 1999;Fledrich et al., 2019). These nerve-derived ligands have been particularly well studied with regard to axonal development and regeneration (Chen et al., 2007;Fledrich et al., 2019), although several studies have shown that they are also important for limb and digit tip regeneration (Kumar et al., 2007;Johnston et al., 2016). These growth factors are thought to be Schwann cell derived, since transplantation of Schwann cells alone is enough to promote CNS axon regeneration (for review, see Bunge, 2016) and murine digit tip regeneration (Johnston et al., 2016). In addition to growth factors, the peripheral nerve provides an extracellular matrix environment that is highly conducive to axonal growth, particularly by contrast to the CNS, where known axon growth inhibitors prevail (Chen et al., 2007;Benowitz et al., 2017). This supportive substrate is also thought to derive in part from Schwann cells, which generate a basal lamina and synthesize ECM proteins and cell adhesion molecules (Muir, 2010;Gardiner, 2011).
These studies all indicate that Schwann cells play an important role in establishing a nerve environment that is supportive of axonal growth. However, the nerve is a structurally-complex tissue containing many different cell types, including vasculature-associated cells, immune cells such as tissue-resident macrophages, and mesenchymal cells of both mesodermal and neural crest origin. In this regard, one recent study identified four transcriptionally and spatially-distinct populations of Pdgfra-positive mesenchymal cells within the injured peripheral nerve, including endoneurial mesenchymal cells that are tightly associated with Schwann cells and axons (Carr et al., 2019). These nerve mesenchymal cells were shown to directly contribute to the repair and regeneration of mesenchymal target tissues including the digit tip, bone, and dermis. Nerve mesenchymal cells have also been shown to play an essential role in forming bridges over gaps in injured nerves (for review, see Cattin and Lloyd, 2016). Together, these findings raise the possibility that mesenchymal cells might also be important for axonal growth in the peripheral nerve.
Here, we provide support for this concept, using an unbiased systems biology approach to define the sciatic nerve ligand environment. We show, using single-cell profiling, that under both homeostatic and injury conditions, mesenchymal cells and Schwann cells are the predominant sources of peripheral nerve ligands, including known and uncharacterized ligands, and that there is induction of ligand expression in both these cell types following injury. Moreover, using mass spectrometry, transcriptional profiling, and computational modeling, we show that peripheral neurons and CNS retinal ganglion neurons express receptors for many of these ligands. Finally, we validate three of these ligands, ANGPT1, CCL11, and VEGFC, as being synthesized and secreted by Pdgfra-positive nerve mesenchymal cells and show that they can promote growth when applied to axons of peripheral sympathetic neurons. Thus, our data support a model where nerve mesenchymal cells and Schwann cells collaborate to establish a generally supportive growth environment in the peripheral nerve.

Animals
All animal procedures were performed in accordance with Canadian Council on Animal Care regulations as approved by the Hospital for Sick Children animal care committee. Sprague Dawley rats (purchased from Charles River) used in this study ranged from embryonic day (E)15 to young adult (six weeks old) and CD1 mice (purchased from Charles River) ranged in age from eight to twelve weeks old. All rats and mice were healthy throughout the duration of the study and had free access to chow and water in a 12/12 h light/dark cycle room. In most cases, rats and mice of both sexes were used with the exception of sixweek-old male rats for sciatic nerve injury microarray experiments. Pdgfra EGFP/1 (B6.129S4-Pdgfra tm11(EGFP)Sor /J; JAX stock #007669; Hamilton et al., 2003) mice were obtained from The Jackson Laboratory and were bred and genotyped as recommended by The Jackson Laboratory. Animals that underwent sciatic nerve injury surgeries were housed individually for recovery purposes.

Sciatic nerve resection surgeries
Sciatic nerve resections were performed on young adult male Sprague Dawley rats (microarray analysis), adult CD1 mice (scRNA-seq analysis) or adult Pdgfra EGFP/1 mice [fluorescence in situ hybridization (FISH) and immunostaining]. Before surgery, animals were anesthetized with 2% isoflurane gas and the surgical site was shaved. Animals were kept under anesthesia for the duration of the surgery. To resect the sciatic nerve, an incision was made along the lateral aspect of the mid-thigh of the right hindlimb, the sciatic nerve was then raised, an ;5to 10- Figure 1. Characterization of ligand expression in the injured sciatic nerve (see also Extended Data Fig. 1-1). A, Images of longitudinal sections of an uninjured adult nerve and a 9 DPI distal sciatic nerve from Pdgfra Egfp/1 mice analyzed for EGFP (green) and immunoreactivity for PDGFRa (red) and S100b (white). Arrowheads denote endoneurial cells positive for both PDGFRa protein and nuclear Pdgfra-EGFP and arrows indicate S100b immunoreactive Schwann cells. Scale bars = 100 mm. B, Venn diagram showing the number of ligands expressed in the uninjured versus 3 and 7 DPI distal sciatic nerves, based on microarray analysis. Ligand mRNAs were defined as expressed if their levels were !Ntf3. C-I, Characterization of ligand expression in injured distal sciatic nerve scRNA-seq datasets. C, t-SNE cluster visualization of 3 DPI sciatic nerve cell transcriptomes analyzed via the computational mm segment was removed, and the distal nerve ending was carefully tucked away (distally) from the injury site to prevent regeneration. The wound was then closed with 4-0 Polysorb sutures (Covidien). Animals were treated subcutaneously with ketoprofen or meloxicam (;2-5 mg/kg) as well as buprenorphine (0.05 mg/kg) before surgery, along with a postoperative treatment of ketoprofen or meloxicam 24 h after surgery. Mice and rats were housed separately following surgery and remained healthy throughout the postoperative period and were monitored twice daily for 3 d following surgery.

Single-cell isolation and myelin removal for Drop-seq analysis
For preparation of the 3 d postinjury (DPI) nerve scRNAseq dataset, young adult CD1 mice underwent unilateral surgical resections as described above, and injured distal sciatic nerve segments were collected 3 d following surgery. For the uninjured nerve and neonatal nerve analyses, bilateral sciatic nerve segments were collected from adult and postnatal day (P)2-P4 CD1 mice, respectively. Freshly dissected nerves were digested in a mixture of collagenase Type XI (1 mg/ml, Sigma) and 0.05% Trypsin-EDTA (Thermo Fisher Scientific) for 30 min at 37°C. Enzymatic digestion was halted by diluting the cell suspension with HBSS (Thermo Fisher Scientific). Following centrifugation (1200 rpm for 5 min) and removal of the supernatant, the cell pellet was resuspended in PBS containing 0.5% BSA and passed through a 70-mm cell strainer (BD Biosciences). For datasets purified with myelin removal beads (3 DPI, neonatal and uninjured nerve; as shown in Figs. 1C,E, 2C,E; referred to as set 2 for the neonatal analyses, where cells were prepared in two ways), myelin debris was removed from the single-cell suspension using Myelin Removal Beads II and a MidiMACS magnetic separator with LS columns (Miltenyi Biotec), according to the manufacturer's instructions. Following myelin removal, the cell suspension was centrifuged (1200 rpm for 5 min), and the supernatant was removed before resuspending the pellet in 0.22-mm sterile-filtered PBS containing 0.01% BSA. For the second neonatal nerve dataset that was purified using fluorescence-activated cell sorting (FACS), a single-cell suspension of dissociated injured nerve cells was prepared as described above. After passing the cells through a 70-mm cell strainer and resuspending them in PBS containing 0.25% BSA, Hoechst 33258 was added to distinguish nucleated cells from myelin debris, in addition to propidium iodide (PI) to exclude dead cells. The Hoechst high and PI-negative cell fractions were FACS purified using a MoFlo XDP cell sorter (Beckmann Coulter) before proceeding with scRNA-seq analysis. In all cases, cells were then resuspended in PBS containing 0.01% BSA, counted with a hemocytometer, and the solution was adjusted to a final concentration of 140,000 cells/ml and run through the Drop-seq apparatus at the Princess Margaret Genomics Facility. Drop-seq, cDNA amplification, library preparation, sequencing, processing of FASTQ sequencing reads, and read alignment steps were all conducted including minor modifications according to previously published protocols (Macosko et al., 2015). For the 3 DPI nerve scRNA-seq analysis (as shown in Fig. 1C), a raw digital gene expression (DGE) matrix was generated from 2500 cell barcodes as described in the Drop-seq Alignment Cookbook (version 1.2, January 2016; http://mccarrolllab. com/dropseq/). Similarly, for the uninjured nerve scRNA-seq analysis (as shown in Fig. 2C), a raw DGE matrix was generated from 2000 cell barcodes and used for all further analyses. In the case of the two neonatal nerve datasets (FACS sorted and bead treated), 2500 and 6200 cell barcodes were used to generate the DGE matrices as described above. DGE matrices described here were used for all subsequent analyses. The previously published DGE matrices for the 9-d injured nerve datasets (both FACS and myelin bead treated; GEO:GSE120678) were described in Carr et al. (2019).

Computational analysis of scRNA-seq data
Drop-seq data (DGE matrices) were analyzed used a previously described custom computational pipeline Carr et al., 2019;Storer et al., 2020; described in detail in Innes and Bader, 2019). Briefly, data were filtered to remove cells with low unique molecular identifier counts, cell doublets, contaminant red blood cells, and cells that contained high mitochondrial gene content. Genes detected in less than three cells were removed. Cell transcriptomes were then normalized as previously described (Lun et al., 2016) using an algorithm in the scran package in R that corrects for differences in sequencing depth by the use of scaling factors within each cell by pooling random subsets of cells, summing their library sizes, and comparing them with average library size across all cells in the group. This is iteratively performed, and the cell-wise scaling factors can be deconvolved from the set of pool-wise scaling factors. Following continued pipeline, with clusters annotated for cell types as identified by marker gene expression. D, t-SNE gene expression overlays on the dataset in C for the Schwann cell marker Sox10, the macrophage marker Aif1, and the mesenchymal cell marker Pdgfra. Relative transcript expression levels are color coded as per the adjacent color keys. E, t-SNE cluster visualization of the combined 3 and 9 DPI distal sciatic nerve cell transcriptomes with clusters annotated for cell types as identified by marker gene expression. F, t-SNE visualization of the dataset in E with cells color coded for their dataset of origin. Numbers correspond to cluster numbers in E. G, Bar graph showing the percentage of the 143 injured nerve ligand mRNAs detectably expressed in the combined 3 and 9 DPI sciatic nerve cell types (shown and annotated in E), including Pdgfra-positive mesenchymal cells (Comb. Mes), Pdgfra-positive endoneurial mesenchymal cells (Endo. Mes), endothelial cells (Endoth.), Schwann cells (SC), VSM/pericyte cells (VSM/Peri.), and immune cells. Ligand mRNAs were considered to be expressed if they were detected in 2% or more of the cell type of interest. H, Stacked bar graph showing the relative percentage of ligand mRNAs expressed at the highest levels in the different injured peripheral nerve cell types shown in E. I, t-SNE gene expression overlays of the combined 3 and 9 DPI sciatic nerve dataset (shown in E) for Fgf7, Il33, Bmp7, and Wnt5a. Cells that detectably express the ligand are colored blue and the numbers correspond to the clusters.
Research Article: New Research normalization, DGE matrices were imported into Seurat (v. 1.4.0.16) in R. Principal component (PC) analysis was then undertaken using highly variable genes, and clustering analysis was performed using top PCs. This analysis was conducted using the shared nearest neighbor (SNN) modularity optimization-based clustering algorithm implemented in Seurat (FindClusters function). Clustering was iteratively performed at increasing resolution until a lower limit of ;10-20 differentially expressed genes [calculated by the Seurat FindMarkers function, p , 0.01 family-wise error rate (FWER), Holm's method] was reached between the most similar clusters. For conservative analysis of all datasets analyzed, clusters were assigned at the lowest resolution that still distinguished distinct cell types, as defined by established marker genes. As a result, clusters were assigned at a resolution of 0.4 for analysis of all datasets.
For the 3 DPI dataset (2075 total cells), nine clusters were identified with 210 differentially expressed genes between most similar clusters (p , 0.01, FWER). Cells were sequenced to an average depth of .70,000 reads/cell. The average number of genes detected per cell was 1027 6 588, and the average number of transcripts was 2257 6 2855. For the 3 and 9 DPI nerve combined dataset, cell transcriptomes from the 3 DPI dataset and the myelin bead removal-treated 9 DPI dataset (from Carr et al., 2019;GEO:GSE120678) were merged using the unique cell identifier barcodes from all cells present in all clusters of the two datasets following pipeline processing. The constructed raw DGE matrices of the combined datasets were then re-run through the pipeline for re-clustering, resulting in 5395 total cells. Twelve clusters were identified with 200 differentially expressed genes between most similar clusters (p , 0.01, FWER). The uninjured nerve dataset was previously analyzed for Pdgfra-positive mesenchymal cells (Carr et al., 2019), but not for other cell types. Reanalysis of this dataset (1841 total cells) identified 11 clusters with 105 differentially expressed genes between the most similar clusters (p , 0.01, FWER). For the combined neonatal dataset, cell transcriptomes from both the FAC-sorted (set 1; Extended Data Fig. 2-1C) and myelin removal bead-treated (set 2; Extended Data Fig. 2-1C) samples were merged using the unique cell identifier barcodes from all cells present in all clusters of the two datasets following pipeline processing. The constructed raw DGE matrices of the combined datasets were then rerun through the pipeline for re-clustering, resulting in 6885 total cells. Ten clusters were identified with 540 differentially expressed genes between most similar clusters (p , 0.01, FWER). For the FAC-sorted sample, cells were sequenced to an average depth of .90,000 reads/cell. The average number of genes detected per cell was 1005 6 728 and the average number of transcripts was 2044 6 1985. For the bead treated sample, cells were sequenced to an average depth of .43,000 reads/cell. The average number of genes detected per cell was 732 6 485, and the average number of transcripts was 1231 6 1095.
For the combined Schwann cell dataset (Fig. 3A,B), the unique cell identifier barcodes from all cells present in Sox10-positive clusters in each of the six described datasets (Figs. 1E, 2C,E; FAC-sorted cells in Carr et al., 2019) were merged. For the combined mesenchymal cell dataset ( Fig. 4A,B), the unique cell identifier barcodes from all cells present in Pdgfra-positive clusters in each of the same six datasets were merged. The constructed raw DGE matrices of the combined datasets were then re-run through the pipeline, resulting in 5331 total Schwann cells and 5416 total mesenchymal cells. Batch correction of these combined datasets as well as data shown in Extended Data Figure 2-1C was performed using the Harmony batch-effect-correction method (Korsunsky et al., 2019) with Seurat V2. Briefly, gene expression data from the combined datasets were transferred to Seurat, where highly variable genes were then used to carry out principal component analysis. The Harmony iterative algorithm was used to integrate datasets and adjust for dataset specific effects based on the top 20 principal components. Iterative clustering was performed using the FindClusters function in Seurat V2, with clusters being assigned at a resolution of 0.4. This resulted in seven and nine clusters in the combined Schwann cell and mesenchymal cell datasets, respectively. t-Distributed stochastic neighbor embedding (t-SNE) visualizations of batch-  In combined datasets, dataset identities were distinguished by using the gg color hue and hcl functions in R. Correlation analysis comparing gene expression between different clusters was performed by averaging the expression of each gene across all cells in the individual clusters to be compared, then Pearson correlation analysis was performed using the Cor function and plotted in R. Genes of interest were then highlighted using Adobe Illustrator (Adobe Systems Incorporated; as in Extended Data Fig.  3-1B). Differential correlation of single-cell transcriptomes as shown in Figure 3D was performed as described previously (Gerber et al., 2018). Briefly, mock bulk transcriptomes were generated for the 9 DPI Schwann cells (both bead and FAC sorted), the uninjured non-myelinating Schwann cells and the neonatal Schwann cells (both bead and FAC sorted) by determining the mean expression of each gene in the total combined cells in each dataset. Each single-cell transcriptome was then correlated with each of the mock bulk transcriptomes. We then determined the differential correlation of each single cell with the bulk uninjured nerve versus the bulk 9 DPI transcriptomes (y-axis) and the differential correlation with the bulk uninjured nerve versus the bulk neonatal nerve transcriptomes (x-axis). Violin plots were generated using the VioPlot package in R. Hierarchical clustering of the batch corrected combined Schwann cell data in Figure 3C was performed based on the top 20 Principal Component (PC) using the BuildClusterTree and PlotClusterTree functions in Seurat. Node numbers were removed from the plot and cluster descriptions and colors were added for clarity using Adobe Illustrator (Adobe Systems Incorporated). The single cell heatmaps were generated (with scaled expression values) using the DoHeatMap function in Seurat at resolution 0.4.

Ligand mRNA expression and Venn diagram analysis
The expression of ligand mRNAs (Table 1) was characterized from the whole nerve microarray analysis using a custom curated ligand-receptor database (modified from Yuzwa et al., 2016). Extracellular matrix proteins and potential ligands without well-defined, receptor-mediated paracrine roles were excluded. The VennDiagram package in R was used to determine overlapping ligands in the uninjured, 3-d injured, and 7-d injured nerve datasets and was modified to show proportional representation of data using the eulerAPE tool (Micallef and Rodgers, 2014).
The combined 3 and 9 DPI nerve scRNA-seq dataset ( Fig. 1E) was analyzed to identify the percentage of cells in each cell type expressing the ligand mRNAs identified by the microarray analysis. For this analysis, Pdgfra-positive mesenchymal cells were separated into endoneurial cells (cluster 5; Fig. 1E) and all other mesenchymal cells. Ligand mRNAs were considered further only if they were detectable in 2% or more cells of at least one cell type. The 143 resultant injured nerve ligands (Table 2) were further analyzed in the other scRNA-seq datasets. Venn diagrams comparing expression of the 143 injured nerve ligands in the various datasets (Figs. 2G, 3E, 4C) were prepared using the VennDiagram package, modified to continued transcriptomes from the injured, uninjured and neonatal nerve datasets (in Figs. 1E, 2C,E; also see Carr et al., 2019) were extracted, combined together, analyzed and batch-corrected using Harmony data integration and clustered in Seurat based on principal components. A, B, t-SNE visualizations of the combined Schwann cells (SCs), showing clusters (A) and datasets of origin (B). The clusters in A are also annotated based on marker gene expression (Extended Data Fig. 3-1A) and datasets of origin. Neonatal and 9 d refer to the FAC-sorted preparations at these timepoints while neonatalset 2 and 9 dset 2 refer to the cells prepared with myelin removal beads. C, A dendrogram showing hierarchical analysis of cell clusters from the combined Schwann cell dataset, with cluster identity numbers, annotations and colors as in A. D, Scatterplot showing differential correlation of single-cell transcriptomes from the combined Schwann cell dataset in A (individual colors represent different datasets) with bulk transcriptomes from uninjured adult versus neonatal Schwann cells on the x-axis and uninjured adult versus 9 DPI on the y-axis. E, Venn diagram showing the overlap of ligand mRNAs expressed in neonatal, injured and uninjured Schwann cell clusters from the combined Schwann cell dataset in A. Ligand mRNAs were considered to be expressed if they were detectable in 2% or more cells in any defined cell type cluster. F, t-SNE gene expression overlays of the combined Schwann cell data in A for Btc, Bdnf, Fgf5, Ucn2, Clcf1, and Gdnf. Cells that detectably express the ligand are colored blue, and the numbers correspond to the clusters. Injured Schwann cell cluster 3 (Inj.) is circled. G, Violin plots showing the relative expression of Bmp1, Gdnf, Shh, and Pdgfa in injured Schwann cluster 3 versus uninjured Schwann cell cluster 4 from the dataset in A. H, t-SNE gene expression overlays of the combined Schwann cell data in A for Igf1, Dhh, and Fgf1. Cells that detectably express the ligand are colored blue and the numbers correspond to the clusters. I, Violin plots showing the relative expression of Dhh and Fgf1 in injured Schwann cell cluster 3 versus uninjured cluster 4 from the dataset in A.  Fig. 4-1). Pdgfra-positive mesenchymal cell transcriptomes from the injured, uninjured and neonatal nerve datasets (in Figs. 1E, 2C,E; also see Carr et al., 2019) were extracted, combined together, analyzed and batch-corrected using Harmony data integration and clustered in Seurat based on principal components. A, B, t-SNE visualizations of the combined mesenchymal cells showing clusters (A) and datasets of origin (B). The clusters in A are also annotated based on marker gene expression (Extended Data Fig. 4-1) and datasets of origin. Neonatal and 9 d refer to the FACsorted preparations at these timepoints, while neonatalset 2 and 9 dset 2 refer to the cells prepared with myelin removal beads. C, Venn diagram showing the overlap of ligand mRNAs expressed in neonatal (cluster 6), injured (cluster 3), and uninjured (cluster 4) nerve endoneurial mesenchymal cells from the combined mesenchymal nerve dataset in A. Ligand mRNAs were considered to be expressed if they were detectable in 2% or more cells in the relevant cluster. D, t-SNE gene expression overlays of the combined mesenchymal cell data in A for Angpt1, Crlf1, Ngf, Sema7a, and Hgf. Cells that detectably express the ligand are colored blue, and the numbers correspond to the clusters. Relevant clusters are circled and annotated, including uninjured endoneurial (Uninj Endo.), injured endoneurial (Inj Endo.), and injured differentiating (Inj Diff.) mesenchymal cells. E, Violin plots showing the relative expression of Angpt1, Crlf1, Lif, and Sema7a in injured endoneurial mesenchymal cell cluster 3 versus uninjured endoneurial cell cluster 4 from the dataset in A. F, t-SNE gene expression overlays of the combined mesenchymal cell data in A for Ccl11, Il33, and Wnt5a. Cells that detectably express the ligand are colored blue, and the numbers correspond to the clusters. G, Violin plots showing the relative expression of Ccl11 and Il33 in injured endoneurial mesenchymal cell cluster 3 versus uninjured endoneurial cell cluster 4 from the dataset in A.
show proportional representation with the eulerAPE tool (Micallef and Rodgers, 2014).

RNA isolation and microarray analysis
Total RNA was extracted from E15 rat dorsal root ganglion (DRG) sensory neurons and neonatal rat superior cervical ganglion (SCG) sympathetic neurons using the RNeasy Micro kit (QIAGEN) according to manufacturer's instructions. RNA was isolated from the distal rat sciatic nerve following injury (3 and 7 d after resection) and contralateral uninjured sciatic nerve as follows. After harvesting, nerves were flash frozen in liquid nitrogen and stored at À80°C until RNA isolation was performed. Nerve tissue was lysed using a Dounce homogenizer with cooled TRIzol (1 ml/50-100 mg of tissue) followed by passing the lysate through a 23.5-gauge needle. The homogenate was then spun at 12,000 Â g at 4°C for 10 min, with the clear supernatant transferred to a new tube and allowed to incubate at room temperature for 5 min, 0.2 ml of chloroform was then added per milliliter of homogenate, shaken vigorously for 15 s, incubated 2-3 min at room temperature, and spun at 12,000 Â g for 15 min at 4°C. The resulting aqueous phase was removed and added to a new tube with equal volume of 70% ethanol before RNA isolation with the RNeasy Micro kit (QIAGEN) according to manufacturer's instructions. Microarray analysis was then performed at the Center for Applied Genomics at the Hospital for Sick Children (Toronto, ON). A total of 250 ng of total RNA was processed using the Affymetrix WT Plus kit to generate cDNA following Bioanalyzer analysis to confirm the quality of the RNA. 5.5 mg of labeled cDNA was hybridized onto Rat Gene 2.0 ST arrays using the Affymetrix FS450_0002 hybridization protocol and scanned using the Affymetrix GeneChip Scanner 3000.

Normalization and differential gene expression analysis of microarray data
Raw probe intensity values were background corrected and then normalized with quantile normalization. Values were transformed into the log2 scale and summarized into probesets using the robust multichip analysis (RMA) algorithm in the oligo Bioconductor package in R. For all of the rat microarray datasets, gene annotation was performed using the "ragene20sttranscriptcluster.db" library in R. Motor neuron microarray data (Kaplan et al., 2014) was obtained from mouse P7 lumbar motor neurons from the GEO database under the accession number GSE52118. For these data, raw probe intensity values were normalized and summarized into probesets as described above, except for the Affymetrix Mouse Genome 430 2.0 Array (n = 3 replicates), and gene annotation was performed using the "mouse4302.db" library in R. The Limma Bioconductor package was used to calculate differential gene expression between the sensory neurons (DRGs, n = 6 replicates) and sympathetic neurons (SCGs, n = 4 replicates). Bayesian statistics were calculated and annotated receptor genes were considered to be differentially expressed if they were !2-fold different with p , 0.05 false discovery rate (FDR; Benjamini and Hochberg correction).

Cell-surface capture mass spectrometry
Cell-surface mass spectrometry was conducted based on modified published protocols (McDonald et al., 2009;Schiess et al., 2009;Yuzwa et al., 2016). Briefly, 6-cm dishes containing sensory and sympathetic neurons were washed with coupling buffer [1Â PBS (pH 6.5) and 0.1% fetal bovine serum (FBS)] following removal of the medium. Cultures were treated with 5 mM NaIO 4 in coupling buffer for 30 min at room temperature in the dark and lysed in buffer containing 20 mM Tris-HCl, 150 mM NaCl, 0.0002% NaN 3 , 1% NP-40 (pH 7.5), and one complete mini-protease inhibitor tablet (Roche) per 10 ml of lysis buffer. Lysates were passed through 23-gauge needles, protein concentrations determined using the Pierce BCA assay kit (catalog #23227, Thermo Fisher Scientific), and equal amounts of total protein (between 0.9 and 1.5 mg for sensory neurons and ,1 mg for sympathetic neurons) were added to 200 ml of Ultralink Hydrazide resin (Pierce) pre-equilibrated with lysis buffer and rotated overnight at room temperature. The unbound protein was removed via centrifugation the following day, washed twice with 8 M urea and three times with 50 mM ammonium bicarbonate (pH 8; ABC). The resin was treated with 50 mM dithiothreitol (DTT; American Bioanalytical) in ABC at 37°C for 60 min, washed once with ABC, and incubated with 65 mM iodacetamide in ABC at room temperature in the dark for 30 min. The resin was washed once with ABC, once with 1.5 M NaCl and three times with ABC and incubated with 40 ng/ml of Trypsin (Worthington) in ABC overnight at 37°C. The following day, the resin was washed three times with 1.5 M NaCl, followed by 80% acetonitrile, methanol, water, and ABC and then incubated with 1300 U/ml of PNGaseF (New England Biolabs) at 37°C overnight in ABC. The following day, the eluted peptides were collected from the resin and the resin washed once with ABC and combined with the eluted peptides. The eluates were lyophilized overnight and prepared for mass spectrometry using C18 reverse-phase ZipTips (EMD Millipore). Peptides were lyophilized and resuspended in 11 ml of 0.1% formic acid following elution from ZipTips. Samples were then analyzed on an Orbitrap analyzer (Q-Exactive, Thermo Fisher) outfitted with a nano-spray source and EASY-nLC nano-LC system (Thermo Fisher); 5 ml of the resuspended peptide mixtures were loaded onto a 75 mm Â 50 cm PepMax RSLC EASY-Spray column filled with 2 mM C18 beads (Thermo Fisher) at a pressure of 800 Bar. Peptides were eluted over 60 min at a rate of 250 nl/ min using a 0-35% acetonitrile gradient in 0.1% formic acid. Peptides were introduced by nano-electrospray into the Q-Exactive mass spectrometer (Thermo Fisher). The instrument method consisted of one MS full scan (400-1500 m/z) in the Orbitrap mass analyzer with an automatic gain control (AGC) target of 1e6, maximum ion injection time of 120 ms, and a resolution of 70,000 followed by 10 data-dependent MS/MS scans with a resolution of 17,500, an AGC target of 1e6, maximum ion time of 120 ms, and one microscan. The intensity threshold to trigger a MS/MS scan was set to 1.7e4. Fragmentation occurred in the HCD trap with normalized collision energy set to 26. The dynamic exclusion was applied using a setting of 8 s. Peptide and protein identification was performed using PEAKS version 8 software (Bioinformatics Solutions Inc.). Peptides and proteins were identified at the 0-0.1% FDR level in the sensory neuron dataset and at the 0.8-1% FDR level in the sympathetic neuron dataset. With both sensory and sympathetic neuron culture datasets, data were pooled across the three samples, and we included all proteins where at least a single peptide was detected in at least one of the samples.

Identification of receptors based on the microarray and mass spectroscopy data
Receptors were identified in both the microarray and mass spectrometry data using the ligand-receptor database described in Yuzwa et al. (2016), with modifications. In addition, proteins that were classified as receptors using the PANTHER (http://pantherdb.org) protein classification system, but not necessarily present in the ligandreceptor database, were included. Manual curation of receptors with known signaling functions (such as the Plxn receptors) as well as other genes identified as receptors by Gene Ontology (GO) terms were also included as an update to the previously published version of the ligandreceptor database. Receptor classifications (Fig. 5B) were based on GO terms and descriptions provided by websites including GeneCards (http://genecards.org) and UniProt (http://uniprot.org).

Computational modeling and pre-processing of data
To define receptors for generation of ligand-receptor models, gene expression values from the injured nerve and from the motor, sensory, and sympathetic neuron microarrays were first averaged across replicates, and in the case of multiple probes corresponding to the same receptor gene, the highest expressed probe was used (this was following averaging). For rat retinal ganglion cells (RGCs), we analyzed RNA-seq data from Blanco-Suarez et al. (2018) as obtained from the GEO database (accession GEO: GSE108484). Data from control RGCs were used (n = 3) and processed as described above for the microarray data. Expression values with fragments per kilobase of transcript per million mapped reads (FPKM) .1 were considered expressed and included for analysis. For the modeling, receptor mRNAs that had expression values exceeding the thresholds as described in the results were included. These receptors and the 143 injured nerve ligands (Table 2) were then analyzed using a custom Python script ("Cellcellinteractnet," Python version 2.7.6) and custom ligand-receptor database (database modified from  to predict ligand-receptor interactions. Models were all constructed from this information using Cytoscape (3.7.0). The Venn diagram comparing predicted interactions in Figure 7C was prepared using the VennDiagram package, modified to show proportional representation with the eulerAPE tool (Micallef and Rodgers, 2014).

Sympathetic neuron cultures
SCGs of newborn (P1-P2) Sprague Dawley rats were dissected, dissociated and cells were plated at eight ganglia per 3.5-cm tissue culture-treated dish for microarrays, ;1 Â 10 6 cells per 6-cm dish for mass spectrometry, or ;1.4 Â 10 5 cells per 13-mm glass coverslip for immunostaining. Neurons were plated on dishes and coverslips coated with poly-D-lysine and laminin (1 mg/ml laminin; VWR). Neurons were cultured in growth medium composed of UltraCULTURE medium (Lonza), 2 mM L-glutamine (Lonza), 100 U/ml penicillin (Wisent) with 100 mg/ml streptomycin (pen/strep; Wisent), and 50 ng/ml NGF (Cedarlane). Neurons were cultured in growth medium containing 3% heat-inactivated FBS for 3 d with inclusion of 7.2 mM cytosine arabinoside (CA) for the first day, and 3.6 mM CA for the second and third days. Cultures were switched to growth medium alone for an additional 3 d, as described previously (Park et al., 2010;Feinberg et al., 2017).

Sensory neuron cultures
Sensory neurons from E15 rat DRGs were cultured as described (Feinberg et al., 2010(Feinberg et al., , 2017. For immunostaining, neurons were plated at a density of 2 Â 10 5 cells per 13-mm glass coverslip precoated with laminin (VWR), and poly-D-lysine (Sigma-Aldrich). Neurons were plated at a density of 1 Â 10 6 cells per 6-cm dish for proteomics. For microarrays, neurons were plated under both conditions. Initially, cells were plated in a medium #1 containing Neurobasal medium (Invitrogen), GlutaMAX (Invitrogen), 50 ng/ml NGF (Cedarlane), B27 supplement (Invitrogen), and pen/strep. The day after plating, cultures were treated for 2 d in medium #2 composed of Eagles basal medium (Invitrogen), ITS supplement (Sigma), 0.2% BSA (Sigma), 4 mg/ml D-glucose (Sigma), GlutaMAX (Invitrogen), 50 ng/ ml NGF (Cedarlane), and pen/strep with 0.8 mM CA to eliminate mitotic non-neuronal cells. Cells were then treated with another cycle of medium #1 for 2 d followed by a final 2-d cycle of medium #2. After the second CA treatment, cultures were grown in medium #1 for two additional days before cell harvest.
Culturing and sorting mesenchymal cells from the sciatic nerve for ELISA analysis Sciatic nerves were harvested from P4 Sprague Dawley rat pups. Each biological sample included four litters of pups (two litters harvested per day over 2 d). Nerve cells from two litters were dissociated and plated on 6-cm dishes in medium consisting of low-glucose DMEM/F12 (3:1, Invitrogen) and 1% pen/strep (referred to as basal medium) with 10% FBS. In order to isolate PDGFRa-positive mesenchymal cells, 1-2 d after plating, cells were treated with an antibody solution containing basal medium with goat anti-PDGFRa antibody (1:250, R&D Systems, catalog #AF1062) and donkey anti-goat Alexa Fluor 488 secondary (1:500, Thermo Fisher Scientific, catalog #A11055) that was preincubated for 30 min at room temperature. Cells were treated for 1 h at 37°C with the antibody solution before being returned to basal medium with 10% FBS overnight. PDGFRa-positive cells were then isolated using a Mo-Flo XDP sorter (Beckman Coulter). Following the sort, between 4-9 Â 10 5 cells were Figure 5. Sensory and sympathetic neuron receptor expression and modeling of their predicted interactions with injured nerve-derived ligands (see also Extended Data Fig. 5-1). A, Images of E15 rat DRG sensory neurons cultured for 9 d (left panel) and neonatal rat SCG sympathetic neurons cultured for 6 d (right panel), immunostained for b III-tubulin (Tubb3; green). Cells were counterstained with Hoechst 33258 to show cell nuclei (blue). Scale bars = 20 mm. B, Bar graphs showing types of receptor proteins detected on re-plated on 6-cm dishes with basal medium and 10% FBS. Approximately 72 h later, medium containing FBS was removed, cells were washed with HBSS, and basal medium without serum was added for conditioning. Cells were treated with cycles of basal medium for conditioning (24-96 h) and basal medium with FBS for a period of 2.5-4 weeks total. Collected conditioned medium was spun down at 1300 rpm for 5 min to removed cell debris and frozen at À80°C until ELISAs were performed. ELISA assays for rat ANGPT1 (catalog #LS-F10827, LSBio), CCL11 (catalog #LS-F11046, LSBio), and VEGFC (catalog #LS-F5482, LS Bio) were performed and results were analyzed according to manufacturer's instructions.

Sympathetic neuron compartment cultures
Campenot cultures were established on 35-mm collagen-coated dishes as previously reported (Campenot et al., 1991); 20 Â 20 mm Teflon chambers (Tyler Research) were grease-sealed and assembled on plates where the substratum was corrugated with 20 parallel tracks 200 mm wide using a pin rake (CAMP-PR, Tyler Research), forming lanes for outgrowth. Neonatal SCG sympathetic neurons were enzymatically and mechanically dissociated and plated in the central chamber at a density of three to four ganglia per compartment in 10 ng/ml NGF (Cedarlane, catalog #CLMCNET-001.1), methylcellulose Ultraculture media (Lonza) supplemented with 2 mM Lglutamine (Lonza), and antibiotics [100 U/ml penicillin (Wisent) and 100 mg/ml streptomycin (Wisent)]; 10 mM CA with 3% FBS was added to the cell bodies for 2-3 d as the axons were permitted to extend into the adjacent compartments containing 10 ng/ml NGF. After this period, central and side chambers were washed with fresh medium, medium in the cell body chamber was replaced with 10 ng/ml NGF, and medium on the axonal chambers was replaced with 100 ng/ml candidate ligand and varying amounts of NGF (0.5 ng/ml for experimental compartments, 50 ng/ml for maximum NGF control). Ligands included human recombinant ANGPT1 (Peprotech, catalog #130-06-5UG), murine recombinant Eotaxin (CCL11, Peprotech, catalog #250-01-5UG), and human recombinant VEGFC (Peprotech, catalog #100-20CD). Images of axon outgrowth were obtained at 3d post-ligand addition using a Zeiss AxioObserver Z1 microscope and outgrowth was analyzed using ImageJ software (NIH). To quantify the density of axonal growth in these compartments, a line was drawn perpendicular to the axis of the outgrowth within the furthest 1 mm of outgrowth where axons were maximally defasciculated and the number of axons crossing the line were quantified.

Immunostaining, imaging, and analysis of cultures
Immunostaining of cultured cells was performed on glass coverslips and tissue culture dishes. Cells were fixed for ;10 min in 4% paraformaldehyde in PBS solution and washed three times (10 min per wash) in PBS before treatment with the primary antibody solution, which contained the antibody, 10% normal donkey or normal goat blocking serum (Jackson ImmunoResearch), and 0.3% Triton X-100 (Fisher) in PBS. After 1-2 h of incubation at room temperature, cells were washed three times (10 min per wash) with PBS. Cells were then incubated in a secondary antibody solution containing antibodies diluted 1:500 in 0.3% Triton X-100 (Fisher) in PBS. After 1 h, cells were then washed three times (10 min per wash) in PBS with Hoechst 33258 (Sigma) added in one of the washes to label nuclei. For glass slides, coverslips were mounted with PermaFluor mounting media (Thermo Fisher Scientific). For tissue culture dishes, mounting medium was added directly on to the cells, and coverslips were placed on top of this. Cells were imaged using either a Zeiss AxioImager M2 microscope with an X-Cite 120 LED light source and a C11440 Hamamatsu camera using Zen acquisition software in the case of glass coverslips, or with a Zeiss AxioObserver Z1 microscope with Zen software in the case of tissue culture dishes. To quantify cells, multiple regions per glass coverslip or dish were imaged and cells were counted. For purified sensory and sympathetic neuron culture quantifications, S100b -positive and b III-tubulin-positive cells were counted in four regions per coverslip across two coverslips per biological replicate and summed before calculating percent positive immunoreactive cells. Fibronectin-positive cells were counted in three regions over two sympathetic neuron replicates. Total cell counts per biological replicate ranged from ;100 to ;600. In cases where cells in dishes were counted, between three and eight regions per dish were used. Images were then processed using Photoshop software (Adobe Systems Incorporated), where brightness and contrast were edited as appropriate.

FISH
Single molecule FISH was performed on 9 DPI adult Pdgfra EGFP/1 nerve sections with probes targeting Vegfc (catalog #492701), Ccl11 (catalog #464031), and Angpt1 (catalog #449271) mRNAs using the RNAscope kit (Advanced Cell Diagnostics), according to the manufacturer's instructions. Briefly, freshly dissected nerves were fixed overnight in 4% PFA at 4°C and cryopreserved in 30% sucrose overnight (at 4°C) under RNase-free conditions using RNase-free reagents. Nerves were cryosectioned sagittally at a thickness of 14 mm. Sections were washed with ethanol, followed by tissue pretreatment (1:10 dilution) for 20 min, probe hybridization, and signal amplification. Positive signal was identified as red punctate dots. Digital images were acquired with a Quorum Spinning-Disk confocal microscope system using Volocity acquisition software (PerkinElmer). Z-stack confocal images were taken with an optical slick thickness of 0.3 mm, and projected Z-stacked images are shown. The scrambled probe provided with the RNAscope kit was used as a negative control.

Statistical analyses
With the exception of analyses related to the microarray and scRNA-seq data, all statistics were performed using GraphPad Prism 8 (GraphPad Software). Standard deviation values are given for the genes and transcripts per cell values for the scRNA-seq data. Scatter plots show standard error of the mean. Statistical significance was considered reached at p , 0.05. All statistical tests and conditions are described in the text. Graphs were generated in both GraphPad Prism 8 and Excel (Microsoft).

Data/software accessibility
Raw scRNA-seq and microarray datasets have been deposited in the GEO database under the ID codes GEO: GSE147285 and GSE146958. The Cellcellinteractnet code, the updated ligand-receptor database, and the raw proteomics data are available on request.

Defining growth factors made within the injured and uninjured peripheral nerve
We hypothesized that peripheral nerves promote axon growth in part due to ligand secretion by resident nerve Pdgfra-positive mesenchymal cells. To test this idea, we analyzed the sciatic nerve, which contains axons from sensory, sympathetic and motor neurons. Initially, we confirmed that Pdgfra-positive mesenchymal cells are found within the nerve endoneurium where axons are located. To do this, we analyzed control and injured sciatic nerve sections from mice carrying a transgene where EGFP is driven by Pdgfra regulatory sequences, and which thus tags Pdgfra-positive mesenchymal cells (Pdgfra Egfp/1 mice; Hamilton et al., 2003). We resected sciatic nerves from these mice, and at 9 DPI immunostained for the PDGFRa protein product and for the Schwann cell marker S100b . For comparison, we performed a similar analysis of the contralateral uninjured sciatic nerve. This Ligands identified in 3 DPI, 7 DPI, and uninjured rat sciatic nerves based on microarray analysis and our curated ligand-receptor database (modified from Yuzwa et al., 2016). Ligands were considered expressed if their expression levels exceeded that of Ntf3 mRNA. Ligands expressed only in the uninjured nerve are indicated by one asterisk, only in the 3 DPI nerve by two asterisks, and only in the 7 DPI nerve by one plus sign. Ligands expressed in the 3 and 7 DPI nerves but not in the uninjured nerve are indicated by two plus signs. Also shown in a separate column are ligands expressed in all populations (intersect). p uninjured only. pp 3 DPI only. 1 7 DPI only. 11 3 DPI and 7 DPI intersect.   (Continued) analysis identified many Pdgfra-EGFP-positive, PDGFRa protein-positive cells within the endoneurium of both the control and injured nerves (Fig. 1A), consistent with a similar recent analysis in Carr et al. (2019).
Having confirmed the presence of many mesenchymal cells in the endoneurium, we obtained an overview of ligands expressed in the sciatic nerve by performing global transcriptomic analysis of rat sciatic nerves that were either uninjured or that were injured by resection 3 or 7 d earlier. In all cases, we isolated total RNA from four independent biological replicates of the distal nerve segment and analyzed samples on Affymetrix GeneChip Rat Gene 2.0 ST Arrays. To define ligand mRNAs, we used a curated database of secreted ligands and their cognate receptors , establishing an expression cutoff by considering only mRNAs expressed at similar or higher levels than neurotrophin 3 (Ntf3) mRNA, which is expressed in and important for peripheral nerves. Ntf3 mRNA was expressed in the top 67.6%, 70.3%, and 71.2% of mRNAs in control, 3-d injured, and 7-d injured distal nerves, respectively. This analysis identified 238, 249, and 258 ligand mRNAs in the uninjured, 3-d injured, and 7-d injured nerves, respectively ( Fig. 1B; Table 1). Most (226) were detected in all three conditions, but 31 ligand mRNAs were found only in the injured nerves, including Areg, Ccl20, Il6, and Wnt7a and five ligands were detected only in the control nerve, including Il21 (Table 1).

Single-cell profiling defines similarities between peripheral nerves at 3 and 9 DPI
The microarray analysis defined the global nerve ligand environment. To ask which nerve cells express these ligands, we performed single-cell RNA sequencing (scRNA-seq), initially analyzing the distal sciatic nerve 3 d following a nerve transection. To do this, we dissociated the distal transected sciatic nerve, removed myelin debris with myelin-removal beads, and sequenced cells using high-throughput, droplet-based scRNA-seq (for details of all sequencing runs, see Materials and Methods). To analyze these individual transcriptomes, we then used a pipeline that incorporates extensive low-level data quality analysis with visualization and clustering methods that use evidence-based parameter selection (Innes and Bader, 2019), as previously described Carr et al., 2019;Storer et al., 2020). Genes with high variance were then used to compute principal components as inputs for projecting cells in two dimensions using t-SNE followed by graph-based clustering (Butler et al., 2018) with a range of resolution parameters.
This analysis defined nine distinct clusters containing 2075 cells in the 3-d injured sciatic nerve dataset (Fig.  1C). To assign cell types to these different clusters, we used well-characterized marker genes (for all cell typespecific marker genes used, see Materials and Methods; Fig. 1D; Extended Data Fig. 1-1A). This approach The combined 3 and 9 DPI scRNA-seq dataset (Fig. 1E) was analyzed for the percentage of cells within the defined nerve cell types that expressed ligand mRNAs as identified by microarray analysis (Table 1). Ligand mRNAs were only considered to be detected if they were expressed in at least 2% of cells within at least one cell type. BT = below threshold and indicates that ,2% of cells detectably expressed the ligand mRNA. Ligands annotated with one asterisk in the leftmost column were expressed in !2% Pdgfra-positive mesenchymal cells and two asterisks indicate ligands with the highest expression in either the epineurial/ perineurial or endoneurial Pdgfra-positive mesenchymal cells. p .2% expression in Pdgfra-positive cells. pp .2% expression and highest expression in Pdgfra-positive cells.   (Continued) identified Pecam1/Cd31-positive endothelial cells, Sox10-positive Schwann cells, Aif1/Iba1-positive macrophages, Trbc2-positive immune cells, and Rgs5-positive VSM and pericyte cells. As previously seen in the 9-d injured sciatic nerve (Carr et al., 2019), we also identified distinct populations of Pdgfra-positive nerve mesenchymal cells, including Dpp4-positive epineurial cells (clusters 1, 4, and 6) and Wif1-positive endoneurial cells (cluster 3).
We next asked whether distal sciatic nerve cells were similar at 3 and 9 d following a transection, taking advantage of a recently published 9 DPI scRNA-seq analysis where nerve cells were also isolated using myelin beads (Carr et al., 2019). To do this, we extracted the relevant raw transcriptomes from the previously published 9 DPI dataset, combined them with the 3 DPI raw transcriptomes we had generated, and put this combined 5395 cell dataset through our computational pipeline, which corrects for any differences in sequencing depth and library size. As seen with the 3 DPI dataset alone, marker gene analysis (  1F) showed that all cell type clusters were comprised of intermingled 3 and 9 DPI cells except for cluster 7, which included Cd19-positive B cells that came almost exclusively (99.4%) from the 9 DPI time point (Extended Data Fig. 1-1B). Similar results were obtained following batch correction. Thus, 3-and 9-d distal transected sciatic nerve cells are transcriptionally similar, although there is an increased proportion of B cells at 9 d.

Schwann cells, mesenchymal cells, vasculature, and immune cells make distinct contributions to the peripheral nerve ligand environment
To understand how the different cell types contributed to the injured nerve environment, we analyzed the combined 3-and 9-d scRNA-seq dataset for expression of ligands identified in the microarray analysis ( Fig. 1B; Table 1), excluding extracellular matrix proteins and ligands without well-defined, receptor-mediated paracrine roles. This analysis identified 143 ligands that were detectably expressed in ! 2% of at least one nerve cell type (Table 2), including many well-known nerve ligands such as the neurotrophins NGF, BDNF, and NT3, the Ret family ligands GDNF, Artemin and Neurturin, DHH, and various Semaphorin and FGF family members. Notably, Pdgfra-positive mesenchymal cells detectably expressed more ligand mRNAs than any other nerve cell type (118/143, or 82.5%; Fig. 1G; Table 2). Of these, 71 were expressed, proportionately, in more mesenchymal cells than any other cell type ( Fig. 1H; Table 2), and some were highly mesenchymally enriched, including Adm, Bmp7, Cxcl13, Fgf7, Fgf10, Gdf10, Hgf, Il33, Ntn1, Pthlh, and Wnt5a ( Fig. 1I; Extended Data Fig. 1-1D). A total of 39 of these mesenchymal ligands were most highly expressed in the endoneurial mesenchymal cells that are closely apposed to Schwann cells and axons ( Fig. 1G; Table 2).
Schwann cells were a second major source of injured nerve ligands, detectably expressing 74 ligand mRNAs (Fig. 1G,H; Table 2). A total of 23 of these were expressed in more Schwann cells than any other cell type, including Artn, Bdnf, Btc, Clcf1, Crlf1, Dhh, Fgf5, Gdnf, Hbegf, Sema3e, Sema4f, Shh, and Ucn2 ( Fig Table 2). Of the other cell types, endothelial and VSM/pericyte cells expressed 18 and 15 ligand mRNAs, respectively, at the highest levels ( Fig. 1G,H; Table 2). These included Bmp4 and Pdgfb mRNAs in endothelial cells and Bmp2 and Ngf in VSM/pericytes (Fig. 2B). The immune cells expressed 16 ligand mRNAs in the highest proportions including well-characterized immune ligands like Osm and Tnf (Figs. 1G,H, 2B; Table 2). Thus, ligands known to be important for axon growth and tissue regeneration are expressed by diverse nerve cell populations after injury.

Single-cell transcriptional profiling of uninjured and developing nerves
We asked whether this cellular profile of ligand expression was exclusive to the injured nerve by analyzing single-cell transcriptional datasets of the uninjured and developing sciatic nerves. For the uninjured nerve, we used a dataset that was previously analyzed for Pdgfrapositive mesenchymal cells (Carr et al., 2019), but not for other cell types. Analysis of the entire 1841 uninjured The uninjured nerve scRNA-seq dataset (Fig. 2C) was analyzed to determine the percentage of cells within a given cell type that detectably expressed (!2%) the 143 injured nerve ligand mRNAs (Table 2). Cell populations were defined as for the injured nerve analysis except that Schwann cells were divided into myelinating (cluster 8) and non-myelinating (cluster 4) cells. BT = below threshold and indicates that ,2% of cells detectably expressed the ligand mRNA. Ligands annotated with one asterisk in the leftmost column were expressed in !2% Pdgfra-positive mesenchymal cells and two asterisks indicate ligands with the highest expression in either the epineurial/perineurial or endoneurial Pdgfra-positive mesenchymal cells. p .2% expression in Pdgfra-positive cells. pp .2% expression and highest expression in Pdgfra-positive cells.   (Continued) nerve cell dataset (Fig. 2C,D; Extended Data Fig. 2-1B) identified clusters comprised of epineurial, perineurial, and endoneurial mesenchymal cells, as well as VSM/pericyte cells, endothelial cells, and a small population (1.9%) of immune cells. There were also two Sox10-positive Schwann cell populations; non-myelinating cells expressing Ngfr/p75NTR, Cdh2, and L1cam, and myelinating Schwann cells expressing high levels of Mbp, Pmp22, and Plp. This latter population was likely relatively reduced in numbers because myelin removal beads were used when isolating the nerve cells. There were very few proliferating cells, as indicated by expression of cell cycle genes like Top2a and Ki67 (Extended Data Fig. 2-1B).
For the developing neonatal nerve, we generated a new dataset, sequencing P2-P4 sciatic nerve cells after isolating them using either flow cytometry or myelin removal beads in two separate preparations. We combined raw transcriptomes from both preparations and analyzed them together. This analysis, which included 6885 total cells, identified 10 clusters containing intermingled cells from both preparations ( Fig. 2E; Extended Data Fig. 2-1C, D). This intermingling was unaffected by batch correction (Extended Data Fig. 2-1C). Cell type-specific marker genes identified endothelial cells, VSM/pericyte cells, macrophages, and epineurial, endoneurial, and perineurial mesenchymal cells (Fig. 2E,F; Extended Data Fig. 2-1D). At this age, there were five Sox10-positive Schwann cell clusters, with one cluster (6) containing proliferative Schwann cells and another (9) myelinating Schwann cells expressing high levels of Mag, Mbp, Mpz, and Pmp22 (Fig. 2E,F).
Analysis of ligand mRNA expression in these datasets showed that cells in the uninjured and developing sciatic nerves expressed many but not all injured nerve ligands. Specifically, of 143 total ligands, 122 and 119 were detectably expressed in uninjured and neonatal nerves, respectively, and 111 were shared in all three conditions ( Fig. 2G; Tables 3, 4). These ligand mRNAs were contributed by all of the different cell types in the uninjured and neonatal nerves (Fig. 2H,I). However, the pattern of ligand expression differed from the injured nerve (compare Figs. 2I, 1H), with endoneurial mesenchymal cells and Schwann cells contributing relatively fewer ligands in the developing and uninjured nerves (endoneurial cells, 11% and 27% in the uninjured vs injured nerves; Schwann cells, 5% and 16% in the uninjured vs injured nerves). Moreover, 13 injured nerve ligands were not detectably expressed in the neonatal or uninjured nerves (Artn, Btc, Ccl25, Crlf1, Fgf5, Gdnf, Gnrh1, Hgf, Rspo3, Sema4f, Shh, Tnfsf8, and Ucn2; Tables 3, 4). Notably, all of these were injured nerve Schwann cell or mesenchymal cell ligands (Table 2). Thus, multiple cell types contribute to the sciatic nerve ligand environment in all conditions, but mesenchymal and Schwann cells become relatively more important following injury.

Injured Schwann cells acquire a unique transcriptional phenotype following injury including upregulation of many growth factor genes
To ask about the apparent injury-associated increase in ligand expression, we analyzed the Schwann cells and Pdgfra-positive mesenchymal cells in more detail. We first combined transcriptomes from all Schwann cell clusters in the neonatal, uninjured, 3 DPI, and 9 DPI nerve datasets [cluster 6 ( Fig. 1E), clusters 4 and 8 (Fig. 2C), clusters 1, 2, 4, 6, and 9 ( Fig. 2E)]. We augmented this combined dataset by including Schwann cells from the FAC-sorted 9 DPI dataset from Carr et al. (2019). Once this combined dataset was put through the pipeline, we used the Harmony batch correction method (Korsunsky et al., 2019) to correct for any technical variation. Analysis of this combined dataset indicated that injured nerve Schwann cells were distinct from both developing and adult uninjured Schwann cells. Specifically, the combined dataset included 5331 Schwann cells in seven clusters (Fig. 3A,B). The differentiating neonatal Schwann cells were present in clusters 0, 1, and 2 with proliferating cells in cluster 2 (Fig. 3A,B; Extended Data Fig. 3-1A). Adult and neonatal myelinating Schwann cells were present in clusters 5 and 6, while the adult uninjured non-myelinating Schwann cells were in cluster 4. By contrast, almost all injured nerve Schwann cells were present in cluster 3.
To better understand these clusters, we performed hierarchical and correlation analyses of average gene expression (Fig. 3C). These analyses confirmed that the injured Schwann cells were distinct from the other populations, and indicated that they were more similar to the differentiating neonatal cells (r = 0.88 for the comparison between clusters 3 and 0) than to the adult non-myelinating Schwann cells (r = 0.76 for the comparison between clusters 3 and 4). To understand these similarities and differences at an individual cell level, we performed single-cell correlation analysis. As comparators for the analysis, we determined average gene expression for uninjured nonmyelinating versus neonatal non-myelinating Schwann The neonatal nerve scRNA-seq dataset (Fig. 2E) was analyzed to determine the percentage of cells within a given cell type that detectably expressed (!2%) the 143 injured nerve ligand mRNAs (Table 2). Cell populations were defined as for the injured nerve analysis. BT = below threshold and indicates that ,2% of cells detectably expressed the ligand mRNA. Ligands annotated with one asterisk in the leftmost column were expressed in !2% Pdgfra-positive mesenchymal cells and two asterisks indicate ligands with the highest expression in either the epineurial/perineurial or endoneurial Pdgfra-positive mesenchymal cells. p .2% expression in Pdgfra-positive cells. pp .2% expression and highest expression in Pdgfra-positive cells.  Continued) cells and for uninjured non-myelinating versus 9 DPI Schwann cells. We then compared each single-cell transcriptome with the averaged bulk transcriptomes using Pearson's correlation and used the resultant correlation coefficients to assign a two-dimensional coordinate to each cell. This analysis (Fig. 3D) showed that (1) virtually all 3 and 9 DPI Schwann cells were more similar to the neonatal Schwann cells than to the uninjured non-myelinating cells, (2) most neonatal Schwann cells were more similar to the injured cells than to the uninjured non-myelinating cells, and (3) despite these similarities, there was very little direct overlap between the injured and neonatal cells.
These data indicate that following nerve injury Schwann cells become more like neonatal Schwann cells, but that they are nonetheless distinct. In this regard, it has been reported that this unique injury state might involve acquisition of mesenchymal-like characteristics (Arthur-Farraj et al., 2017;Clements et al., 2017). To explore this idea further, we compared injured nerve Schwann cells and Pdgfra-positive mesenchymal cells from the combined 3 and 9 DPI nerve dataset (Fig. 1E). Correlation analysis showed that the injured Schwann cells were very distinct from both endoneurial and epineurial cells in the injured nerve (Extended Data Fig. 3-1B; r = 0.74-0.78). Thus, after injury, Schwann cells acquire a unique transcriptional profile that is similar but not identical to neonatal Schwann cells.

Upregulation of ligands in endoneurial mesenchymal cells following nerve injury
We performed a similar analysis of nerve mesenchymal cells, combining transcriptomes from the neonatal, The combined Schwann cell scRNA-seq dataset (Fig. 3A,B) was analyzed to determine the percentage of cells within the different Schwann cell clusters that detectably expressed (!2%) the 143 injured nerve ligand mRNAs (Table 2). Also shown is the difference, expressed as fold change, in the percentage of positive cells in the injured versus uninjured, non-myelinating Schwann cell clusters. BT = below threshold and indicates that ,2% of cells detectably expressed the ligand mRNA. Also shown are the absolute values, since these were used to determine the fold changes. Neonatal includes cells in clusters 0, 1, 2, and 5, uninjured myelinating includes cluster 6 cells, uninjured non-myelinating includes cluster 4 cells, and injured includes cluster 3 cells. Asterisks indicate ligand mRNAs with the highest expression in the injured Schwann cell cluster. pInjured, highest expression.   Continued) uninjured and injured nerve Pdgfra-positive clusters [clusters 1, 3, 5, and 10 ( Fig. 1E), clusters 1, 2, 6, and 10 ( Fig. 2C), clusters 3 and 5 ( Fig. 2E)], as well as the Pdgfra-positive mesenchymal transcriptomes of the FAC-sorted 9 DPI cells from Carr et al. (2019). Once this combined dataset was put through the pipeline, we used Harmony batch correction (Korsunsky et al., 2019) to correct for any technical variation. Analysis of this combined dataset showed that, as published previously (Carr et al., 2019), some mesenchymal populations were transcriptionally altered by injury, while others were largely unaffected. Specifically, the combined dataset included 5416 cells in nine Pdgfra-positive clusters (Fig. 4A,B). The injured and uninjured epineurial cells were coclustered, as were the injured, uninjured, and neonatal perineurial cells (Fig. 4A,B; Extended Data Fig. 4-1). By contrast, the uninjured, injured, and neonatal endoneurial cells were all segregated from each other. The other segregated clusters included neonatal epineurial cells (cluster 1) and the injured nerve differentiating bridge cells (cluster 0).
We used this combined dataset to ask about injury-induced ligand induction in mesenchymal cells. This analysis indicated that the endoneurial mesenchymal cells were largely responsible for this induction. Specifically, 102 of the 143 injured nerve ligands were detectably expressed in endoneurial mesenchymal cells (Fig. 4C), and, of these, 49 were expressed in at least 3-fold more injured versus uninjured cells, with 26 detectably expressed only in the injured cells (Table 6). These upregulated ligand mRNAs included Angpt1, Ccl9, Crlf1, Cxcl2, Inhbb, Lif, Sema7a, and Ngf (Fig. 4D,E; Table 6). In addition to this endoneurial cell response, some ligands were highest in the injured bridge cells, such as Bdnf, Cxcl9, and Hgf ( Fig. 4D; Table 6). By contrast, many ligand mRNAs were expressed to a greater or lesser degree in all nerve mesenchymal cell populations regardless of nerve injury, including for example Adm, Bmp1, Ccl11, Cxcl12, Il33, Pthlh, Fgf18, Pdgfa, Tgfb3, Vegfa, and Wnt5a (Fig. 4F,G). Thus, injury induces expression of many ligand mRNAs in endoneurial mesenchymal cells, but many ligands are also expressed under homeostatic conditions in uninjured nerve mesenchymal cells.

Identification of growth factor receptors on peripheral sympathetic and sensory neurons
To determine which of these ligands are likely to be important for axonal growth, we characterized growth factor receptors on sensory and sympathetic neurons which project their axons via the sciatic nerve. To do this, we coupled cell-surface proteomics and transcriptome profiling on purified neuronal cultures. For sensory neurons, we cultured E15 rat DRG neurons for 9 d in medium containing NGF. Immunostaining showed that these cultures were comprised of relatively pure b III-tubulin-positive neurons with 2-3% contaminating S100b -positive Schwann cells (Fig. 5A; Extended Data Fig. 5-1A). For sympathetic neurons, we isolated cells from the neonatal rat SCG and cultured them for 6 d in NGF. These cultures contained b III-tubulin-positive neurons and low percentages of Fibronectin-positive fibroblasts and S100b -positive Schwann cells ( Fig. 5A; Extended Data Fig. 5-1A).
Initially, we characterized the neuronal cell-surface proteomes, taking advantage of the fact that many cell-surface proteins are glycosylated. Specifically, we performed periodate oxidation of cell-surface glycans, bound modified proteins on a hydrazide resin after cell lysis, digested the bound proteins with trypsin and PNGase F, and identified peptides by mass spectrometry. For each sample, we analyzed three independent biological replicates. This The combined mesenchymal cell dataset (Fig. 4A,B) was analyzed to determine the percentage of cells within the different mesenchymal cell clusters that detectably expressed (!2%) the 143 injured nerve ligand mRNAs (Table 2). Also shown is the difference, expressed as fold change, in the percentage of positive cells in the injured versus uninjured endoneurial mesenchymal cells. BT = below threshold and indicates that ,2% of cells detectably expressed the ligand mRNA. Also shown are the absolute values, since these were used to determine the fold changes. Neonatal epineurial includes cluster 1 cells, neonatal endoneurial includes cluster 6 cells, inj/uninjured epineurial includes clusters 2 and 5 cells, neonatal/uninjured/injured perineurial includes cluster 8 cells, uninjured endoneurial includes cluster 4 cells, injured endoneurial includes cluster 3 cells, injured differentiating includes cluster 0 cells, and neonatal/injured proliferating includes cluster 7 cells. AceIII Acp2 Ache Ache Adam22 Acp2 Acp2 Cpd Enpp4 (Continued)  Ctsl L1cam    Tage4  Itfg1  Tenm3  Itga1  Tenm4  Itga3  Tfrc  Itga4  Thbs1  Itga5  Thsd7a  Itga6  Thy1  Itga7  Timp1  Itga9  Tll2  Itgal  Tm9sf3  Itgam  Tmed4  Itgav  Tmed7  Itgb1  Tmed9  Itgb8  Tmeff1  Jag1  Tmem106b  Jkamp  Tmem132a  Kcnc4  Tmem63b  Kdelc2  Tmem63c  Kiaa0319  Tmem87a  Kirrel3  Tmem87b  L1cam  Tpp1  Lama4  Trpv2  Lama5  TSLC1  Lamb1  Tspan3  Lamc1  Tspan6  Lamp1  Tspan8  Lamp2 Ttyh3      Tage4  Tapbp  Tctn1  Tctn2  Tenm2  Tenm3  Tenm4  Tfrc  Tgfb2  Tgfbr3  Thbs1  Thsd7a  Thsd7b  Thy1  Timp1  Tm2d1  Tm9sf3  Tmed4  Tmed7  Tmed9  Tmeff1  Tmem106b  Tmem132a  Tmem132c  Tmem132e  Tmem158  Tmem181  Tmem2  Tmem200c  Tmem231  Tmem255a  Tmem63b  Tmem63c  Tmem87a  Tmem87b  Tmem9b  Tmtc4  Tmx3  Tor1aip2  Tor2a  Tpbg  Tpcn1  Tpp1  Trhde  Trpv2  Tspan13  Tspan3  Tspan6  Tspan7  Tspan8  Ttyh3  Tuba1a  Tuba1b  Tuba1c  Tuba3a  Tubb2a Tubb2b (Continued) analysis identified 608 and 271 unique proteins on sensory and sympathetic neurons, respectively, with 219 of these common to both populations (Extended Data Fig.  5-1B; Table 7). The lower number of unique proteins on sympathetic neurons is likely due to decreased protein that was isolated (samples averaged %1100 vs 300 mg/ml for sensory vs sympathetic neurons). PANTHER classification identified most of these proteins as receptors, transporters and hydrolases, indicating appropriate enrichment for cell-surface proteins (Extended Data Fig.  5-1C). We then used the ligand-receptor database and manual curation to identify 102 proteins as receptors of various types, including G-protein-coupled receptors, receptor tyrosine kinases and phosphatases, cytokine receptors, and ligand-gated ion channels ( Fig. 5B; Table 8). Among these were well-characterized receptors found on both types of neurons such as TrkA (encoded by Ntrk1), BMP receptor 2, RET, gp130 (encoded by Il6st), and IGF2 receptor, receptors identified only on sensory neurons such as GFRa3 and receptors identified only on sympathetic neurons such as ALK and GFRa2. Cell-surface proteomics is relatively insensitive, and it ha5s previously been shown that more sensitive transcriptomic profiling can also be used to identify biologically relevant paracrine interactions (Johnston et al., 2016;Yuzwa et al., 2016;Voronova et al., 2017). We therefore complemented the proteomics by analyzing six and four independent biological replicates of cultured DRG and SCG RNA, respectively, on Affymetrix GeneChip Rat Gene 2.0 ST Arrays. To analyze these data, we defined an expression cutoff based on the proteomics data. Specifically, we identified the cell-surface receptor Gene symbols of proteins identified using cell-surface capture mass spectrometry on sensory neurons (column 1), sympathetic neurons (column 2), and both neuron types (column 3; intersect). Total numbers of proteins are indicated. Proteins included in this list were annotated by the terms "cell membrane" and/or "secreted" by the UniProtKB database (http://uniprot.org) and were verified by manual curation.  (47) Adra1b (Continued)   Epor Fzd9 Flt1 proteins with the lowest mRNA expression on the microarray for each neuron type and used those as the cutoffs. For sensory and sympathetic neurons, these were Itgam and Sorcs3 mRNAs, respectively (expressed at 87% and 81% of total mRNAs; Extended Data Fig. 5-1D). When these thresholds were applied to the microarray data, there were 321 and 297 receptor mRNAs in sensory and sympathetic neurons, respectively (Table 8). Importantly, there was good correspondence between the proteomics and microarray data; TrkA/Ntrk1, Bmpr2, Ret, and Igf2r mRNAs were similarly expressed in both populations of neurons; Rtn4r, Gfra3, and Acvr2a mRNAs were enriched in sensory neurons (2.8-, 10-, and 3.1-fold enriched, respectively); and Alk mRNA was 9-fold enriched in sympathetic neurons (p , 0.05 FDR for differences).

Computational modeling predicts that ligands deriving from multiple types of nerve cells, including mesenchymal cells, act on peripheral neurons
We performed computational modeling with the 143 injured nerve ligands and the sensory and sympathetic neuron receptors we had defined to predict how the injured nerve environment might regulate peripheral axon biology. This modeling predicted 122 and 125 potential unidirectional paracrine interactions between the injured nerve and sympathetic and sensory neurons, respectively (Fig. 5C,D; Table 9). Of these, cell-surface receptor protein expression was detected for 49 and 60 sympathetic and sensory neuron predicted interactions (Fig. 5C,D, blue boxes). Many predicted interactions involved known peripheral nerve ligands such as the neurotrophin and GDNF families. Notably, all but three ligands (GNRH1, CXCL1, and CXCL2) were predicted to act on both sympathetic and sensory neurons. The receptors for these predicted interactions were also largely the same, except for several sensory neuron receptors; the Erbb4 receptor for EGF/neuregulin family ligands, KDR for the VEGF family, ACVR1C for the activin/BMP family, and the CXCR5 chemokine receptor (Table 9). Thus, the injured nerve is predicted to produce ligands that act on both  Tgfbr2  Ror1  Tgfbr3  Ror2  Thbd  Rorb  Thra  Rtn4r  Thrap3  Rtn4rl1  Tnfrsf10b  Rxfp1  Tnfrsf11a  Rxrg  Tnfrsf11b  Ryr1  Tnfrsf12a  Ryr2  Tnfrsf13c  Sctr  Tnfrsf14  Sdc4  Tnfrsf17  Sfrp1  Tnfrsf18  Sfrp2  Tnfrsf1a  Slc1a5  Tnfrsf1b  Sorcs3  Tnfrsf25  Sort1  Tnfrsf4  Sstr1  Tnfrsf8  Sstr2  Tnfrsf9  Sstr3  Tshr  Sstr4  Unc5b  Sstr5  Unc5c  Tek  Uts2r  Tgfbr2  Vipr1  Tgfbr3  Vldlr  Thbd  Vtn  Thra  Xcr1  Thrap3  Tnfrsf10b  Tnfrsf11a  Tnfrsf11b  Tnfrsf12a  Tnfrsf13c  Tnfrsf14  Tnfrsf17  Tnfrsf18  Tnfrsf1a  Tnfrsf1b  Tnfrsf25  Tnfrsf4  Tnfrsf8  Tnfrsf9  Tshr  Unc5b Unc5c (Continued)  Gene symbols of proteins identified by mass spectrometry exclusively on sensory neurons (DRGs), exclusively on sympathetic neurons (SCGs), or on both neuron types, which were identified from the protein lists shown in Table 7. Also shown are receptor mRNAs identified by microarrays as expressed by DRG and SCG neurons, defined using the updated ligandreceptor database (modified from Yuzwa et al., 2016). Only receptor mRNAs that had expression exceeding the cutoffs for each neuron type (DRGs; Itgam, 87% and SCGs; Sorcs3, 81%) are included. The total numbers of receptors in each column are indicated.                (Continued) sympathetic and sensory neurons, largely through the same receptors.
We then used the scRNA-seq data to define the nerve cell types that made the ligands involved in the predicted interactions (Figs. 5C,D, 6A,B; Extended Data Fig. 6-1A, B). Almost half (57) of 122 shared predicted sympathetic and sensory neuron interactions involved ligands that were expressed in the highest proportions in Pdgfra-positive mesenchymal cells including FGF10, FGF18, HGF, SEMA3D, BMP7, IL33, and PTHLH (asterisks in the models denote ligands made by 4-fold more of the indicated cell relative to all other cell types). Moreover, 33 of these were highest in the endoneurial mesenchymal cells. Another 22 ligands were highest in Schwann cells, including ARTN, BTC, DHH, FGF5, GDNF, SEMA3B, SHH, and UCN2. The remaining 43 predicted interactions were split almost equally between endothelial cells (17), VSM/pericyte cells (14), and immune cells (12), and included wellcharacterized nerve ligands such as NGF, NT3 (Ntf3), and IGF2. Notably, while some of these ligands were predicted to signal via dedicated receptors (Fig. 6A,B; Extended Data Fig. 6-1A,B; Table 9), many others shared receptors or coreceptors, suggesting the potential for convergent signaling. As examples, CLCF1, CRLF1, and LIF were all predicted to share the receptor components LIFR and gp130, and ANGPT1, ANGPT2, and ANGPT4 were all predicted to act by interacting with TEK.

Modeling predicts that many nerve ligands have the capacity to act on both PNS and CNS neurons
This analysis predicts that peripheral nerve cells produce ligands that could act on at least two populations of peripheral neurons. To test the idea that this might reflect a ligand environment that is generally supportive of axonal growth, we asked about motor neurons, which also project axons via the sciatic nerve, and RGCs, CNS neurons that normally do not regenerate in the CNS, but will regenerate into peripheral nerve grafts (Politis and Spencer, 1986; for review, see Benowitz et al., 2017).
For motor neurons, we analyzed previously published microarray data from microdissected P7 mouse lumbar motor neurons (Kaplan et al., 2014), identifying receptor mRNAs using the same thresholding cutoff as for the sensory neurons (top 87% of mRNAs). Of 322 receptor mRNAs defined in the motor neuron dataset using this approach, 272 were also expressed by sensory and sympathetic neurons (Table 10). Computational modeling with the motor neuron receptors and the 143 injured nerve ligands showed that of the 122 shared sympathetic and sensory neuron interactions, 121 were also predicted for motor neurons, with the endoneurial ligand EDA the only exception (Fig. 7A,C; Table 9). There were also two predicted nerve to motor neuron interactions involving GRP-GRPR and GNRH1-GNRHR that were not shared with both sympathetic and sensory neurons.
We obtained similar findings when we analyzed RGCs using a bulk RNA-seq dataset generated from P5-P7 rat RGCs that were cultured for 12 h (Blanco-Suarez et al., 2018). Using an FPKM of 1 as a threshold for expression, we found 320 receptors expressed by RGCs, with 258 of them also expressed by the three peripheral neuron types (Table 10). Modeling of the unidirectional paracrine interactions between the injured nerve and RGCs identified 126 predicted interactions that included all 122 shared sensory and sympathetic neuron interactions (Fig. 7B, shared interactions are outlined in purple; Table 9). There were also four interactions that were unique to the RGC model involving CCL9, CXCL1, CXCL2, and CXCL16. Thus, the injured peripheral nerve and in particular the Schwann cells and endoneurial mesenchymal cells are predicted to provide a ligand environment that acts on multiple populations of neurons.

The communication networks identify mesenchymalderived ligands that regulate peripheral axon growth
These models predict that, in addition to Schwann cells, nerve mesenchymal cells are key growth factor sources in the injured nerve. To validate this concept, we focused on two predicted endoneurial ligands that have not previously been explored within a nerve context, ANGPT1 and CCL11 (Fig. 8A,B). We also analyzed VEGFC, which is expressed at a low level in injured nerve mesenchymal cells (Fig. 8A), but that had a coreceptor (NPR1) that was validated at the protein level in both sensory and sympathetic neurons.
We first asked whether these three ligands were secreted by nerve mesenchymal cells. To do this, we isolated mesenchymal cells from rat sciatic nerves using antibody-based Predicted unidirectional ligand-receptor interactions between the injured sciatic nerve ligands and receptors on sympathetic (SCGs) neurons, sensory (DRGs) neurons, motor neurons (MNs), and retinal ganglion cells (RGCs). Directionality of predicted paracrine interactions are indicated by the source cell (ligand) and target cell (receptor) column designations. Some ligands are predicted to act on more than one receptor. The interactions shown here accompany the models presented in Figures 5-7 and Extended Data Figure 6-1.
flow sorting for cell-surface PDGFRa. We cultured and expanded these sorted mesenchymal cells for 2.5-4 weeks, at which point the cultures were comprised of over 90% PDGFRa-positive mesenchymal cells, with the remaining cells being S100b -positive Schwann cells (Fig. 8C,D). After three further days in culture, we added defined, serum-free medium for 24-96 h, collected this conditioned medium, and performed ELISAs. ANGPT1, CCL11, and VEGFC were all detected in three independent preparations of nerve mesenchymal cell conditioned medium (Fig. 8E).
Having confirmed that these ligands were secreted by nerve mesenchymal cells in culture, we asked whether Angpt1, Ccl11, and Vegfc mRNAs were expressed in endoneurial cells of the injured nerve. To do this, we analyzed injured distal nerve sections from the Pdgfra Egfp/1 mice. We resected sciatic nerves and at 9 DPI performed single molecule FISH. Angpt1 and Ccl11 mRNAs were detectable in many Pdgfra-EGFP-positive mesenchymal cells within the injured nerve endoneurium (Fig. 8F). Consistent with the scRNA-seq data, many but not all EGFP-positive cells were positive for these mRNAs (Fig.  8A,B; 24% and 87% of endoneurial cells express detectable Angpt1 and Ccl11 mRNAs, respectively). There were also some Angpt1 or Ccl11-positive cells that were  Figure 5 (Schwann cell ligands in gray and endoneurial mesenchymal cell ligands in yellow). Receptors are shown on either side of the ligand column and also include coreceptors that are well-characterized components of receptor complexes. Receptors that were observed at both the transcriptomic and proteomic levels are colored green while those defined only at the transcriptomic level are colored blue. Arrows indicate directionality of interactions. Note that many ligands interact with multiple receptors and, conversely, that multiple ligands are sometimes predicted to share receptors.
Pdgfra-EGFP-negative. For Angpt1 mRNA these are likely VSM/pericyte cells, while for Ccl11 mRNA, they could be VSM/pericytes or immune cells (Fig. 8B). Vegfc mRNA was also detected in scattered endoneurial Pdgfra-EGFPpositive cells but consistent with the scRNA-seq analysis (Fig. 8A), there were fewer double-labeled cells and the FISH signal was low (Fig. 8F).
Having validated their expression, we asked whether ANGPT1, CCL11, or VEGFC could promote growth of peripheral axons. To do this, we used compartmented cultures of neonatal SCG sympathetic neurons (Campenot et al., 1991;Park et al., 2010). In these cultures, cell bodies are physically separated from axons so that ligands can be applied just to the axons, as would occur in the peripheral nerve (Fig. 8G). We established sympathetic neurons in these compartmented cultures with 10 ng/ml NGF, their obligate survival factor, in all compartments. Three days later, when axons had crossed into the side compartments, we replaced the side compartment medium with medium containing 0.5 ng/ml NGF plus 100 ng/ml of ANGPT1, CCL11, or VEGFC for three additional days. As a positive control, we added 50 ng/ml NGF into side compartments, and as a baseline control, we maintained axons in 0.5 ng/ml NGF alone. To quantify the density of axonal growth in these compartments, we drew a line perpendicular to the axis of the outgrowth within the furthest 1 mm of outgrowth where axons were maximally Tek Sstr4  Tgfbr1  Sstr5  Tgfbr2  Tek  Tgfbr3  Tgfbr1  Thbd  Tgfbr2  Thra  Tgfbr3  Thrap3  Thbd  Tnfrsf10b  Thra  Tnfrsf11a  Thrap3  Tnfrsf11b  Tnfrsf10b Tnfrsf12a (Continued) Tnfrsf11a  Tnfrsf13b  Tnfrsf11b  Tnfrsf13c  Tnfrsf12a  Tnfrsf14  Tnfrsf13b  Tnfrsf17  Tnfrsf13c  Tnfrsf18  Tnfrsf14  Tnfrsf1a  Tnfrsf17  Tnfrsf1b  Tnfrsf18  Tnfrsf25  Tnfrsf1a  Tnfrsf4  Tnfrsf1b  Tnfrsf8  Tnfrsf25  Tnfrsf9  Tnfrsf4  Trhr  Tnfrsf8  Tshr  Tnfrsf9  Unc5b  Trhr  Unc5c  Tshr  Uts2r  Unc5c  Vipr1  Uts2r Vipr2 Vipr1 Vldlr Vipr2 Vtn Vldlr Xcr1 Vtn Xcr1 Receptor mRNAs identified by microarray for motor neurons (MNs; column 1) and by RNA-seq for RGCs (column 3), as defined using the updated ligand-receptor database (modified from Yuzwa et al., 2016). Only included are receptor mRNAs that had expression levels exceeding a cutoff of the top 87% of mRNAs for motor neurons and FPKM of 1 for RGCs. Also shown are receptor mRNAs that overlap between sensory, sympathetic and motor neurons (column 2) or between all four populations of neurons (column 4). The sensory and sympathetic neuron receptor mRNAs are shown in Table 8. The total numbers of receptor mRNAs in each column are indicated. defasciculated. This analysis showed that 50 ng/ml NGF caused a robust increase in the density of sympathetic axons relative to the 0.5 ng/ml NGF baseline control (Fig.  8H,I). Notably, all three mesenchymal ligands modestly but significantly enhanced axonal density, although to a lesser degree than maximal NGF (Fig. 8H,I). Thus, at least some of the mesenchymally derived ligands predicted in our models were bioactive on sympathetic axons.

Discussion
In the present study, we have characterized the ligand environment of the uninjured, injured, and developing sciatic nerves using bulk and single-cell transcriptional profiling. We have identified receptor proteins on the surface of two types of peripheral neurons (sensory and sympathetic) and made predictions of ligand-receptor paracrine interactions between the injured nerve and peripheral neurons. We then go on to show, based on these predictions, that mesenchymal cells express factors that are capable of augmenting growth of peripheral axons in vitro, indicating that at least some of these ligands may directly contribute to the positive axonal growth environment of the developing and regenerating peripheral nerves.
Peripheral nerves provide a highly supportive environment for axonal growth during development and following injury (Chen et al., 2007;Cattin and Lloyd, 2016;Fledrich et al., 2019) and promote the repair and regeneration of innervated tissues (Kumar and Brockes, 2012;Johnston et al., 2013Johnston et al., , 2016Mahmoud et al., 2015). Many nerve-derived growth factors have already been well studied, including NGF, BDNF, NTF3, GDNF, and cytokines of the LIF/CNTF family (for review, see Terenghi, 1999). These factors are generally assumed to be Schwann cell derived, although macrophages express factors like GAS6 that promote proper function of Schwann cells in the regenerating nerve (Stratton et al., 2018). To identify other important factors that might be involved in providing a supportive peripheral nerve environment, we used a modeling strategy based on transcriptomic analysis and cellsurface mass spectrometry, as has been previously done for embryonic cortical development and digit tip regeneration (Johnston et al., 2016;Yuzwa et al., 2016). However, while bulk transcriptomic profiling was used in these earlier studies, here we added single-cell transcriptional profiling, thereby providing a level of cellular resolution previously not possible for complex tissues. This approach allowed us to define a previously unappreciated role for mesenchymal cells in establishing the nerve paracrine environment, to identify new nerve ligands, and to predict that many nerve ligands will act on both PNS and CNS neurons, thereby potentially providing an explanation for why peripheral nerves can promote growth of CNS axons.
Our study defined many growth factor mRNAs induced in Schwann cells following nerve injury. Some of these encoded previously studied factors like Artn, Bdnf, Gdnf, Pdgfa, Shh, and Lif, while others encoded factors not well studied in this regard, including Ucn2, Fgf5, and the CNTF-like cytokines Clcf1 and Crlf1. Previous studies have proposed that this ligand induction is part of a unique Schwann cell "repair" phenotype (Jessen and Mirsky, 2019) that is important for axonal regeneration in the case of ligands like BDNF, GDNF, and LIF, and for tissue repair, in the case of PDGFA (Johnston et al., 2016). What is this repair phenotype? Previous work has shown that following injury Schwann cells display altered morphology and gene expression that is thought to be conducive to promoting axonal regeneration Jessen and Mirsky, 2019). Repair Schwann cells have also been reported to have enhanced epithelial to mesenchymal transition gene expression and TGFb signaling that contributes to the establishment of an invasive, "mesenchymal-like" phenotype (Arthur-Farraj et al., 2017;Clements et al., 2017). Our findings also show that following nerve injury, Schwann cells induce ligand genes that are not expressed at detectable levels in the uninjured neonatal or adult nerves. Moreover, our transcriptional comparisons expand on previous work and show that repair Schwann cells are more similar to neonatal than to adult uninjured nerve Schwann cells but that they are nonetheless distinct. In this regard, our developmental comparison was limited to the neonatal nerve when myelination is ongoing, and it would be interesting to use single-cell transcriptional profiling to see how similar repair Schwann cells are to embryonic nerve Schwann cells before myelination has commenced.
An important finding of this work is that Pdgfra-positive mesenchymal cells are a major source of ligands in the developing, adult, and injured nerves and that they express well-characterized nerve growth factors like NGF, HGF, and LIF. Of particular interest is the high ligand expression by endoneurial mesenchymal cells, which are neural crest derived (Joseph et al., 2004) and are scattered throughout the endoneurial space in close apposition to Schwann cells and axons. These endoneurial mesenchymal cells are thus ideally positioned to regulate axon and Schwann cell biology, and, like Schwann cells, they display increased expression of many ligand mRNAs following injury, including well-studied ligands like Crlf1, Ngf, and Lif and less-studied ligands such as Angpt1, Ccl9, and Sema7a. Equally intriguing was the observation that endoneurial cells express many little-studied ligands under homeostatic conditions, including Adm, Bmp7, Il33, Pthlh, and Wnt5a. Since mesenchymally derived continued by computational modeling using the ligand-receptor database, and then manually curated for well-validated interactions. Nodes represent ligands that are color coded to identify the injured nerve cell type with the highest expression of the ligand mRNA based on the scRNA-seq analysis. Arrows indicate directionality of interactions. Predicted interactions that were shared among all four injured nerve-neuron models are indicated by a purple box around the corresponding ligand nodes. C, Venn diagrams showing the overlap of predicted ligand-receptor interactions between injured nerve mesenchymal cells (left) or Schwann cells (right) and sympathetic neurons (SCGs), sensory neurons (DRGs), and motor neurons (MNs). Figure 8. Identification and characterization of ligands expressed in nerve mesenchymal cells that locally promote sympathetic axon growth. A, tSNE gene expression overlays of Angpt1 and Vegfc on the combined 3 and 9 DPI nerve dataset shown in Figure  1E. Cells that detectably express the ligand are colored blue, and the numbers correspond to the clusters. B, Heatmap showing expression of Ccl11 mRNA in single cells within clusters of the combined injured nerve scRNA-seq dataset shown in Figure 1E. Each column line represents the level of expression in a single cell. Gene expression represents scaled expression values using Seurat's scaling function and is color coded as per the adjacent color key, where yellow indicates the highest expression. Cluster numbers are on the bottom, and the cell types in that cluster are annotated on the top. Mes. = mesenchymal, Endoth. = endothelial, Endo ligands include both well-studied nerve growth factors such as NGF as well as many ligands with unknown roles in nerve biology, it is likely that mesenchymal cells express ligands that are critical components of the regenerative response of the injured nerve and/or the growth program of the developing nerve. Our validation studies indicate that at least some of these endoneurial cell ligands are active on axons, but they might be equally important for other nerve cell types and/or for the tissues they innervate. As one example, PTHLH and nerve innervation are both important for bone homeostasis and repair (Elefteriou et al., 2014;Ansari et al., 2018), and endoneurial cells, which express Pthlh, migrate into the injured bone where they directly contribute to bone repair (Carr et al., 2019). As a second example, the vasodilator peptide adrenomedullin (ADM) was previously shown to stimulate cAMP accumulation in endothelial cells and Schwann cells (Dumont et al., 2002), suggesting that endoneurial cell-derived ADM might be important for nerve vasculature and/or Schwann cell biology. As a final example, BMP7 inhibits myelin gene expression in Schwann cells (Liu et al., 2016) and promotes mammalian digit tip regeneration (Yu et al., 2010), suggesting that endoneurial cell-derived BMP7 might play multiple important roles.
Another important finding is that most injured nerve ligands are predicted to act on all three populations of PNS neurons as well as RGCs. In this regard, the nerve could promote axonal growth and regeneration in two somewhat disparate ways. In one model, Schwann cells and endoneurial mesenchymal cells would produce different ligands depending on the axons that they are currently or have previously interacted with, thereby tailoring the nerve environment to the axons that need to grow or regenerate. Support for this model comes from studies showing that denervated Schwann cells of motor versus sensory nerves provide ligands specific to different types of axons (Höke et al., 2006;Brushart et al., 2013). In a second model, during development or following nerve damage Schwann cells and mesenchymal cells could express a broad repertoire of ligands, thereby ensuring that growth of all types of PNS axons would be supported. Our findings support this latter model, since we find broad injury-induced ligand expression and a relatively broad repertoire of receptor expression on four different types of neurons, culminating in many similar predicted paracrine interactions. Such a mechanism would provide maximum flexibility and would explain why peripheral nerve grafts promote regeneration of multiple types of injured CNS neurons, which do not normally project in peripheral nerves. Nonetheless, our findings are still consistent with the finding of differential ligand expression in different nerve subtypes (Höke et al., 2006;Brushart et al., 2013), since we have only defined the ligand landscape in a mixed nerve.
How predictive are these models? Previous studies using bulk transcriptomics and/or cell-surface mass spectroscopy predicted three factors important for embryonic cortical neurogenesis (IFNg , Neurturin, and GDNF; Yuzwa et al., 2016) and one for oligodendrogenesis (Fractalkine; Voronova et al., 2017), and two factors important for digit tip regeneration (PDGFA and OSM; Johnston et al., 2016). In those studies, the cell of origin for each ligand was identified either by isolating cells or by performing single-cell resolution morphologic approaches. Here, we have instead used single-cell transcriptional profiling to provide the necessary resolution, an approach with much broader applicability. The validity of the resultant models is attested to by our finding that almost all ligands previously shown to be important for peripheral nerve regeneration were independently assigned in our models, including many ligands known to be expressed in Schwann cells. Nonetheless, to ensure the validity of these models, we also examined three ligands that were predicted to be made by nerve mesenchymal cells, ANGPT1, CCL11, and VEGFC. Two of these factors, ANGPT1 and VEGFC, are well-known angiogenesis factors, while the third, CCL11 or eotaxin, is a chemokine involved in eosinophil recruitment (Jose et al., 1994). None of the three has been studied as a positive factor within the context of the injured nerve, although ANGPT1 has been shown to promote growth of cultured sensory neurons (Kosacka et al., 2006), VEGFC promotes growth of developing motor neurons in zebrafish (Kwon et al., 2013), and CCL11 inhibits Schwann cell differentiation (Büttner et al., 2018). Based on our predictive models, we tested these ligands and found that all three (1) were expressed by endoneurial mesenchymal cells in the injured nerve, as shown by both scRNA-seq and FISH analyses; (2) were synthesized and secreted by cultured nerve-derived PDGFRa-positive mesenchymal cells; and (3) enhanced sympathetic axon outgrowth when applied locally in the presence of minimal NGF. While we recognize that additional studies will be required to show that these three ligands are secreted by mesenchymal cells in vivo, our data indicate that they are highly expressed following nerve injury, raising the possibility that they are important for nerve repair. In this regard, ANGPT1, CCL11, and VEGFC were not as potent as NGF in promoting growth of sympathetic axons in culture, but they did enhance growth and thus could be important factors for peripheral axon growth in a regenerating nerve context. We therefore feel that our studies validate the predictive value of the modeling approach and provide support for the idea that mesenchymal cells within the nerve are important ligand sources, particularly within the context of nerve repair and regeneration.
The data presented here reinforce the importance of Schwann cells as sources of growth promoting factors and provide evidence that mesenchymal cells also play an important role in determining the ligand environment of the developing and injured nerve. Notably, some wellknown nerve regeneration ligands such as GDNF and BDNF were expressed at the highest levels in Schwann cells, while others, such as NGF and HGF, were instead highest in mesenchymal cells. In addition, both cell types express ligands that are not well-characterized as nerve injury ligands, including BTC and UCN2 for Schwann cells, and ANGPT1, CCL11, and VEGFC for mesenchymal cells. While the relative contributions of growth factors from these two cell types to nerve growth and repair remain unknown, our study does highlight mesenchymal cells as a previously overlooked reservoir of growth factors for axon growth and potentially for nerve tissue regeneration, an area we are only now starting to understand (for examples, see Parrinello et al., 2010;Cattin et al., 2015). The data presented here thus provide an important step toward defining nerve paracrine interactions not only with regard to axon growth and peripheral nerve regeneration, but also with regard to the paracrine roles of the nerve during repair and regeneration of target tissues.