Generation of iPSC-Derived Human Peripheral Sensory Neurons Releasing Substance P Elicited by TRPV1 Agonists.

Neural crest stem cells (NCPCs) have been shown to differentiate into various cell types and tissues during embryonic development, including sensory neurons. The few studies addressing the generation of NCPCs and peripheral sensory neurons (PSNs) from human induced pluripotent stem cells (hiPSCs), generated sensory cells without displaying robust activity. Here, we describe an efficient strategy for hiPSCs differentiation into NCPCs and functional PSNs using chemically defined media and factors to achieve efficient differentiation, confirmed by the expression of specific markers. After 10 days hiPSCs differentiated into NCPCs, cells were then maintained in neural induction medium containing defined growth factors for PSNs differentiation, followed by 10 days in neonatal human epidermal keratinocytes- (HEKn-) conditioned medium (CM). We observed a further increase in PSN markers expression and neurites length after CM treatment. The resulting neurons elicited action potentials after current injection and released substance P (SP) in response to nociceptive agents such as anandamide and resiniferatoxin. Anandamide induced substance P release via activation of TRPV1 and not CB1. Transcriptomic analysis of the PSNs revealed the main dorsal root ganglia neuronal markers and a transcriptional profile compatible with C fiber-low threshold mechanoreceptors. TRPV1 was detected by immunofluorescence and RNA-Seq in multiple experiments. In conclusion, the developed strategy generated PSNs useful for drug screening that could be applied to patient-derived hiPSCs, consisting in a powerful tool to model human diseases in vitro.

Introduction

Human induced pluripotent stem cells (hiPSCs) are used to generate different neuronal types to study biology, pharmacology and to screen for potential new drugs. Developing human peripheral sensory neurons (PSNs) are of great interest, not only to advance pain research aiming at finding new analgesics, but also to the cosmetic industry to develop alternative in vitro test methods (Leist et al., 2012; McNeish et al., 2015; Yekkirala et al., 2017). In the field of skin irritation, the activity of transient receptor potential (TRP) channels are of valuable importance, since these channels detect a range of topical nociceptive irritants, such as capsaicin and mustard oil (Nilius and Szallasi, 2014). Capsaicin has been used as a reference in clinical dermatology for skin neurosensitivity and reactivity, with minimal inter individual variation within different ethnic backgrounds (Jourdain et al., 2009; Gueniche et al., 2014).

Transient receptor potential channels comprise a diverse family of ligand-gated, mostly non-selective, cation channels that are robustly expressed in sensory systems throughout species (Nilius and Szallasi, 2014). Of these, TRPV1 is the most well-studied and is considered to be the prototypical TRP channel present in somatosensory neurons (Basbaum et al., 2009). TRPV1 can be directly gated by external molecules such as capsaicin and resiniferatoxin, and also modulated positively or negatively via activation of other receptors and second messenger systems, such as PIP2 hydrolysis and PKC phosphorylation (Julius, 2013). One of the receptors that seem to inhibit TRPV1 activation is the cannabinoid 1 receptor (CB1), also present in somatosensory neurons (Julius and Basbaum, 2001). However, an endogenous agonist of CB1, anandamide, is also a TRPV1 agonist, albeit with an EC50 one order of magnitude higher in the latter (Zygmunt et al., 1999).

Some attempts were made to generate somatosensory neurons from hiPSCs, but so far, these cells have been proved to be difficult to generate robustly. For instance, Chambers and collaborators developed a differentiation protocol based on a small molecule screen and successfully obtained neurons that express the main expected markers (Chambers et al., 2012). Nevertheless, although most neurons exhibited those markers, few showed important functional readouts, such as TRPV1 activation, as only about 1–2% of the cells responded to capsaicin. Other studies were able to demonstrate the expression of canonical peripheral markers, such as BRN3A, peripherin and β-tubulin III, in sensory neurons derived from hiPSCs, but evidence of specific functional activity was still missing (Menendez et al., 2013; Young et al., 2014; Blanchard et al., 2015; Eberhardt et al., 2015). Young and collaborators reported measuring TRPV1 activity induced by capsaicin, but only after 6 weeks in medium containing growth factors (Young et al., 2014). Another group showed general electrophysiological activity (i.e., action potentials, sodium currents) of human sensory neurons in culture, but did not show responses to capsaicin (Wainger et al., 2015). In the same paper, neurons obtained from mice with a similar method showed robust TRPV1 activity.

Keratinocytes and sensory neurons have an extensive interplay during development and within mature skin. For instance, keratinocytes release neurotrophic factors that induce branching of free nerve endings and neurite outgrowth toward the skin surface (Albers and Davis, 2007). They also release inflammatory mediators involved in responses to tissue damage and hypersensitivity reactions, as well as responses to cold and heat, through receptors of the TRP family of cation channels (Chung et al., 2004). On the other hand, sensory endings do not only transduce sensory signals, but play an active role in the cutaneous metabolism and homeostasis, through the secretion of pro-inflammatory neuropeptides and inflammatory mediators that control vascularization and tissue renewal (Roosterman et al., 2006). Particularly, TRPV1-positive nociceptors may also regulate skin longevity and metabolism, as well as the immune response over aging, as shown in TRPV1 knock-out mice (Riera et al., 2014).

In the glabrous skin, the sensory information is processed by four types of mechanoceptors: Merkel cells, Pacinian corpuscles, Meissner and Ruffini terminations. In the hairy skin, tactile transmission is associated with the lanceolate Aβ-LTMR terminations. Hair follicles are innervated by the longitudinal lanceolate terminations C-LTMR, Aδ-LTMR, and Aβ-LTMR. The intercommunication between nerve endings and epidermal keratinocytes occurs through neuropeptides involved in nociceptive and pruriceptive sensitization (Abraira and Ginty, 2013).

Substance P (SP) is a neuropeptide, member of the tachykinin family, synthesized by sensory neurons that emit their extensions from the dorsal root ganglion (DRG) to the more superficial layers of the skin, mediating the communication between peripheral neurons and epidermal keratinocytes (Ribeiro-da-Silva and Hokfelt, 2000). Most of the neurons that release substance P are sensitive to capsaicin, highlighting the importance of TRPV1 expression and sensory neurons-keratinocytes interplay.

Literature analysis reveals the efforts made by several groups to develop more efficient and cost-effective protocols for obtaining functional peripheral neurons from hiPSCs (Menendez et al., 2013; Young et al., 2014; Blanchard et al., 2015; Eberhardt et al., 2015). However, none of these studies have shown results with human cells or investigated the role of the interaction between human epidermal keratinocytes and the transcriptional profile of mature neurons in the presence of factors released by these cells.

In the present work, we developed a new protocol for the direct differentiation of hiPSCs to PSNs in the presence of conditioned medium (CM) obtained from HEKn. These neurons showed some expected functional responses, such as neuropeptide release. The analysis of differential gene expression confirms the transcriptional profile of putative human somatosensory peripheral neurons.

Materials and Methods

Cell Culture

Human induced pluripotent stem cells (GM23279A from Coriell Institute for Medical Research) were cultured on matrigel-coated plates with mTeSR (Stemcell Technologies, Vancouver, BC, Canada). The cells were checked for pluripotency markers expression such as Nanog, Sox2, Tra-1-60, and Tra-1-81 (data not shown), and were negative for the neuronal markers Islet1, BRN3A, peripherin, and TRPV1 (Supplementary Figure ).

When approximately 60–70% confluence was reached, cells were split using a 0.5 mM EDTA solution in PBS (Life Technologies, Rockville, MD, United States) (d0, Figure ). When the same confluence was achieved again, which could take between 1 and 2 days, the differentiation was initiated by switching to 3N medium (DMEM-F12, Neurobasal Medium, Glutamax, N2 Supplement, B27 Supplement, non-essential amino acids and β-mercaptoethanol) containing 500 nM LDN-193189 (Stemgen, Cambridge, MA, United States) and 10 μM SB431542 (Sigma, St. Louis, MO, United States). In the following day (d2), the same two factors were added again, together with 3 μM CHIR-99021 (Tocris Bioscience, Minneapolis, MN, United States) (d0–d5, Figure ). On day 3 (d3), the medium was changed to have LDN removed and on day 4 (d4) SB was not added, maintaining just CHIR until d10, with medium changes at odd days. Treatment with these Smad pathway inhibitors leads to differentiation of hiPSCs cells toward neural crest progenitor cells (NCPCs) (Mica et al., 2013) (Figure ). On day 11 (d11) of differentiation, NCPCs were plated onto Polyornithine (PLO)/Laminin-coated plates and maintained in 3N medium for neuronal differentiation with addition of 10 ng/ml FGF-2 (Life Technologies, PHG0263, Rockville, MD, United States) and 10 ng/ml EGF (Life Technologies, PHG0313, Rockville, MD, United States) for approximately 2 days for expansion. Then, the cells were split again and plated to be differentiated into neurons, when 80% confluency was reached (between 1 and 2 day). From about d13 to d33 of differentiation, the cells were cultured in 3N medium supplemented with 10 ng/ml BDNF (R&D systems, 248-BD-025, United States), 200 μM ascorbic acid (AA) (Sigma Aldrich, A4034, St. Louis, MO, United States), 10 ng/ml GDNF (R&D Systems, 212-GD-010, United States), 10 ng/ml NGF (R&D Systems, 256-GF-100, United States), 10 ng/ml NT-3 (R&D Systems, 267-N3-025, United States) and 0.5 mM cAMP (Sigma-Aldrich, D0260-100MG, St. Louis, MO, United States), resulting in the formation of immature PSNs. After these 35 days, cells could be replated after enzymatic passage with Accutase (Thermo Fisher Scientific, Waltham, MA, United States) and incubated with the same neuronal medium supplemented with 10 μM iRock (Millipore). Subsequently, after 1–2 days, neurons were kept for additional 2, 5, or 10 days in CM from neonatal human epidermal keratinocytes (HEKn), obtained as described in the following session.

Peripheral sensory neurons (PSNs) differentiation protocol time course. Schematic model of the differentiation protocol to generate sensory neurons from hiPSCs. Day 0 (d0): differentiation starts with hiPSCs presenting all main pluripotency markers and being replated to initiate differentiation; d5: formation of neural tube-like structures; d10: cell commitment to neural crest stem cell (NCPC) phenotype; d13: replacement to 3N media supplemented with brain derived neurotrophic factor (BDNF), ascorbic acid (AA), glial derived neurotrophic factor (GDNF), neurotrophin 3 (NT-3), and cyclic AMP (cAMP) to induce NCPC differentiation into sensory neurons; d13-d33: maturation of sensory neurons. d34-d44: treatment with conditioned medium (CM) from human neonatal epidermal keratinocytes (HEKn).

HEKn-Conditioned Medium

Plates were previously treated with gelatin (Sigma, St. Louis, MO, United States) and HEKn (Thermo Fisher Scientific, Waltham, MA, United States) were cultured in Epilife medium (Life Technologies, Rockville, MD, United States). Cells were split when 70% confluence was reached. After 48 h conditioning, medium was collected, centrifuged to remove debris and dead cells and was added fresh to the neurons at 75% with 25% of neuronal media.

Immunocytochemistry

Neural crest stem cells and sensory neurons were cultured in 96-well plates and fixed with 4% paraformaldehyde, permeabilized with Triton X-100 (0.3%) and blocked with 3% fatty acid-free bovine serum albumin (BSA) (Sigma, St. Louis, MO, United States). Cells were incubated for 2 h with primary antibodies diluted in 3% BSA. After washing with PBS, fluorescently labeled secondary antibodies were added for 40 min in the dark, washed thoroughly with PBS followed by a 5-min incubation with DAPI (4′,6-diamidino-2-phenylindole) for nuclear staining. After rinsing with PBS and water, 50 μl of glycerol were added as mounting medium in each well and the plates were sealed with aluminum sticker before analysis. The primary antibodies used were directed against: nestin (1:100, Sigma, St. Louis, MO, United States), β-tubulin Class III (TuJ1, 1:200, Merck-Millipore, Merck Millipore, Darmstadt, Germany), Islet1 (1:1000, Abcam, Cambridge, United Kingdom), TRPV1 (1:1000, Abcam, Cambridge, United Kingdom), BRN3A (1:250, Abcam, Cambridge, United Kingdom) and peripherin (1:250, Santa Cruz Biotechnology, Dallas, TX, United States). Secondary antibodies conjugated with Alexa Fluor 488 and Alexa Fluor 594 (1:400, Life Technologies, Rockville, MD, United States). Images were acquired with a High-Content Screening microscope, Operetta (PerkinElmer, Waltham, MA, United States) using a high numerical aperture 20x objective lens and specific filters for fluorescence excitation and emission. Image analyses were performed using high-content image analysis software Harmony 5.1 (PerkinElmer, Waltham, MA, United States) for morphology analysis and fluorescence quantification. Briefly, for quantification of immunostaining and cell counting, cells were segmented by proprietary algorithm (PerkinElmer, Waltham, MA, United States) in a constitutive marker channel (e.g., TuJ1 for neurons or Nestin for NCPCs and DAPI for nuclei). The segmented areas in each image constituted a template used to measure fluorescence in the same image but in the channels related to the staining of interest. Three independent experiments were performed for immunofluorescence imaging and quantifications. Each marker was probed in triplicate in each experiment and twenty different unbiased microscopic fields were imaged in each replicate. Neurite length measurements were performed using proprietary algorithms as described above, however, TuJ1 staining was used as the primary marker for identifying, segmenting and measuring neuronal projections. The number of cell in each image was counted with DAPI staining and the ratio of TuJ1+ cells and DAPI stained nuclei was calculated. Fluorescence intensity measurements for Islet1 were performed using DAPI staining as the constitutive marker. Raw data of immunostaining quantifications are available on Supplementary Datasheet of the Supplementary Material. Negative controls for immunocytochemistry were obtained by withholding the primary antibodies and resulted in no staining (Supplementary Figures ).

Electrophysiology

Whole-cell patch clamp recordings were performed on neurons plated on PLO and laminin-treated coverslips, as previously described (Acquarone et al., 2015). The coverslip was transferred to a 0.5 ml chamber continuously perfused with 165 mM NaCl, 5 mM KCl, 2 mM CaCl2, 10 mM dextrose, 5 mM HEPES, 2 mM NaOH, pH 7.35, at room temperature (23°C). The microelectrode (intracellular) solution contained 120 mM K-gluconate, 45 mM KCl, 0.5 mM CaCl2, 5 mM EGTA, 10 mM HEPES, 2 mM ATP-Mg, 1 mM MgCl2, pH 7.35. Electrode resistance was 2.9–3.7 MΩ and corrections were applied for the estimated -13 mV liquid junction potential. Cells randomly selected for recordings had clear neuronal morphology, with 2–3 thin processes emerging from the cell body, a visible nucleolus and smooth surface under Hoffmann modulation contrast optics. Cells were first tested under current clamp for action potentials evoked by 500 ms depolarizing currents steps. Next, under voltage clamp, membrane potential was held at -70 mV and voltage-sensitive currents were evoked by 50 ms depolarizing pulses in 10 mV increments from -50 to +20 mV, preceded by a 20 ms hyperpolarization to -90 mV. Leak currents were subtracted with a P/4 protocol. Data were recorded and analyzed using WinWCP software, available at http://spider.science.strath.ac.uk/sipbs/software_ses.htm (John Dempster, University of Strathclyde, United Kingdom).

Calcium Assays

Neurons were plated onto PLO and Laminin-treated coverslips where they were treated with a HEKn CM for 10 days. On the day of the experiment, the medium was removed and cells were incubated with fresh medium (without factors) containing 5 μM Fura-2 AM (Molecular Probes, Eugene, OR, United States) and 0.04% Pluronic acid for 1 h at 37°C. This solution was then replaced by fluorimetry buffer, composed of (in mM): 145 NaCl, 5 KCl, 1.2 NaHPO4, 1.5 CaCl2, 1 MgCl2, 10 D-glucose, 5 HEPES, pH 7.4.

After 30 min at 37°C, the coverslips were placed in an Eclipse Ti inverted microscope (Nikon, Minato, Tokyo, Japan) equipped with an Evolve 512 EMCCD camera (Photometrics, Tucson, AZ, United States), a Lambda DG-4 light source (Sutter Instrument, Novato, CA, United States) and a fluid flow system. The images were acquired every 500 ms via MetaFluor software (Molecular Probes, Eugene, OR, United States) and the levels of emitted fluorescence measured at 510 nm after stimulation with 340 and 380 nm were converted to a ratio (340/380 nm) and represented in pseudocolor. All drugs were freshly prepared from stocks and diluted in the fluorimetry buffer, with the control stimulation containing the vehicle (0.1% ethanol).

ELISA

Substance P (SP) levels were determined using a commercially available ELISA kit (Cayman, Pickerington, OH, United States), following the manufacturer’s instructions. Sensory neurons were plated in 24-well plates at high confluence. After 24 h, cells were rinsed with Hanks’ balanced salt solution (HBSS) buffer, composed of (in mM): 145 NaCl, 5 KCl, 1.2 NaHPO4, 1.5 CaCl2, 1 MgCl2, 10 D-glucose, 5 HEPES, pH 7.4. Following washing, the cells were incubated with different concentrations of TRPV1 agonists, other agents, or 45 min in the same buffer. When used, antagonists were applied 15 min before the agonist addition. Subsequently, the supernatants were collected and immediately assayed with the ELISA kit, following the manufacturer’s instructions. All experiments were performed in triplicate and absorbance levels were measured in a microplate reader.

Anandamide (AEA), rimonabant (SR141716A) and resiniferatoxin were obtained from Cayman Chemical, Co. (Ann Arbor, MI, United States) and hydrogen peroxide from Calbiochem.

Transcriptomic Analysis

Total RNA from three independent cultures of hiPSC-derived neurons, matured in the absence (C) or presence of human epidermal keratinocyte conditioned medium (CM) were extracted using QIAGEN miRNAasy kit (Invitrogen, Carlsbad, CA, United States), following the manufacturer’s instructions. RNA concentrations were measured by spectrophotometry using Nanodrop (Thermo Fisher Scientific, Waltham, MA, United States). The integrity of the RNA was analyzed using Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States) and the RNA integrity number (RIN) was obtained by identifying the ribosomal subunits in the samples. Only the samples with RIN greater than eight were sequenced, with an average coverage of more than 50 million reads per library. RNA sequencing was performed at Life Sciences Core Facility (LaCTAD, State University of Campinas, São Paulo, Brazil) using HiSeq 2500 platform (Illumina) in a paired-end mode (Supplementary Table ).

Quality analysis was performed with Fastqc and Trimmomatic 0.36. Reads passed in the quality control were used to calculate transcript expression profile with the software Kallisto 0.43.1 (Bray et al., 2016), using bootstrap 30 and kmer-size 31 with ENSEMBL GRC38 Homo sapiens release 88 transcriptome. Kallisto output was used as input in the software Sleuth (Pimentel et al., 2017) to perform the differentially expressed genes (DEG) analysis. Only transcripts with an adjusted p-value < 0.01, using Benjamini–Hochberg method, were considered significant. The Venn diagrams were generated using the Interactive Venn tool available at http://www.interactivenn.net/ (Heberle et al., 2015). The transcriptional profile of the PSNs were compared with Genotype-Tissue Expression (GTEx) public data using Heat∗seq tool (Devailly et al., 2016) available at http://www.heatstarseq.roslin.ed.ac.uk. Correlation heatmap was generated from average RNA-Seq data obtained from six biological replicates of PSNs cultures, comparing with gene expression levels of neuronal and skin tissues. Functional protein association network of hiPSC-derived PSN expressed genes was built using STRING version 10.5, available at https://string-db.org/ (Szklarczyk et al., 2015).

Statistical Analysis

Graph results are presented as average ± SEM, and differences were tested for significance as indicated in the figure legends. All statistical tests were performed with GraphPad Prism 6. For transcriptomic evaluation, DEGs analyses were performed using Sleuth (Pimentel et al., 2017), with qval output equivalent to false discovery rate (FDR) with adjusted p-value using the Benjamini–Hochberg test.

Results

Neural Crest Progenitor Cells (NCPCs) Differentiation From Human iPS Cells

Human induced pluripotent stem cells were differentiated into NCPCs after Smad inhibitors treatment for 10 days (d10, Figure ). Cells were then dissociated, cultured in 96-well plates and characterized for the expression of specific markers. NCPCs were positive for Nestin (Figures ), TRPV1 (Figure ), Peripherin (Figure ), BRN3A (Figure ), Pax6 (Figure ) (Zhang et al., 2010), and negative for Islet1 (Figure ), merged images shown in Figures . The expression of these markers was quantified as shown in Figure . Similar results were also obtained with a different hiPSCs line (Supplementary Figure ; Paulsen et al., 2012).

Characterization of NCPCs with specific markers. On day 10, cells presented neural progenitor cell’s morphology and positive staining for Nestin (A,D,G,J,M), TRPV1 (B), Peripherin (E), BRN3A (H), Pax6 (K), and negative for Islet1 (N). Nuclei were stained with DAPI and superimposed with staining for nestin and TRPV1 (C), Peripherin (F), BRN3A (I), Pax6 (L), and Islet1 (O). Calibration bar = 100 μm. Quantification of positive cells for each marker relative to DAPI-stained nuclei, n = 2–3 independent experiments (P).

Neural Induction

On day 10, NCPCs growth medium was replaced with fibroblast growth factor (FGF-) and epidermal growth factor (EGF-) containing medium for 2 days. Subsequently, cells were plated onto Poly-Ornithine (PLO)/Laminin-coated plates. Cells were maintained in 3N medium with addition of brain derived neurotrophic factor (BDNF), AA, glial derived neurotrophic factor (GDNF), nerve growth factor (NGF), neurotrophin-3 (NT-3) and cyclic AMP (cAMP), resulting in the formation of PSNs (Figure ). Of note, neurons tended to form ganglion-like structures after 7 days in 3N medium, which is also described for mouse and rat primary cultures (Sweetnam et al., 1982). On day 33 (d33) of differentiation, cells were enzymatically detached and replated onto 96-well plates.

Keratinocytes are important in the PSNs maturation process (Mantyh et al., 2011), thus, HEKn CM was added to neuronal culture at 75% final concentration for 2, 5, and 10 days before experimental assays. Neuronal maturation was quantified by β-tubulin III staining at these different time points. The neurites emanating from the ganglion-like structures increased robustly over this time, increasing approximately 157% in length after 5 days in culture and 542% after 10 days (Figures ). The presence of HEKn CM accounts for a 17.4% increase in cell differentiation, when compared to control cultures maintained in regular 3N medium (Figure ).

Sensory neurons differentiation. Beta-tubulin class III (green)-stained neurons after 2, 5, or 10 days in neuronal medium conditioned with 75% HEKn medium. Nuclei were stained with DAPI. (A–C) Quantification of neurite length in PSNs was based on beta-tubulin Class III (green) immunostaining. Computer assisted image segmentation identified all neurites and measured the total length at each time point. Neurites increased by 328% in length after 5 days and 514% after 10 days in culture with 75% HEKn CM (D). One-way ANOVA followed by Tukey’s multiple comparison test, ∗p < 0.05 and ∗∗p < 0.01, n = 3 independent experiments individual experiments from three separate differentiations. Percentage of differentiated neurons was increased by HEKn-CM treatment by 17.4% after 10 days of culture, as measured by β-tubulin III-staining (E). Unpaired t-test, ∗∗∗p < 0.001, n = 3 individual experiments from three separate differentiations. Scale bar = 100 μm.

Neuronal induction was further verified by electrophysiological recordings (Figure ). Out of 10 cells from a single differentiation batch and sampled from different coverslips, nine fired at least one action potential under current clamp and all showed fast-inactivating inward currents under voltage clamp, typical of voltage-dependent sodium channels (mean maximal amplitude: -537.2 ± 57.9 pA). Slowly-activating, non-inactivating potassium outward currents ranged from 0 to 1504 pA in response to a depolarizing step from -90 to 20 mV (mean: 611.2 ± 160.2 pA).

Patch-clamp recordings in hiPSC-derived neurons. (A) Photomicrography illustrating the patched neuron whose action potentials are shown in (B). (B) Action potentials fired during a 30 pA positive current injection (baseline: –70 mV) and (C) inward and outward voltage-gated currents evoked in the same representative cell in response to incrementing depolarizations from –90 mV to –40, –30, –20, –10, 0, 10, and 20 mV.

There are several markers that make up the distinctive expression profile of somatosensory neurons. We assessed the expression of TRPV1 (Figures ), Peripherin (Figures ), and Islet1 (Figures ). All markers showed increased expression over the observed time points (Figures ). TRPV1 increase was the most noteworthy, enhancing approximately five times. These results suggest that CM promoted the maturation and growth of sensory neurons, as they exhibited higher expression of differentiation markers, hallmarks of the transition into sensory neurons maturity. The expression of all these markers, except Islet1, was also detected among RNA-Seq transcripts (see Figure ).

Peripheral sensory neurons maturation. Neurons cultivated in 75% HEKn cells media presented time dependent increase in expression of TRPV1 (A,C,E,G), Peripherin (H,J,L,N) and Islet-1 (O,Q,S,U). Beta-tubulin class III (green) and DAPI (blue) were co-stained in all conditions (B,D,F,I,L,K,M,P,R,T). TRPV1 fluorescence intensity increased by 328 and 514% after 5 and 10 days, respectively (G); peripherin showed 118 and 165% increase after 5 and 10 days, respectively (N) and Islet-1 also presented 205 and 276% increase after 5 and 10 days, respectively (U). One-way ANOVA followed by Tukey’s multiple comparison test, ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗∗p < 0.0001, n = 3 individual experiments from three separate differentiations.

Activation of Human PSNs by Nociceptive Agents

The TRPV1 channel is a non-selective cation-permeable channel, with significant calcium permeability (Szallasi et al., 2007). For this reason, intracellular calcium measurements have been used as a proxy to TRPV1 activation. Here, we attempted to quantify the increases in calcium mediated by TRPV1 activation using dynamic images of calcium transients in individual cells. In these conditions, capsaicin-induced activity was poorly detected (Supplementary Figure ). Moreover, depolarization of neurons induced by a high potassium solution also showed weak calcium responses (Supplementary Figure ). These results suggest that in terms of intracellular calcium increases, these neurons have not fully matured.

Since other authors have also observed very few calcium responses in hiPSCs-derived PSNs, attributing this lack of response to neurite activity rather than cell soma (Chambers et al., 2012), we decided to utilize another readout of functional activity. PSNs respond to irritants, such as capsaicin, by releasing neuropeptides at the epidermal end and glutamate and neuropeptides in the dorsal horn of the spinal cord (Basbaum et al., 2009). To verify whether the sensory neurons were functional with respect to neuropeptide release, after d33 they were plated onto 24-well plates and treated 10 days with HEKn CM, as described above. Basal Substance P (SP) release was approximately 3 pg/mL/h. This amount was tripled after 300 nM Capsaicin treatment; however, this increase was not statistically significant. A higher concentration of capsaicin was applied (3 μM) but did not lead to an increased release of SP, probably due to TRPV1 desensitization (data not shown). Nevertheless, 300 nM resiniferatoxin effectively increased SP release, as did 20 nM bradykinin (Figure ). Interestingly, 100 μM anandamide evoked the highest SP release by the tested TRPV1 agonists with these cells. Moreover, when a high potassium solution (70 mM) was added to depolarize neurons, we obtained an increase in substance P release, though smaller then resiniferatoxin, bradykinin, and anandamide (Figure ). In addition, the cells responded to 0.25% hydrogen peroxide by releasing the highest amount of SP.

Substance P release after in vitro treatment with nociceptive agents. (A) Substance P release measured in the supernatant of sensory neurons culture after incubating for 45 min at 37°C in the presence of ligands. CAP: capsaicin (300 nM); RES, resiniferatoxin (300 nM); AEA, anandamide (100 μM, unless stated otherwise); BK, bradykinin (20 nM); KCl, high potassium solution (70 mM) and H2O2 (0.25%). (B) Concentration-dependent effect of anandamide in substance P release. (C) Anandamide effect (30 μM) is partially inhibited by capsazepine (CPZ) (30 μM). (D) Anandamide effect (100 μM) is not significantly changed by rimonabant (SR) (10 μM). (E) Transcript levels of the putative receptors involved in substance P release by the tested agonists. TRPV1, transient receptor potential vanilloid 1; CB1, cannabinoid receptor 1; CB2, cannabinoid receptor 2; BDKRD1, bradykinin receptor 1; BDKRD2, bradykinin receptor 2. One way ANOVA followed by Tukey’s multiple comparison test, ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗∗p < 0.0001. n = 2–3 individual experiments from three separate differentiations.

The anandamide effect, albeit at relatively high concentrations, was concentration-dependent in a range compatible with TRPV1 agonism rather than CB1 activation (Figure ). Furthermore, AEA-increase of SP release was partially inhibited by capsazepine, a TRPV1-selective antagonist (Figure ). Rimonabant, a CB1-selective antagonist did not block the AEA effect, on the contrary, it tended to increase it (Figure ). These results indicate that functional nociceptive responses are present in these neurons, particularly TRPV1-mediated ones. The expression of the potential receptors involved in these functional responses are plotted in Figure and include TRPV1, CB1 and CB2, and both bradykinin receptors, BDKRD1 and BDKRD2. Interestingly, CB2 was detected in negligible amounts.

Transcriptomic Analyses

RNA-Seq was used to determine the transcriptional profile of the cells obtained using the differentiation protocol aforementioned. RNA from control PSN (PSN-C) and treated with HEKn-conditioned media (PSN-CM) during maturation in culture have both shown expression of all main expected neuronal markers (Figure ). Transcripts coding the same markers that were detected by immunostaining, such as Nestin (NES), Pax6 (PAX6), Peripherin (PRHP), BRN3A (POU4F1). TRPV1 was the most expressed TRP channel of sensory neurons, confirming the immunofluorescence staining. Also, the presence of the most expressed neuronal marker genes like PIEZO2 (piezo type mechanosensitive ion channel component 2), TH (tyrosine hydroxylase), TUBB3 (tubulin beta 3 class III), Neurotrophic receptor tyrosinase 2 (NTRK2, also known as TrkB), SCN3A and SCN3B (Nav 1.3 or sodium voltage-gated channel alpha and beta subunit 3) suggest that the obtained cell type has a C-low threshold mechanoreceptor (C-LTMR) profile, EGF (epidermal growth factor), RUNX1T1 (RUNX1 translocation partner 1), CALCB (CGRP receptor component), LDHB (lactate dehydrogenase B), MDK (midkine, neurite growth-promoting factor 2), MSN (moesin) and P2RX3, P2RX4 (purinergic receptor P2X ligand-gated ion channel 3 and 4), were found at different levels in both PSN-C and PSN-CM transcriptional fingerprints, in three independent cultures per condition (Figure and Table ). RNA-Seq global expression data is available on Supplementary Datasheet of the Supplementary Material.

Neuronal markers expressed and regulated in hiPSC-derived PSNs. (A) Heatmap (blue: low expression, red: highly expressed) of transcripts related to neuronal marker genes expressed in cultures of control PSN (PSN-C) and PSN treated with HEK-conditioned media (PSN-CM), obtained from RNA-Seq data. Transcripts coding marker genes Nestin (NES), Peripherin (PRHP), BRN3a (POU4F1). Moreover, mRNA encoding the cannabinoid receptor 1 (CNR1) were found in higher levels compared to TRPV1, which was also expressed but in lower levels. PIEZO2 (piezo type mechanosensitive ion channel component 2), TH (tyrosine hydroxylase), TUBB3 (tubulin beta 3 class III), Neurotrophic receptor tyrosinase 2 (NTRK2, also known as TrkB), SCN3A and SCN3B (Nav 1.3 or sodium voltage-gated channel alpha and beta subunit 3), EGF (epidermal growth factor), RUNX1T1 (RUNX1 translocation partner 1), CRCP (CGRP receptor component), LDHB (lactate dehydrogenase B), MDK (midkine, neurite growth-promoting factor 2), MSN (moesin) and P2RX3, P2RX4 (purinergic receptor P2X ligand-gated ion channel 3 and 4) were highly expressed in all cultures. (B) Venn diagrams were generated using InteractiVenn tool showing unique and shared genes expressed in PSN-C and PSN-CM (treated with HEKn-CM).

Venn diagrams (Figure ) show about 10,325 genes significantly expressed (TPM > 10), with 9,610 genes expressed in both conditions, 266 expressed only in PSN-C and 449 genes exclusively expressed in PSN-CM condition, which could be interesting targets to be investigated to increase PSN maturation. Most of the transcripts that changed were downregulated after the treatment with HEK-CM (Supplementary Figure ). Moreover, mRNA encoding the cannabinoid receptor 1 (CNR1) was found in higher levels compared to TRPV1 (Figure ). This could mean that the AEA-induced SP release via TRPV1 could be underestimated because CB1 is also stimulated by anandamide and might inhibit TRPV1. In fact, neurons pretreated with rimonabant, a CB1-selective antagonist, show a tendency to increase SP release with anandamide (Figure ). Interleukin 1 α and β (IL1α, IL1β) transcripts were found only in PSN-CM cultures. Principal component analysis (PCA) was performed using ClustVis and provided an accurate verification of the applied method, as the first principal component clearly separated the samples by condition (C vs. CM) and the second principal component separated the samples by independent PSN cultures differentiation (Supplementary Figure ).

The transcriptome comparison between PSN and GTEx RNA-Seq dataset from human tissues shows higher scaled correlation with brain tissues like hypothalamus, amygdala, frontal cortex, hippocampus (scaled correlation 0.950; 0.936; 0.935, and 0.935 respectively) (Figure ), suggesting that PSN present a transcriptional profile closer to neuronal than other human tissues or skin cells, like fibroblasts, which were the source of the cells that were reprogrammed to generate the hiPSC used in this study.

Correlation heatmap comparing PSN transcriptional profile with human GTEx gene expression data of neuronal and skin tissues. Heatmap was generated using Heat∗seq tool (Devailly et al., 2016). uscal replicates of PSNs cultures, showing higher scaled correlation with neuronal human tissues such as hypothalamus, amygdala, frontal cortex, hippocampus (scaled correlation 0.950; 0.936; 0.935, and 0.935 respectively).

With the RNA-Seq it is possible to identify transcripts uniquely expressed in neuronal subtypes, directly related to functionality and the identification of neuronal types. The relevance of these markers was previously validated in vivo with identified population-specific genes in high fraction of positive cells (Usoskin et al., 2015). Neuronal population expressing PIEZO2 and an exclusive C-LTMRs Slc17a8, codifying the vesicular glutamate transporter 3 (Vglut3), involved in mechanical pain and pleasant touch (Seal et al., 2009; Coste et al., 2010; Li et al., 2011; Abraira and Ginty, 2013).

The transcriptomic analysis confirmed the expression of the main neuronal markers expected for somatossensorial DRG neurons, with predominantly transcriptional profile compatible with C-LTMR, at least with the available data from rodents, as showed in Supplementary Figure .

Discussion

According to literature, several groups have reported generation of PSNs from pluripotent stem cells. However, none was able to show robust TRPV1 activity. Chambers et al. (2012) showed TRPV1 expression via quantitative PCR and correlated with only 1–2% of the cells responding to capsaicin. Even the responding cells seem to exhibit a very low activity. Others showed electrophysiological activity of human sensory neurons directly reprogrammed from fibroblasts, but not capsaicin-elicited signals (Wainger et al., 2015). In the same work, mouse neurons obtained with a similar method showed robust TRPV1 activity. One could speculate that the authors were not able to detect TRPV1 activity in human neurons consistently as well. Another work reported detecting TRPV1 activity induced by capsaicin after 6 weeks in media containing growth factors but did not further describe this response (Young et al., 2014). PSNs generated by other groups that did not demonstrate TRPV1 activity but focused on sodium currents and action potentials (Menendez et al., 2011; Lee et al., 2012; Reinhardt et al., 2013). To our knowledge, this is the first report of human PSN obtained by any method shown to release SP in response to nocifensive agents.

TRPV1 was detected by immunostaining in ∼90% of neurons, but the RNA levels and substance P release induced by TRPV1 agonists were relatively small. One hypothesis to explain this discrepancy between the TRPV1 staining and function is that, although present, TRPV1 does not reach full activity due to modulation of the presence of this channel on the cell surface, via interaction with other proteins (Supplementary Figure ). For instance, CB1 (CNR1) is overexpressed (Figure ) relatively to TRPV1. This overexpression could at least partially explain the low detectable activity of TRPV1, as CB1 is known to negatively modulate this channel through dephosphorylation (Supplementary Figure ) (Ross, 2003; Yang et al., 2013). Accordingly, acute pretreatment with a selective CB1 inverse agonist tended to augment AEA-induced SP release, presumably through reduced inhibition of TRPV1. However, it remains to be established whether chronic suppression of CB1 activity could effectively enhance TRPV1 function.

In addition, in the case of sensory neurons expressing TRPV1 channels, these proteins are packaged into CGRP- or SP- and VAMP1-containing vesicles and delivered to the plasma membrane, involving the formation of SNARE complexes composed of SNAP-25, syntaxin 1, and VAMP1, as well as Munc18–1 (Meng et al., 2016). This delivery can be enhanced by several factors, including TNFα. Since SP release was modestly elicited by capsaicin, it is possible that this aspect of nociceptor development is not completely mature at the differentiation level presented here and therefore one could speculate that TRPV1 is not adequately being shipped to the cell surface.

Furthermore, other proteins known to negatively modulate TRPV1 function that actually are binding partners are calmodulins, which are very highly expressed (data not shown and Supplementary Figure ) (Ho et al., 2012). Again, by negatively modulating the levels of these proteins, a more robust TRPV1 activity could arise.

We have shown that, although a negligible number of cells respond to capsaicin measured by intracellular calcium increases, we could detect substance P release under similar conditions. Release of SP or any other neuropeptide is thought to be dependent of intracellular calcium increases, which would allow the proper fusion of synaptic vesicles to the presynaptic membrane thereby releasing its contents in the synaptic cleft. These two findings would therefore seem irreconcilable. However, Purkiss et al. (2000) proposed that there are two mechanisms by which SP is released by rat nociceptors: one being calcium dependent and the other independent. It is possible, therefore, that the calcium-independent mechanism is responsible for the observed SP release. Alternatively, calcium increases might have been predominant in neurites rather than in cell bodies, as observed by Chambers et al. (2012), which is harder to detect with the imaging technique.

The fact that Young et al. (2014) only observed TRPV1 activity after 6 weeks in the presence of growth factors corroborates this notion. However, maintaining cells for another 2 weeks in vitro does not guarantee that this activity would be satisfactory in the end, as mentioned. Therefore, we must consider additional mechanisms to enrich the neuronal population in TRPV1+/TrKA+/CGRP+ cells. First, the amount of NGF necessary might be higher, or perhaps a gradient of this neurotrophin might be necessary. NGF is present in high concentrations within the epidermis, where nerve fibers with free endings (i.e., nociceptor fibers) will project (Davies et al., 1987). In murine models, it is known that several steps are necessary to differentiate TrkA+/Met+/CGRP+ neurons, the ones with putative TRPV1 activity (Lallemend and Ernfors, 2012). Our transcriptomic results show a high expression of TrkB, which was reported to be present at high levels in a subset of murine medium-diameter DRG neurons (Li et al., 2011). It was unexpected to find high expression of TrkB and expression of TRPV1, based on extensive rodent data. However, we cannot rule out these transcripts are present in different cells. In addition, previous works have compared human and mouse DRG gene expression and have found an overlap of approximately 50–70% (Parisien et al., 2017; Ray et al., 2018). Consequently, it is possible that not all transcripts found in mouse DRG neurons will also be found in humans’.

Although individual cell gene expression variability is influenced by differentiation conditions, our PSN expressed the five sensory-neuronal genes marker genes, except SCN9A (Schwartzentruber et al., 2017). TH, which catalyzes the production of L-DOPA from tyrosine in the catecholamine biosynthesis pathway, is a defining feature of low threshold mechanoreceptors (LTMRs) in adult DRGs (Brumovsky et al., 2006) and was found to be also highly expressed in hiPSCs-generated PSN. The presence of many individual transcripts characteristic of DRG and somatosensory neurons was confirmed. These include RET, TH, PRPH, and LDHB. However, some specific LTMR transcripts were not detected, such as VGLUT3, Nav 1.8 and TRPA1. PIEZO2 and GFRA2, markers of LTMR were expressed and P2RX3, a purinergic ion channel expressed in C-fiber non-peptidergic neurons was also detected (Yin et al., 2016) as shown in Table . CACNA1H was highly expressed, being a voltage-gated calcium channel (CaV) subtype 3.2, uniquely expressed in unmyelinated C-LTMR (Reynders et al., 2015). Together with CACNAH1, the high expression of TrkB, TH, and LDHB suggest the PSN generated here have a LTMRs profile of peripheral neurons involved in sensory perception (Usoskin et al., 2015). Our RNA-Seq results suggest that the hiPSC-derived PSN present a transcriptional profile compatible with C-LTMR predominantly (Supplementary Figure ). However, to precisely define the proportion of this neuronal type in the culture further investigations have to be performed using single cell RNA sequencing.

In summary, although some reports have described the generation of human PSNs, none have successfully demonstrated a robust and useful TRPV1 activation. The present work shows detectable TRPV1 activity, the release of SP mediated by resiniferatoxin and anandamide that could be used for screening of nociceptive agents and possibly analgesics. The transcriptional data obtained from RNA-Seq show the main neuronal marker genes expressed, in different levels, in independent PSN cultures. The strategy of treating PSN with CM to investigate whether it could increase cell maturity and TRPV1 functionality had opposed outcomes. Although neurite growth increased, which is an indication of neuronal maturity, there were no conclusive changes in relevant neuronal markers expression in RNA-Seq data, although TRPV1 immunostaining increased. Nevertheless, as TRPV1 activity would be useful for the purpose of screening for new nociceptive/irritant agents, it is appealing to further improve the robustness of this activity and to develop reliable in vitro models using these differentiated human neurons.

Data Availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author Contributions

MG, RDV, JS, BP, RdC, NC, LB, and SR conceived and designed the experiments. MG, RDV, GV, JS, BP, IL, RdC, and NC performed the experiments. MG, RDV, BP, RdC, FRdS, and NC analyzed the data. MG, RDV, BP, RdC, FRdS, NC, LB, and SR wrote the paper.

Conflict of Interest Statement

LB and RDV are employees of L’Oréal Research & Innovation. The remaining authors have received sponsored research support from L’Oréal Research & Innovation.

Table 1

RNA-Seq data of the main sensory neuronal markers expression in transcripts per million (TPM).

Gene	Full gene description	C1	C2	C3	CM1	CM2	CM3	Mean
LDHB	Lactate dehydrogenase B	938.44	882.70	908.97	897.17	901.21	890.95	903.24
MDK	Midkine, neurite growth-promoting factor 2	125.85	114.15	119.24	122.63	111.19	137.75	121.80
TH	Tyrosine hydroxylase	20.45	24.49	30.11	25.72	22.07	25.99	24.81
CACNA1A	Calcium voltage-gated channel subunit alpha1 a	12.94	12.77	14.63	10.42	9.08	9.96	11.63
NTRK1	Neurotrophic receptor tyrosine kinase 1 (trka)	0.39	0.57	0.25	0.26	0.41	0.46	0.39
NTRK2	Neurotrophic receptor tyrosine kinase 2 (trkb)	876.38	929.35	842.48	802.70	846.97	832.39	855.05
NTRK3	Neurotrophic receptor tyrosine kinase 3 (trkc)	25.82	25.30	26.50	21.46	24.21	22.77	24.35
RET	Ret proto-oncogene	2.26	2.06	2.24	2.00	2.82	2.29	2.28
G-protein coupled receptors and associated proteins
GNB4	Guanine nucleotide binding protein (g protein), beta 4	41.82	53.73	37.25	43.15	35.51	41.15	42.10
GNA14	Guanine nucleotide binding protein, alpha 14	0.17	0.16	0.21	0.17	0.13	0.14	0.16
MRGPRD	mas-related gpr, member d	0.03	0.00	0.00	0.00	0.00	0.00	0.01
MRGPRX1	mas-related gpr, member x1	0.00	0.00	0.00	0.00	0.00	0.00	0.00
NPY1R	Neuropeptide Y receptor Y1	1.49	1.71	1.19	1.27	0.88	1.21	1.29
GALR1	Galanin receptor 1	0.15	0.25	0.14	0.21	0.20	0.14	0.18
LPAR3	Lysophosphatidic acid receptor 3	0.17	0.21	0.30	2.87	2.28	2.52	1.39
RGS8	Regulator of G-protein signaling 8	16.46	16.30	18.15	12.08	10.77	15.40	14.86
GPR35	G protein-coupled receptor 35	1.39	2.12	1.91	2.01	1.82	2.15	1.90
LPAR5	Lysophosphatidic acid receptor 5	0.14	0.12	0.11	1.07	0.98	0.79	0.53
SSTR2	Somatostatin receptor 2	13.29	13.54	12.83	9.08	8.06	8.67	10.91
PTGDR	Prostaglandin D receptor	0.00	0.02	0.00	0.00	0.00	0.00	0.00
BDKRB1	Bradykinin receptor B1	1.57	2.30	1.14	1.89	2.29	2.66	0.00
BDKRB2	Bradykinin receptor B2	0.93	1.02	0.95	1.68	1.93	1.92	0.00
CNR1	Cannabinoid receptor 1	20.26	22.07	20.90	14.09	15.79	15.69	0.00
CNR2	Cannabinoid receptor 2	0.10	0.07	0.06	0.13	0.13	0.03	0.00
Hox genes transcription factors involved in anteroposterior axis development
HOXA6	Homeobox A6	17.16	18.08	20.44	16.58	15.9	19.54	17.95
HOXA7	Homeobox A7	57.88	52.75	61.19	48.32	49.54	47.9	52.93
HOXA9	Homeobox A9	224.0	239.8	206.6	235.2	202.6	193.7	217.01
HOXA10	Homeobox A10	57.0	53.78	69.09	57.31	52.23	54.04	57.24
HOXB3	Homeobox B3	78.75	90.21	90.8	64.6	67.77	73.61	77.62
HOXB5	Homeobox B5	63.72	61.63	72.71	60.33	61.81	57.38	62.93
HOXB6	Homeobox B6	128.5	142.3	140.5	124.4	121.3	134.8	131.97
HOXB7	Homeobox B7	119.3	115.4	118.8	103	94.44	90.08	106.85
HOXC6	Homeobox C6	94.16	92.87	105.3	89.83	89.83	89.76	93.62
HOXC8	Homeobox C8	53.43	53.32	51.39	41.72	41.61	41.59	47.17
HOXC9	Homeobox C9	76.41	65.56	79.64	65.68	72.56	66.36	71.03
HOXC10	Homeobox C10	89.2	80.74	97.64	73.03	78.39	82.41	83.57
HOXD8	Homeobox D8	20.67	18.29	14.75	18.45	18.39	17.81	18.06
HOXD9	Homeobox D9	7.317	7.226	8.353	7.036	6.515	6.71	7.19
HOXD10	Homeobox D10	26.98	24.63	28.09	21.83	22.94	22.57	24.51
Ligand and voltage gated ion channels
P2RX3	Purinergic receptor P2X, ligand-gated ion channel, 3	77.21	77.75	85.73	52.13	64.43	67.41	70.77
SCN8A	Sodium channel, voltage-gated, type VIII, alpha; Nav1.6	6.32	6.87	7.38	5.11	5.04	5.93	6.32
SCN11A	Sodium channel, voltage-gated, type XI, alpha; Nav1.9	0.79	0.67	0.67	0.60	0.63	0.75	0.68
TRPM8	Transient receptor potential cation channel, subfamily M, 8	0.34	0.29	0.53	0.24	0.61	0.24	0.37
TRPA1	Transient receptor potential cation channel, subfamily A, 1	0.31	0.20	0.14	0.32	0.31	0.14	0.24
TRPV1	Transient receptor potential cation channel, subfamily V, 1	5.75	4.39	5.27	4.74	5.13	4.56	4.97
GRIK1	Glutamate receptor, ionotropic, kainate 1	4.98	4.73	4.86	4.23	3.46	3.83	4.35
KCNH6	Potassium voltage-gated channel, subfamily H (eag-related), 6	0.38	0.42	0.38	0.25	0.27	0.43	0.35
KCNF1	Potassium voltage-gated channel, subfamily F, member 1	10.54	11.22	14.08	10.15	12.84	12.12	11.83
TRPC3	Transient receptor potential cation channel, subfamily C, 3	0.36	0.46	0.50	0.55	0.51	0.66	0.51
CACNA1I	Calcium channel, voltage-dependent, alpha 1I subunit	4.53	4.22	4.71	3.43	3.84	3.69	4.07
Transporters and channel interacting proteins
CHRNB2	Cholinergic receptor, nicotinic, beta polypeptide 2	6.60	6.89	6.61	4.32	5.18	4.66	6.60
KCNIP2	a-Type potassium channel modulatory protein 2	24.03	23.33	23.93	16.01	16.91	18.71	20.49
PIEZO2	Piezo type mechanosensitive ion channel component 2	25.61	28.53	23.84	26.83	22.85	26.48	25.69
SCN3B	Sodium channel, voltage-gated, type III, beta	39.68	43.94	43.96	31.32	34.51	33.59	37.83
SLC37A1	Solute carrier family 37 (glycerol-3-phosphate transporter) 1	11.96	11.11	14.14	8.80	10.81	10.81	11.27
SLC4A11	Solute carrier family 4, sodium bicarbonate transporter-like 11	3.63	3.24	2.92	4.03	2.76	3.53	3.35
SLC17A8	Vesicular glutamate transporter 3 (VGluT3)	9.09	9.50	8.23	5.09	5.63	4.76	9.09
Growth factors, matrix and glycoproteins, semaphorins, and cytokine receptors and modulators
CADM1	Cell adhesion molecule 1; Cadm1	127.29	135.69	124.18	120.75	115.50	119.27	123.78
CD55	CD55 antigen; complement decay-accelerating factor	1.47	1.49	1.40	13.13	9.48	13.49	6.74
GFRA2	Glial cell line derived neurotrophic factor family receptor alpha 2	5.59	6.34	6.30	4.79	4.68	4.88	5.43
NELL1	NEL-like 1	2.96	3.28	3.86	2.70	2.23	2.19	2.87
Synaptic vesicle proteins and associated proteins
SYNPR	Synaptoporin	25.71	26.81	20.51	21.09	20.54	20.72	22.56
SYT7	Synaptotagmin VII	14.30	16.11	16.42	13.56	15.31	16.94	15.44
NRSN1	Neurensin 1	25.37	26.77	24.56	21.22	17.69	20.02	22.61
SYT9	Synaptotagmin IX	3.43	3.77	3.16	2.13	3.12	3.42	3.17
NPTX2	Neuronal pentraxin 2	23.62	23.30	26.38	17.35	17.76	17.09	20.92
ARHGAP26	Rho GTPase activating protein 26	6.38	7.37	6.93	5.76	7.79	6.92	6.86
Membrane, structural, and Rho/Rac proteins
RASGRP1	ras Guanyl releasing protein 1	2.64	1.90	2.64	2.15	2.24	1.84	2.23
PRPH	Peripherin	13.12	13.47	14.27	12.71	13.68	16.82	14.01
RAB27B	rab27b. member ras oncogene family	9.78	13.13	9.73	16.07	13.73	13.92	12.73
TUBA1A	Tubulin alpha 1a	5465.0	5242.5	5408.9	4256.1	4360.1	4470.4	4747.6
CDH9	Cadherin 9	13.38	15.50	13.63	11.26	11.00	11.36	12.69
Neuropeptides and cytokines
TAC1	Tachykinin peptide hormone family (substance P, neurokinin A)	1.26	1.35	1.93	1.37	0.57	1.12	1.26
CALCB	Calcitonin related polypeptide beta	1.74	1.67	2.30	1.94	1.06	1.19	1.65
ADCYAP1	Adenylate cyclase activating polypeptide 1; PACAP	46.45	49.57	50.18	41.51	47.80	49.25	47.46
Transcription factors, mRNA-binding proteins and transcription related and splicing factors
LDB2	LIM domain binding 2	4.17	3.07	3.16	6.84	4.27	6.60	4.69
CELF6	Bruno-like 6, RNA binding protein (Drosophila)	15.41	13.49	16.38	11.07	13.98	13.88	14.03
RUNX1	Runt related transcription factor 1	30.17	32.53	35.22	23.10	23.94	25.78	28.46
MYT1	Myelin transcription factor 1	14.40	15.93	15.19	13.56	13.37	14.13	14.43
PAQR5	Progestin and adipoQ receptor family member V	1.33	1.74	1.03	4.56	2.73	4.43	2.64
Kinases, phosphatases, proteases, Mmps, phosphodiesterases and other enzymes and enzyme inhibitors
PLCB3	Phospholipase C, beta 3	12.87	12.40	12.11	11.84	13.15	13.70	12.68
CAMK2A	Calcium/calmodulin-dependent protein kinase II alpha	7.10	8.11	8.09	6.10	7.13	6.50	7.17
PLCXD3	Phosphatidylinositol-specific phospholipase C, domain 3	1.16	1.84	1.49	0.96	1.02	1.17	1.27
DGKI	Diacylglycerol kinase, iota	9.99	10.21	11.17	7.47	7.23	9.45	9.25
HS6ST2	Heparan sulfate 6-O-sulfotransferase 2	137.20	160.55	127.81	116.13	108.99	112.14	127.14
RNF182	Ring finger protein 182	5.88	6.27	5.47	4.86	4.14	4.65	5.21
PDE11A	Phosphodiesterase 11A	1.03	0.96	0.96	1.02	1.04	0.91	0.98
Typical fibroblast expressed genes (Schwartzentruber et al., 2017)
CD63	Transmembrane 4 superfamily (tetraspanin family)	378.82	349.53	396.76	361.73	387.02	391.03	377.48
RPLP1	Ribosomal protein lateral stalk subunit P1	10.28	11.05	11.77	10.80	11.72	11.85	11.24
MSN	Moesin	98.13	98.60	97.33	89.68	97.73	101.48	97.16
VIM	Vimentin	2264.2	2191.6	2118.1	2248.6	2157.1	2292.1	2211.9
ITGB1	Integrin subunit beta 1	131.16	135.96	106.59	155.41	124.98	144.77	133.15