NOCICEPTRA: Gene and microRNA Signatures and Their Trajectories Characterizing Human iPSC‐Derived Nociceptor Maturation

Abstract Nociceptors are primary afferent neurons serving the reception of acute pain but also the transit into maladaptive pain disorders. Since native human nociceptors are hardly available for mechanistic functional research, and rodent models do not necessarily mirror human pathologies in all aspects, human induced pluripotent stem cell‐derived nociceptors (iDN) offer superior advantages as a human model system. Unbiased mRNA::microRNA co‐sequencing, immunofluorescence staining, and qPCR validations, reveal expression trajectories as well as miRNA target spaces throughout the transition of pluripotent cells into iDNs. mRNA and miRNA candidates emerge as regulatory hubs for neurite outgrowth, synapse development, and ion channel expression. The exploratory data analysis tool NOCICEPTRA is provided as a containerized platform to retrieve experimentally determined expression trajectories, and to query custom gene sets for pathway and disease enrichments. Querying NOCICEPTRA for marker genes of cortical neurogenesis reveals distinct similarities and differences for cortical and peripheral neurons. The platform provides a public domain neuroresource to exploit the entire data sets and explore miRNA and mRNA as hubs regulating human nociceptor differentiation and function.


DOI: 10.1002/advs.202102354
harm. Their cell bodies are located in the dorsal root (DRG) or trigeminal ganglia (TG) and their developmental fate is specified by gene expression networks governed by key genes including neurotrophic receptor tyrosine kinase A (TrkA), B (TrkB), and C (TrkC) together with Neurogenin 1 (Ngn1) and 2 (Ngn2). [1][2][3] Although comprehensive transcriptomic and proteomic analyses of mature DRGs have been performed in animal disease models and patient post-mortem tissues, human nociceptor expression signatures related to human pain disorders are scarce. [4][5][6][7][8] Notable differences between mice and men challenge the applicability of rodent model systems for the investigation of human pain pathogenesis. [9,10] Humancell-based model systems, such as induced pluripotent stem cells (iPSCs), which allow directed cell-type specific programming into virtually any cell type by applying specific differentiation protocols, [11][12][13][14][15] offer the opportunity to overcome this apparent translational paresis. Human iPSCs differentiated into nociceptor-like cells express fundamental nociceptor-like characteristics and emerge as a valuable tool to explore mechanisms underlying hereditary pain disorders, such as erythromelalgia or migraine, and even to develop and select personalized treatments. [16][17][18][19][20][21][22] However, the molecular equipment of iPSC derived nociceptors and its similarity to native human nociceptors are not yet sufficiently documented, which may limit the utilization of these for mechanistic studies.
Developmental gene expression signatures and their regulation are currently only available for iPSC-derived cortical neurons. [23,24] Protein-coding genes and long non-coding RNAs (lncRNAs) have been the focus of recent gene expression studies during cell differentiation. [25,26] At the same time, a critical role for a class of transcriptional regulators called microRNAs as hub regulators of pluripotency, balancing of cell fate decisions, as well as guiding of cell-cycle progression and self-renewal properties of human embryonic stem cells are emerging. [27][28][29][30][31][32][33][34] Specific miRNA expression patterns, time trajectories, and interactions with proposed target genes during nociceptor differentiation are so far unavailable but, as a decisive layer of regulatory complexity, may be highly important for neuronal cell fate decisions, morphology, target finding and finally function and dysfunction.
In order to tackle these issues, we performed unbiased mRNA::miRNA co-sequencing followed by analyses of timecourse specific gene expression signatures of human iPSCderived nociceptors (iDNs) at six different timepoints. We identified five specific sub-stages of nociceptor differentiation as well as their control by hub genes and miRNAs, highlighting the critical importance of mRNA and miRNA expression trajectories. Finally, we provide the integrative online platform NOCICEPTRA to both retrieve the presented experimental data as well as to query custom gene expression data based on specific disease and pathway enrichments during nociceptor differentiation.

Human Nociceptor Differentiation Involves Highly Correlated Gene Expression Trajectories
Temporal gene expression patterns during the differentiation of iDNs from iPSCs were assessed at six different timepoints in three different iPSC cell lines ( Figure 1A). Differential gene expression analysis deployed on all samples with timepoints as independent variables revealed 22 380 differentially expressed genes (DEGs). The vast majority of DEGs were protein coding genes, lncRNAs, and processed pseudogenes, that showed pronounced intra-timepoint clustering as indicated by principal component analysis (PCA) and hierarchical clustering (Figure 1B,C; for analysis of individual cell lines see Figure S1A, Supporting Information). Enrichment analysis highlighted significantly enriched neuron-related pathways and entities indicating the transition from pluripotent cells to iDN morphology and function such as growth cone, axon guidance, dendritic spine development, synapse assembly, and glutamatergic synapse (Table  S1, Supporting Information).
Using weighted gene correlation network analysis (WGCNA) on a consensus inner joint dataset of expressed genes from all three cell lines (20 887 genes), we found 49 consensus gene modules with an average preservation of 0.89 and highly shared preservation between cell lines as indicated by the temporal trajectories ( Figure 1D,E, Table 1, Figures S1B,C, S5-S7, Supporting Information). Fourteen modules showed significant enrichment for neuron-related pathways associated with neurite outgrowth, synapse development, and ion channels, of which two modules (yellow and red) were highly enriched for neuronal related pathways such as synapse, neurite outgrowth, and development and were therefore involved in all three aspects ( Figure 1E,F; Table S2, Supporting Information). Networks of the top 30 hub genes of the respective consensus modules are presented based on Pearson correlations of gene expression data ( Figure 1F; Table S3, Supporting Information).
Sensory neuron fate determination and neural crest cell development were initiated by two networks of highly connected genes (lightyellow module, magenta module), containing NEUROG1, NEUROD2, NEUROD4, PRDM12, and SSTR2 as significant hub-genes ( Figure 1E,F). Analysis of important pluripotency, neural crest, and nociceptor markers revealed a general timedependent decrease of pluripotency and increase of neuronal and nociceptor markers. The neurogenic differentiation markers NEUROG1 and NEUROG2 followed highly correlated expression trajectories (Pearson R = 0.84, p-value < 0.05) during differentiation with peak expression at day 9, although only NEU-ROG1 was abundantly expressed (≈200 TPM versus ≈4 TPM at D9, respectively, Figure 1H). This finding supports successful reprogramming into small-diameter sensory neurons. [35,36] NTRK1 and RUNX1, two master regulators of nociceptor fate determination, showed highly correlated transient expression up to day 16, followed by an inverse-correlative expression pattern. Together with the abundant expression of the heat sensitive transducer ion channel TRPV1 and the voltage-gated sodium channel gene SCN9A, this suggested further specialization into nociceptors during maturation ( Figure 1F,H). In addition to RNA sequencing, we validated protein expression of SCN9A, panTRK (NTRK1, NTRK2, NTRK3), peripherin (PER), Tuj1 (TUBB3), and PIEZO2 at three different timepoints (D16-D36) during iDN differentiation, and protein expression fluorescence intensity trajectories resembled the RNA sequencing trajectories ( Figure  S2A,B, Supporting Information). In accordance with the temporal RNA trajectories, the vast majority (≈90%) of differentiated cells were positive for panTRK at D16 and positive for SCN9A (>80%) at D36, while negative for NANOG, Ki-67, TMEM119, and GFAP (Figures S2C and S3A-C, Supporting Information).
The magenta consensus module likely represents the network involved in neural tube and neural crest cell development and several genes known for their contribution in neural crest/neural tube development, such as PAX3, SOX9, TAF9B, were identified as hub genes within this module ( Figure 1F). Interestingly, EFNB3 coding for the trans-synaptic organizing protein ephrin B3 was identified as a magenta module hub gene, while the roles of other hub genes within this module, such as SPA17 or TKTL1 (transketolase like 1), are largely unattended in neurons to date (Table S3, Supporting Information). RUN and FYVE domain containing 3 (RUFY3) and AP-2 associated kinase (AAK1) as hub genes for the red module are well known neurodevelopmental proteins and, most interestingly, RAB7A, prolyl 4-hydroxylase (P4HB), and CDC42 as identified hub genes for protein-protein interactions are implicated in synapse development, neuronal growth, and survival. [37][38][39][40][41][42] Tropomodulin-1 (TMOD1), a tropomyosin-binding protein involved in the structuring of actin filaments in sensory neurons, and the neuronal survival factor and neuroprotectant Tomoregulin-2 (TMEFF2) were identified as hub genes for the yellow module. [42][43][44] Two more modules (saddlebrown and grey60) were implicated in both synapse and ion channel development. Importantly, several sodium and potassium channels (e.g., SCN2A, SCN7A, SCN10A, SCN11A, and KCNQ2) emerged as significant hub genes and may in addition to their functional roles qualify as critical indicators of neuron as well as synapse development ( Figure 1F,H, Tables S3,S4, Supporting Information for top hub genes).
Since several ion channel mRNAs were expressed, functionality of iPSC-derived sensory neurons was assessed using single cell whole-cell patch-clamp recordings ( Figure S4, Supporting Information). All cells exhibited large peak inward currents indicating activation of voltage-gated sodium channels and sustained outward currents through voltage-gated potassium channels ( Figure S4A, Supporting Information). Amplitudes of peak inward currents were significantly larger in AD2 cells compared to AD3 cells and this well reflects higher SCN1A mRNA log-normalized counts for each sample revealed high inter-timepoint specific separation and low intra-timepoint specific variation. C) Overall RNA expression similarity was determined using hierarchical clustering using correlation as similarity metric. D) WGCNA coexpression analysis dendrogram with module affiliations. Soft-threshold power [16] was determined assuming scale free topology. Topological overlap matrix was constructed and preservation between the three cell lines calculated and visualized in preservation as well as bar plots. E) Preservation plots between the three cell lines and enrichment analysis revealed biologically conserved modules. Enrichment analysis was performed for each module applying a hierarchical best per-parent algorithm on the g:Profiler (FDR-corrected p-value < 0.05) results to obtain most specific pathways enriched, which revealed neuronal enriched modules of four entities (synapse development, ion-channel development, neurite outgrowth, and development). F) Top 30 hub gene networks for multifunctional neuronal modules, ion channel and synapse development related modules, and neural tube development related modules. Hub-gene networks were determined using the module affiliation indicated by the k-module Eigenvector and the corresponding p-value. Hub gene networks between top 30 hub genes were constructed calculating the pairwise Pearson correlations and selectively choosing a hard threshold of r > 0.7 for interconnectivity. G) Top 30 hub gene network for the lightyellow module implicated in neural crest cell development. H) Literature-based pluripotency marker, neural crest marker, and nociceptor marker time trajectories were drawn from z-scored vst-normalized counts that were averaged across the three cell lines.  Figure S4A,E, Supporting Information). The resting membrane potential of ≈−55 mV for AD3 cells and ≈−47 mV for AD2 cells at D36 and their capability of firing action potentials indicated a mature state of iDNs but also cell line variabilities ( Figure S4B-D, Supporting Information). Common peptidergic neuron marker expression trajectories such as TAC1 (Substance P, SP) and CALCB (Calcitonin-gene related peptide, GRP2) were assessed with immune fluorescence microscopy and this revealed high expression of CALCB at D16 and TAC1 at D36 (Figures S1D,E, S2, Supporting Information).
To further support our model, we queried important markers of nociceptor development, recently determined in mice as well as transcription factors found enriched in human DRG. [26,45] Most of the development related markers were expressed from days 5-9 (HES6, TCF4, EZH, TFAP2, and FOXP4), indicating their importance in human iDN differentiation ( Figure S1G, Supporting Information developmental marker). Other murine nociceptor markers, such as SYT13, GAL, TH, and PRDM2, peaked around days 9-16 ( Figure S1F, Supporting Information, nociceptor marker). All human DRG transcription factors such as DRGX, PIRT, POU4F1, RBFOX3, SKOR2, and TLX23 were expressed in iDNs and steadily increased throughout iDN differentiation with a mean TPM above 5, supporting the sensory neuron phenotype and their role in human sensory neuron differentiation ( Figure  S1H, Supporting Information). We further compared our iDN signatures to human native DRG and cultured human DRG data sets and found that iDNs in general were more similar to cultured human DRGs compared to native DRG. iDNs became transcriptomically more similar to human DRG during differentiation from D16 on, however, since DRG contain a mix of different cell types including sensory neurons, Schwann cells, macrophages, and cells of the vasculature, a full overlap of gene sets cannot be expected, [46,47] Figures S1J,K, Supporting Information).
Given that many of the key genes are well known checkpoints in murine sensory neuron differentiation and human sensory neurons, [26] our analysis indicates a largely conserved sensory neuron developmental program from mouse to human.

Transcriptomic Trajectories Identify Five Distinct Stages of Human Nociceptor Differentiation
To further investigate iDN developmental stages, modules with similar gene expression time trajectories were merged by agglomerative hierarchical clustering, and singular value decomposition (SVD) was performed. Individual genes were projected onto the first two SVD-components, which resulted in a transcriptomic clock indicating different developmental stages (Figure 2A). Cluster analysis revealed six gene expression super-clusters with unique time trajectories that resembled pluripotency (2255 genes, mainly associated with D0), early differentiation (6057 genes, D0-D9), early neural progenitor (2016 genes, D5-D9), late neural progenitor (564 genes, D16-D26), and finally the nociceptor stage (8027 genes, transient expression until D36), as well as one super-cluster showing peak expression both during pluripotency and maturation (1085 genes, D0, and D16-D36) ( Figure 2B-D). Enrichment analysis indicative of the functional relevance for each respective stage showed that neuron-related pathways were induced as early as in the Figure 2. mRNA trajectory-based stages of nociceptor differentiation. A) Singular value decomposition (composition of three matrices) performed along the gene expression axis using averaged z-scored vst-normalized counts and projected onto the first two components revealed a transcriptomic clock indicative of different nociceptor stages. Dot color represents the clusters affiliation identified by merging WGCNA module with similar expression profiles via hierarchical clustering using the "average" agglomerative algorithm and mean z-scored variance-stabilized counts. B) Pie-chart of the distribution of mRNAs for each stage of nociceptor differentiation with numbers of genes. Outer circle represents the main stages during nociceptor development and inner circle shows the numbers of WGCNA modules affiliated with each main stage as well as the overall number of genes associated with each main stage and module. C,D) Mean z-scored vst-normalized curves for each module affiliated with the representative main-stage. E) Top 3 g:Profiler enrichments for each main stage, ranked based on the negative log10 p-value for 3 annotation spaces (GO:BP, GO:CC, GO:MF). The best-per parent approach was used to identify specific enriched pathways as described in Zeidler et al. (2020). www.advancedsciencenews.com www.advancedscience.com early neural progenitor phase (e.g., neuron differentiation). Later stages included neurite outgrowth related pathways (e.g., neuron projection) and synapse-related pathways (e.g., glutamatergic synapse, synaptic membrane; Figure 2E, Table S5, Supporting Information) and several neuronal related diseases, such as cognitive disorders, neurodegenerative disease, and nervous system disease (Table S6, Supporting Information). We further identified hub genes for each differentiation stage based on functional connectivity using StringDB network analysis. This revealed MYC and also RHOA as hub regulators in the pluripotency phase, and AKT1 and SRC encoding serine/threonine kinases as well as UBC in the nociceptor stage (Table S7, Supporting Information). We found that SOX9 (magenta module) and MED13 (grey60 module), which are strongly involved in neural tube/neural crest cell development, were most abundantly expressed between days 5-9, ( Figure 1F, Tables S3,S7, Supporting Information).

microRNAs Show Differentiation Stage-Specific Expression
As numerous protein-coding genes are regulated by miRNAs, we performed small RNA co-sequencing from the same samples to identify timepoint specific expression trajectories of differentially expressed miRNAs (DEmiRs) during nociceptor differentiation. PCA and hierarchical clustering deployed on all samples revealed 979 DEmiRs that similarly to DEGs showed pronounced intratimepoint clustering although with slightly higher variances ( Figure 3A,B; for analysis of individual cell lines see Figure S8A-C, Supporting Information). WGCNA analysis of a consensus inner joint miRNA dataset from all three cell lines (1647 miRNAs) revealed 17 highly correlated miRNA consensus modules as indicated by temporal trajectories with an average preservation of 0.88 ( Figures S8E, S9, Supporting Information). SVD resulted in five main differentiation stages based on miRNA expression trajectories that overall resembled the gene expression trajectories ( Figure 3C-E). Interestingly, no miRNA trajectory cluster was found that resembled the pluripotency/maturation stage. The top regulated miRNAs per super-cluster were hsa-miR-302a-5p for the pluripotency stage, hsa-miR-25-3p for the early differentiation stage, hsa-miR-363-3p and hsa-miR-96-5p for early and late neural progenitor stage, and hsa-mir-183-5p for the nociceptor stage ( Figure 3F). miRNA qPCR validation of hsa-miR-302a-5p, hsa-miR-25-3p, hsa-miR-103a-3p, hsa-miR-355-5p and hsa-miR-146-5p revealed similar time-dependent trajectories as found mRNA. miRNAs were chosen based on temporal trajectories as well as expression strength to cover low expressed (hsa-miR-103a-3p) but also highly expressed miRNAs (hsa-miR-302a-5p, Figure S8F, Supporting Information).
The majority of the consensus miRNA modules was affiliated with the pluripotency stage (7 modules) and the nociceptor phase (4 modules), indicating that the miRNAs within these modules may function as critically important master switches controlling large gene sets essential for these stages ( Figure 3G).
To investigate which pathways were regulated by miRNAs during iDN development, a miRNA::mRNA target edge database was constructed that included a combined target score based on the DIANA microT-CDS prediction score, the correlation between miRNA and mRNA expression (vst counts) in our samples, the percentile-ranked miRNA base expression as well as the target validation status (using StarBase, miRTarBase, and Tarbase). All targeted genes (correlation < −0.7, target-score > 70) for each module were determined and reverse enrichment analysis was performed. Predominantly miRNA modules belonging to the pluripotency and early differentiation phases were enriched for neuron related pathways such as growth-cone, axon guidance (turquoise module, pluripotency phase), synapse (purple module, pluripotency phase), and glutamatergic synapse (red module, early differentiation phase). In particular, the miR-17≈92 cluster was highly expressed and associated with the early differentiation of iPSC-derived nociceptors (D5-D9). In addition, hsa-miR-20a-5p, hsa-miR-92a, hsa-miR-19b, and hsa-miR-18a were highly regulated and associated with glutamatergic synapse formation, neuronal cell body, and endomembrane system organization pathway enrichments ( Figure 3F,G).

miRNAs Regulating Hub Genes
Since hub genes (k-module eigenvector (kME) > 0.8) were the main regulators per module, only those were considered for miRNA::mRNA target analysis. An adjacency matrix was calculated and UMAP dimensional reduction on a sparse target-score matrix followed by k-Means cluster analysis (n clusters = 8) performed to identify miRNA::mRNA interaction clusters associated with nociceptor differentiation stage ( Figure 4A). Cluster analysis revealed four miRNA::mRNA interaction clusters associated with the nociceptor stage ( Figure 4A), which were significantly enriched for neuron-related pathways such as axon guidance (cluster 7), synapse (cluster 1), presynapse (cluster 2), and glutamatergic synapse (cluster 6) (Table S8, Supporting Information).
Of the consensus gene modules, especially the red, yellow, turquoise, violet, and purple modules were highly targeted by miRNAs belonging to the miR-302 family, the miR-17≈92 cluster, the miR-106≈25 cluster, hsa-miR-31, and the miR-200 family ( Figure 4B). The anti-correlative miRNA::mRNA expression trajectories revealed two distinct patterns: clusters 1 and 7 showed a transient increase in gene expression and decrease in miRNA expression, whereas clusters 2 and 6 showed no change until D9, followed by a sustained increase in gene expression and decrease in miRNA expression ( Figure 4C). This documents the critical importance of miRNA inhibitory control over certain gene sets, which is released by decreasing miRNA expession during iDN development.
Overall, 586 miRNAs were highly expressed during pluripotency and expression decreased with proceeding differentiation. To further dissect the role of specific miRNAs in the development of iDNs, we performed functional enrichment analysis of the hub-genes of the consensus gene modules in the same interaction cluster, which are therefore targeted by similar miRNAs. This revealed a strong post-transcriptional control of synapserelated pathways (mainly in the red and yellow consensus gene modules, Figure 4D, Table S9, Supporting Information). Within the red module, several hub genes targeted by the miR-302 family were significantly enriched for synaptic vesicle related pathways (VAMP2, SCD5, SYBU, SYNJ1, cluster 1). Hub genes that were predominantly targeted by hsa-miR-367-3p, hsa-miR-25-3p, and hsa-miR-31-5p were enriched for presynapse/ion-channel complexes (SCN9A, KCNK3, GAP43, CACNA1B, cluster 2) ) is annotated at the corresponding axis and inter-timepoints as well as inter-cell line variance investigated. B) Hierarchical clustering using the "average" agglomerative algorithm based on similar miRNA expression revealed high intra-timepoint clustering with low inter cell-type variances. C) Singular value decomposition (SVD) of each individual miRNA based on expression projected onto the first two SVD components. Each individual miRNA is depicted by a dot and colored by the cluster affiliation determined applying the "average" hierarchical clustering algorithm on the main averaged temporal WGCNA module trajectories, which revealed a transcriptomic clock of miRNA expression with 5 main stages modelling the differentiation of iPSCderived nociceptors. D,E) Resulting clusters were ranked based on the transcriptomic clock and average (cell line) z-score vst-normalized expression patterns representing time trajectories of each individual WGCNA module belonging to the main stage and are visualized with line plots. F) miRNA abundance analysis for each individual differentiation stage based on baseMean expression (normalized DEseq2 counts) showing the Top 5 most abundant miRNAs for each main stage, and the temporal expression patterns depicted in a heatmap (z-score standardized vst(counts)). G) Reverse miRNA enrichment profiling for each WGCNA module ranked according to the main differentiation stages was determined using the so constructed miRNA_edge.csv database target-score. Only miRNA::gene interactions with a target-score above 70 (most of them highly predicted and validated) and an inverse correlation below r < −0.7 were considered, and top regulated pathways of targeted genes were determined using g:Profiler, ranked by the negative log10 p-value and visualized in the heatmaps.
Several miRNAs such as hsa-miR-25-3p targeted genes within several clusters and modules and could therefore serve as master switches suppressing genes associated with mature neuron machinery and function ( Figure 4D,F). Enrichment analysis of all high confidence targets revealed an important function in nervous system and postsynapse development (Table S10, Supporting Information). Furthermore, the synaptic genes SYT1 and SNAP91, which are associated with synapse development but also implicated in central nervous system disorders such as autism spectrum disorder or epilepsy, emerged as highconfidence targets of hsa-miR-25-3p ( Figure 4G) highlighting this miRNA as an interesting novel target for neuron-related disorders.

Neurite Outgrowth Modulated by miRNAs
The formation and outgrowth of fine processes such as dendrites and axons is a prerequisite to connect sensory neurons to their target tissue and to postsynaptic neuronal assemblies, [48] but also to allow for nerve regeneration following peripheral nerve injury. [49] Of the seven consensus gene modules related to neurite outgrowth ( Figure 1E), three (royalblue, violet, and floralwhite) were strongly enriched for neurite outgrowth, cell migration, and cytoskeleton regulation (Figure 5A,B, Table S2, Supporting Information). All three modules showed a strong increase in gene expression after D9 and contained well known hub genes that are mechanistically involved in neuronal outgrowth (royalblue: COL1A1, MCAM, RHOC, RHOB, RHOJ, DUSP6; violet: TFBR1, FAM114A1, ROCK2, POSTN, floralwhite: AMIGO1, CLSTN1, PHACTR1, MAPK1, PTPN1; Figure 5A,B, Table S4, Supporting Information).
Based on their expression trajectories and target prediction scores emerging from our pipeline, we selected several interesting miRNAs that could be critically involved in these processes as regulators of larger gene sets spanning multiple clusters and modules. Hub genes such as RHOC, MBNL2, TACC1 from the royalblue and ROCK2, DST, and RAB's from the violet module were predominantly targeted by the human miR-302 family (cluster 1, Figure 5C,D). The resulting regulation of RHO GTPase activity and integrin binding putatively suppressed neurite outgrowth and cell migration until D5. Since the miR-302 family shows a diverse target-spectrum, an analysis of all highconfidence targets (target-score > 0.75, r < −0.7) was performed to pinpoint functionality even further. This revealed a significant enrichment for neurite outgrowth related pathways microtubule binding, axonal transport, growth-cone, and dendrite (Table S11, Supporting Information).
Interestingly, the hsa-miR-103a-3p target space within cluster 7 appeared to be implicated in the suppression of neurite outgrowth and axonal elongation by targeting genes belonging to the floralwhite dendrite module but also to clusters enriched for phosphatidylinositol-3-phosphate binding (violet) and cell morphogenesis (royalblue) ( Figure 5C-E, Table S12, Supporting Information).
For a third relevant miRNA, hsa-miR-31-5p, high-confidence targets were found in all 3 modules in most interaction clusters, suggesting a general control via this miRNA and therefore an important role in neurite outgrowth and axonal development (Figure 5E, Table S13, Supporting Information). When performing pathway enrichment for all anti-correlated high-confidence targets, it became evident that hsa-miR-31-5p might further exert its function by suppressing expression of RAB GTPases RAB14, RAB1A, RAB21, RAB2B ( Figure 5D,E), which regulate membrane trafficking. [50] We further identified hsa-miR-155-5p as significantly enriched for neurite outgrowth/synapse related pathway predominantly targeting hub-genes of cluster 7, which is crucial for the suppression of axonal elongation in mice [51] (Figure 5D,E, Table S14, Supporting Information).

Ion Channel Expression as a Signature of Functional iDN Maturation
Ion channels are essential for the general functioning of cells and in particular of excitable mature neurons including nociceptors. Therefore, we investigated the expression trajectories of . miRNAs regulate pre-and post-synapse development via module hub genes. A) A miRNA::hub-gene target-space interaction map (adjacency matrix using the target-score between each miRNA::gene interaction) was constructed by using the miRNA_edge.csv database with a target-score above 70 and an inverse correlation below −0.7 between the gene target and the miRNA. UMAP dimensional reduction was performed to reduce the resulting adjacency matrix along the targeted genes axis, which revealed higher local connectivity of genes targeted by the same miRNAs as well as high local connectivity of genes representing the same main stage. Each individual dot following UMAP construction was colored according to the affiliation to the differentiation main stage, which revealed separation between early and late iDN differentiation based on miRNA targeting. k-Means clustering (with n_cluster = 8) of the adjacency matrix revealed 8 clusters of targeted genes targeted by different miRNA groups. K-Means affiliation is projected onto the UMAP and k-mean clusters affiliated with late-stage nociceptor development were chosen further analysis. B) k-Means clustering revealed 4 clusters corresponding to the late stages of iDN differentiation, WGCNA module affiliation for each gene was determined and the number of targets within each WGCNA module for each k-Means cluster was determined as well as the number of genes targeted by each individual miRNA. Results were ranked descending to evaluate miRNA as well as module importance for each k-Means cluster. C) Four clusters showed neuronal enriched gene pathways in the g:Profiler analysis; miRNA (lightred) and gene expression (skyblue) trajectories for each of the k-Means clusters were represented as z-scored vst-normalized temporal expression of each gene and miRNA trajectory belonging to the k-Means cluster; an averaged gene (red) and miRNA (blue) expression curve was calculated and Pearson correlation calculated between the averaged curves. D) Top 30 miRNA targeted hub gene networks for each WGCNA module per k-Means cluster were constructed using the number of miRNA interactors as node-size as indicator of suppression strength within the module. E,F) miRNAs with the highest number of miRNA targets using groupby and counting functions provided by python for the red and yellow modules were determined for each k-Means cluster and z-scored data visualized for each of the 4 clusters. G) Protein-protein interaction network (StringDB) of all hsa-miR-25-3p targets and the main differentiation stage affiliation indicated by the arc color visualized using a chord interaction plot. all currently assigned ion channels ( Figure 6A). Six consensus gene modules were significantly associated with ion channel activity, including the brown and cyan modules ( Figure 6B; brown: top hub genes: ESRP1, KCNK6; enrichments: gated channel activity, integral component of the plasma membrane; cyan: top hub genes: AL365203.2, CYB561D1, SLCO2B1, NGFR, BIN1; enrichments: plasma membrane, ion channel regulator activity, Table S2, Supporting Information). NGFR and BIN1 are well established regulators of ion channel trafficking to the membrane and both were highly correlated with the expression of the DRG specific voltage-gated sodium channels SCN10A and SCN11A (r > 0.86, p < 0.05). [52] A large number of ion channels was already expressed at the iPSC stage and early differentiation stages, and underwent a drastic reduction of expression during differentiation (Figure 6A). This included potassium channels KCNQ1, KCNQ5, KCNS1, and KCNK6, voltage-gated sodium channels SCN4A, SCN5A as well as HCN4, CLCN2, CACNA1G, and TRPM7. As expected, several ion channels were predominantly expressed in the nociceptor stage (e.g., SCN1A, SCN2A, SCN3A, SCN9A; KCNQ2, KCNJ6, KCNJ13, KCNJ16; CLCN3, CLCN4; HCN1, HCN2) ( Figure 6A). The two hub genes SCN10A (grey60 module) and SCN11A (cyan module) peaked around D16, and SCN8A a voltage-gated sodium channel which is essential for neuron Figure 6. miRNAs regulate ion channels via hub genes. A) Heatmap of ion channel expression during iDN development based on vst-normalized counts averaged across the cell lines, ranked based on main differentiation stage affiliation and iDN developmental progress was calculated. B) Top 30 hub genes of the brown module highly enriched for ion channel related gene pathways and expressed in iPSC-stage and top 30 hub genes of the cyan module with ion channels marked in red as top hub regulators. C) Time trajectories of ion channels (blue) expressed in the early and late stages of iDN differentiation and the mean projection curve, as well as miRNA expression trajectories (rose) found to regulate ion channels (right, top). D) miRNA::mRNA network analysis considering the target-score as edge weight, miRNAs (blue) and genes (red) with more than 3 edges are annotated and node size represents the Eigenvector centrality. E) The number of genes targeted by a miRNA (red barchart) and the number of targeted genes by miRNAs (blue barchart) represented in ascending orders. www.advancedsciencenews.com www.advancedscience.com function in mammalian neuronal tissues and in the pathogenesis of neuropathic pain [53][54][55] peaked already at D9 at the transition from pluripotency to neuronal differentiation ( Figure 1F, Figure 6B). This suggests a critical role of SCN8A, SCN10A, and SCN11A not only for nociceptor function but also differentiation and maturation. In general, these findings support the idea that neuronal excitability and action potential discharge may not only be required for neuronal function, but also critically affect cell fate decisions and establishment of target tissue innervation and neuronal circuitry. [56]

miRNAs Regulating Ion Channel Expression
Analysis of miRNAs targeting ion channels revealed two highly connected communities of miRNA::ion channel networks that could be separated based on their differentiation timepoint affiliation. Ion channels and corresponding miRNAs revealed highly significant inversely correlated trajectories early (D0-D5, Figure 6C) and later during differentiation (D16-D36; Figure 6C), as indicated by the mean trajectories drawn for miRNA and mRNA vst counts.
The gap junction related ion channels GJC1 and GJA1 were among the most targeted ion channels in the iPSC stage (Figure 6D,E). Also, the transient receptor potential-melastatin-like 7 channel (TRPM7) was expressed mostly in the iPSC stage, and appeared to be post-transcriptionally controlled by at least 17 miRNAs during differentiation ( Figure 6D,E). At the more mature stages, several ion channels emerged, whose relevance for the transduction of painful stimuli, excitability, or action potential propagation is well accepted. Two important nociceptor ion channels SCN2A and SCN9A [57,58] were highly and inversely correlated with hsa-miR-20b-5p, hsa-miR-15b-5, hsa-miR-18-5p, hsa-miR-423-3p, hsa-miR-93-5p, and others, which decreased during the differentiation, while SCN2A and SCN9A expression increased. Also, the chloride transporters and channels CLCN4-6 were strongly targeted by miRNAs, with CLCN6 being controlled by 14 miRNAs and CLCN5 being targeted by 23 miRNAs. Interestingly, all three chloride channels were targeted by members of the hsa-miR-302 family. In general, the miR-302 family and hsa-miR-372 were found to suppress more than ten ion channels ( Figure 6E).

NOCICEPTRA Resource
For general use and utilization of our data set, we provide the online tool NOCICEPTRA based on the Python "streamlit" framework, to analyze and visualize time trajectories of genes and miR-NAs of interests. Both DESeq2 variance-stabilized counts as well as TPMs can be queried for each gene and miRNA present in the dataset, and a correlation matrix for the queried genes calculated. We further implemented an integrated exploratory data analysis tool, to browse all consensus gene and miRNA modules, enrichments, and the top 30 hub gene regulators, as well as miRNA interactions with the hub regulators. In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) (www.genome.jp/kegg) and disease pathways (www.disgenet.org) can be browsed, and standardized Pearson residuals calculated to identify potential windows of susceptibility during iDN differentiation, together with an integrated g:Profiler enrichment analysis (www.biit.cs.ut.ee/ gprofiler) and StringDB PPI analysis (http://www.string-db.org/, Figure 7A). Also, single miRNA target spaces can be queried with adjustable parameters (e.g., target-score, correlation) to determine the putative function of miRNA during iDN differentiation indicated by StringDB PPIs and enriched pathways ( Figure 7A) and single genes or gene sets miRNA target spectra determined.
Since a number of severe human pain disorders are strongly associated with gain or loss-of-function mutations of ion channels in neurons, [59][60][61][62][63] we applied the NOCICEPTRA tool to exemplarily investigate possible susceptibility windows during iPSC differentiation using the custom-gene set function for the two disease entities "neuropathic pain" and "migraine with Aura" evaluating the risk association score (association score > = 0.2) obtained from the OpenTarget.org database. The overall number of genes within each super-cluster was used as background to identify disease linked genes, and standardized Pearson residuals were calculated to determine the contribution of each cell to the contingency table statistics ( Figure 7B, Table S15, Supporting Information). Supercluster enrichments for neuropathic pain and migraine with aura associated genes revealed a significant enrichment during the nociceptor stage for migraine with Aura, while for neuropathic pain the largest number of genes were associated with the nociceptor stage as expected, such as ORPM1, CHRNA4, CACNA1B and the metabotropic glutamate receptors GRM5, however no significant enrichment was detected (Figure 7B, contingency table statistics > 0.05; Table S19/S20, Supporting Information).
Interestingly, several genes were abundantly expressed already during iDN differentiation and associated with migraine with aura such as HTT, TOR1A but also DRD2 which suggests a putative implication in the development of neural progenitor cells ( Figure 7B). Additionally, ion channels associated with migraine with aura but also neuropathic pain were abundantly expressed in the nociceptor stage such as SCN1A, SCN2A, SCN3A, SCN9A, KCNA1, HCN2, CLCN3, and KCNA2 but also serotonergic receptors such as HTRA2 and glutamatergic receptors such as GRM5 and GRIN3A, proving suitability of the model to study both diseases ( Figure 7B, Table S20, Supporting Information).
Mapping miRNAs associated with neuropathic pain, chronic pain, or migraine to diseases using the human miRNA Disease Database (HMDD v.3.2, www.cuilab.cn) revealed two groups of miRNA which were associated with disease entities ( Figure 7B, righ panel). These were significantly expressed in early differentiation (D0-D5) and in nociceptor stage (D16-36). microRNAs associated with early differentiation putatively imply a high regulatory control of iDN in the differentiation and maturation. Interestingly, hsa-miR-155-5p was significantly associated with both disease entities "neuropathic pain" and "migraine" and further was associated with neurite outgrowth related genes in this study. We identified common targets associated with migraine and the hsa-miR-155-5p high confidence target space such as PLP1, TRAK1, SYNJ1, and CREB5 ( Figure 7B).
Based on these findings, it can be anticipated that impairments of migraine associated genes and associated miRNA regulation can be assessed experimentally already during early differentiation/neural progenitor phases of iPSC derived nociceptor models. The same appears to apply for genes associated with neuropathic pain questioning the role of those genes not only in mature Figure 7. The NOCICEPTRA tool and disease enrichment analysis. A) Schematic diagram of all databases incorporated in the NOCICEPTRA tool as well as analysis outputs. B) Susceptibility windows of diseases determined using the NOCICEPTRA tool with OpenTarget.com risk associations of migraine with aura and neuropathic pain. Contingency table analysis as well as Pearson residuals were determined for "neuropathic pain," "migraine with aura" and a chord plot with all genes implicated in the queried disease with PPI confidence scores (StringDB) as edge weights and the main differentiation stage affiliation as arc color was constructed. Standardized Pearson residuals were used to determine potentially enriched/depleted disease-associated genes for each differentiation stage (Left Panel). Human miRNA disease database (HmDD) was used to detect miRNAs associated with migraine disorders, chronic pain, and neuropathic pain, and temporal trajectories of associated miRNAs were determined and clustered using hierarchical clustering of z-scored vst(counts) and visualized as a heatmap along the time.
www.advancedsciencenews.com www.advancedscience.com nociceptor function but also nociceptor development and potential genetic impairments affecting both.
Hence, NOCICEPTRA was queried and compared against CORTECON expression patterns of important developmental marker genes during cortical neurogenesis which revealed distinct similarities and expected differences. As a first example, FGF receptors (FGFR1, FGFR2, and FGFR3) contribute significantly to corticogenesis and peak around the early corticogenesis in the CORTECON dataset, [23] However, expression was only conserved for FGFR2 and FGFR3 in iDNs peaking around D5 which expectedly indicates neural crest cell generation common to both neuron types. By contrast, FGFR1 trajectories steadily decreased, suggesting that FGFR1 depression is necessary for nociceptor genesis only ( Figure S10A, Supporting Information). Second, WNT signaling is crucial for corticogenesis/cortex formation, and WNT7B and WNT8A expression increases during cortical specification. [23] This was not the case for iDNs where expression of WNT7B and WNT7A genes was low at all stages. Querying all WNT genes to elucidate potential susceptibility windows revealed that the majority of WNT genes (WNT1, WNT4, WNT10B, WNT3A, and WNT5B) appeared to contribute to neural tube formation and neural crest cell development around D5 and D9 and that especially WNT1 was abundantly expressed (Figure S10B,C, Supporting Information). [64,65] Both complex NOCI-CEPTRA and CORTECON data sets and tools can therefore be exploited for mechanistic studies dissecting differential developmental regulation of neuron subtypes in the peripheral and central nervous system.

Discussion
Fueled by the urgent need for better and accessible human-based model systems, increasing efforts have been made to understand principal transcriptomic signatures responsible for the development of various tissue types, including cortical neurons by developing iPSC-derived differentiation protocols. [23,66] The majority of studies to date focus primarily on expression profiles of protein coding genes thereby disregarding the critical importance and application potential of miRNAs for the regulation of cell development and differentiation. In order to fill this gap, we performed the first unbiased temporal transcriptomic analysis of paired mRNA::miRNA expression during iPSC-derived nociceptor (iDN) differentiation in three different iPSC cell lines and detected several entities where miRNAs controlling an entire set of target genes apparently serve critical roles as regulatory master switches.

mRNA and miRNA Signatures Characterize Five Distinct Stages of Nociceptor Differentiation
Different iPSC donor lines and clones can react to differentiation protocols in a highly variable manner thereby impeding the identification and validation of conserved developmental pathways and molecular disease signatures. [67,16,24] All three iPSC lines in our study developed into a nociceptor like phenotype and were functionally capable of firing action potentials with largely conserved transcriptomic profiles and highly correlative inter-module expression patterns. Expression trajectories with an increased expression of RUNX1, a reduction in TrkA expression but an increased expression of TAC1 were consistent with previous literature, and the applied protocol predominantly generated peptidergic sensory neurons. [3,19] Thereby, five independent stages of sensory neuron differentiation were characterized by representative marker genes such as POU5F1 (pluripotency stage), NEUROG1 (early differentiation stage), NTRK1 (earlydifferentiation-neural progenitor stage), SOX10 (neural progenitor stage), and CLCN3/SCN9A (nociceptor stage), as well as gene pathway enrichments for each stage. [68][69][70] Characteristic trajectories of mRNAs and miRNAs emerged for the three relevant features of maturating nociceptors: neurite growth, synapse, and ion channel expression.

Hub Genes and miRNAs Orchestrating Differentiation and Fate Determination
miRNA-driven posttranscriptional gene regulation is pivotal during neurogenesis and loss of the miRNA-synthesizing enzyme Dicer leads to impaired neural stem cell differentiation and abnormal brain development. [71,72] Several miRNAs and their target genes are well accepted regulators of neuronal development and differentiation at defined stages of nociceptor differentiation, such as members of the miR-302, the miR-125b, and the let-7 families. [73] In particular, the miR-17≈92 cluster, was highly expressed and associated with the early differentiation of iPSCderived nociceptors (D5-D9). The miRNAs hsa-miR-20a-5p, hsa-miR-92a, hsa-miR-19b, and hsa-miR-18a were highly regulated and associated with glutamatergic synapse formation, neuronal cell body, and endomembrane system organization pathway enrichments. The strong association of high confidence targets with neurite outgrowth and neurotransmitter systems, moves the suppression of miR-302 into focus as a critically important mechanism for neuronal development and fate decision. [74][75][76] However, also miRNAs sparsely investigated in neuron development, such as miR-106a/b, miR-199b, and miR-504, were abundantly expressed, and presumably involved in suppressing pluripotency related genes, thereby inducing and maintaining the neuronal phenotype.

Neurite Formation Related Pathways Controlled by Hub Genes and miRNAs
Neuronal differentiation is accompanied by significant morphological changes including development of filopodia and outgrowth of neurites/processes particularly in the neural progenitor stage (D9-D16). miRNAs appeared to be involved in certain aspects of neurite formation. The most abundantly expressed intragenic miRNAs hsa-miR-363-3p, hsa-miR-20b-5p, hsa-miR-106a-5p, which are all embedded in the lncRNA AC002407.2, [77] targeted a variety of pathways implicated in axon guidance and showed strong correlation with NEUROG1, supporting their importance for neuron morphology and neurite outgrowth. [78,79] Interestingly, miR-20b/miR-106 directly suppress NEUROG2, [80] thereby indirectly favoring NEUROG1 expression and nociceptor differentiation.

www.advancedsciencenews.com www.advancedscience.com
Half of all consensus modules associated with neuronal development were implicated in neurite development and neurite growth. Within these, two classes of genes were overrepresented: DUSPs-dual phosphatases, and RABsbrain-associated Ras small G-proteins, which are well associated with neurite outgrowth. [81][82][83] The antiproliferative miRNAs hsa-miR-31-5p but also hsa-miR-103a-3p were highly abundant early in differentiation. [84] Since they specifically target DUSP and RAB genes, these two miRNAs could act as main drivers balancing neurite growth in sensory neurons during nociceptor differentiation and maturation.
Neurite outgrowth is also regulated by Rho-related pathways, and RHOC, RHOB, and RHOJ were highly correlated with CDK8, for which regulation of actin-cytoskeleton assembly is well established. This pathway appears to be unique for sensory neuron differentiation in humans: whereas RHOA is highly expressed up to postnatal day 1 in rodents, RHOA in human sensory neurons gradually declines, while RHOC increased in iPSC-derived nociceptors, suggesting opposing biphasic regulatory control of cytoskeletal assembly during differentiation. [85,86]

miRNAs Regulating Ion Channels in iPSCs and Maturing iDNs
Our results further support the importance of the miR-17≈92 cluster, since it may act as a general ion channel regulator by putatively downregulating an entire group of important ion channels, including the voltage-gated sodium channels SCN2A and SCN9A, the pacemaker channels HCN1 and HCN2, the voltagegated potassium channels KCNA1, KCNB1 and KCNH5, the inward-rectifier KCNJ6 as well as TRPM4, most of which are dysregulated in neuropathic pain or implicated in neuronal development, thereby putatively protecting from excitotoxcity. [87,88] In particular, the target-space of miR-20a (a member of the miR-17≈92 cluster) included multiple genes specifically regulating neurite outgrowth, making it a critically important hub miRNA for nociceptor differentiation. However, we also detected new sets of miRNAs with a putative role in targeting ion channels. miR-182 and miR-26b targeted the gap-junction related channels GJC1 and GAJ1, while miR-20b-3p (and the miR-17≈92 cluster) and the miR-302 family were found to target SCN9A, TRPM7, or CLCN5-CLCN7 channels. GJC1 and GJA1 might be of particular importance for retaining pluripotency as they facilitate the reprogramming efficiency toward iPSCs. [89] Also, the expressed TRPM7 channel appeared to be strongly controlled by at least 17 miRNAs during differentiation and this stresses its role in neuronal differentiation in particular of mechanosensitive sensory neurons through, however, unknown mechanisms. [90,91] At the more mature stages, several ion channels emerged as targets of miRNAs, whose relevance for the transduction of painful stimuli, excitability, or action potential propagation is well accepted. Two important nociceptor ion channels, SCN2A and SCN9A, [57,58] were suppressed in immature stages by highly expressed miRNAs of the miR-17≈92 cluster which decreased during differentiation, while SCN2A and SCN9A expression increased. Also, the chloride transporters and channels CLCN4-6 were strongly targeted by miRNAs especially by the miR-302 family. In particular the miR-302 family suppressed ten ion channel targets. This moves the miR-302 family into focus as a critical and efficient hub regulator mechanism suppressing entire clusters of channels and other machineries that are silenced in pluripotent cells.

Differentiation of the Synaptic Release Machinery
Although native primary sensory afferents do not form autapses, [92] iDNs express synaptic proteins. [18] This is surprising but may be related to the unique feature of nociceptors to release neuropeptides and glutamate from differing vesicle pools at peripheral and spinal terminals. [93] Two paralogs of the calcium-dependent activator protein for secretion are priming factors for synaptic vesicles containing glutamate and large dense-core vesicles containing neuropeptide. Essential components of the synaptic release machineries such as syntaxin, synapsin, SNAP91 as well as different RAB proteins were subject to the control by hub miRNAs. While hsa-miR-302a-5p is essential for neuronal differentiation and neural tube closure the present study for the first time hints toward the implication of the miR-25 family in nociceptor development and suppression of synapse machinery related genes. [94][95][96][97] Yet, it is unclear how the exocytotic release machinery from these two vesicle types in nociceptive neurons is differentially regulated, wherein fast synaptic transmission and neuropeptide release are equally important. [98] To further investigate timepoint specific mRNA and miRNA regulation, we developed the open access NOCI-CEPTRA tool to extract expression patterns and time trajectories of miRNAs and their target gene space for mechanistic studies.

Querying Gene Trajectories Using the NOCICEPTRA Tool
Based on disease enrichment data, we found that a set of genes and miRNAs contributing to nociceptor development was also associated with chronic pain and migraine. To further extend the current knowledge on disease-relevant gene and miRNA trajectories during nociceptor differentiation, we developed the NO-CICEPTRA tool to visualize stage specific gene enrichments for diseases and KEGG pathways, as well as for single miRNAs. The NOCICEPTRA resource allows to explore pain and other sensory neuron-related disorders through the discovery of disease onset and susceptibility windows. The NOCICEPTRA platform is provided as a containerized online tool that enables the scientific community to retrieve the entire dataset of the experimentally determined expression trajectories, as well as to query custom gene sets of interest for pathway and disease enrichments. Querying our resource to differentiate expression patterns of important developmental marker genes identified during cortical neurogenesis (CORTECON) revealed distinct similarities and dissimilarities. As an example, the FGF receptor FGFR1 significantly contributed to corticogenesis but not nociceptor development, while FGRF2 and FRGR3 were found important for both nociceptor genesis and corticogenesis. This pinpoints the added value of CORTECON and NOCICEPTRA tools to generate novel mechanistic insight into neuron subtype programming.
Together, our findings suggest that miRNA::mRNA interactions act as critically important hubs for suppressing mature www.advancedsciencenews.com www.advancedscience.com protein coding gene sets during pluripotency and for silencing pluripotency genes when neurons serve their mature function in the nervous system. This posttranscriptional regulation via miR-NAs emerges as effective and of particular importance for phenotype decisions, neurite growth, and target finding as well as synaptic processes for nociceptor interaction with target tissues and second order neurons. In order to fathom the complexity of the modules involved in these unique functions we provide an open resource for the scientific community that can be used to query pathways and miRNA-controlled gene sets.

Conclusion
In summary, our paired long and short RNA co-sequencing approach allowed us the complex mapping of a developmental gene and miRNA expression atlas discovering major hub genes and miRNAs significantly contributing to the development of iDNs. Thereby, miRNAs and hub genes with specific roles in synapse development, neurite outgrowth as well as ion channel development were identified and the resource NOCICEPTRA generated allowing general exploration of gene and miRNA expression, hub modules and pathways as well as disease susceptibility windows.
Time-Course Paired RNA and Small RNA Sequencing: Paired long and short RNA sequencing were performed to identify temporal changes of gene expression during the development of nociceptor differentiation. Data were generated for six different time-points (i.e., D0, D5, D9, D16, D26, and D36) chosen based on morphological rearrangement such as major treatment changes during the protocol and neurite outgrowth in all three different iPSC cell lines (i.e., AD2, AD3 and 840; N = 18 per cell line).
To explore gene and miRNA expression patterns during the development of iPSC-derived nociceptors, cells were collected in triplicates for each cell line. To harvest cells, phosphate buffered saline (PBS) was used to wash off the remaining medium and to reduce the amount of cell debris. Cells were detached with 1× Tryp-LE Express (Sigma) for 3 min. Tryp-LE was inactivated with PBS in a 1:1 mixture and cells were centrifuged at 500 g for 5 min. Subsequently, supernatant was depleted and cells were snap frozen in liquid nitrogen and stored at −80°C until further processing. Long and short RNA co-sequencing were performed by IMGM Laboratories (Martinsried, Germany). To avoid batch effects, all samples were processed in parallel. Briefly, in order to isolate RNA and short RNA from frozen samples, the miRNeasy Mini Kit (QIAGEN) was used according to the manufacturer's instructions. RNA concentration and purity were assessed using NanoDrop ND-1000 spectral photometer (Peqlab) and the quality of the RNA was determined using RNA 6000 Nano LabChip Kits (Agilent Technologies) on a 2100 Bioanalyzer (Agilent Technologies). RNA quality based on RNA integrity (RIN) is presented in Table S17, Supporting Information. Furthermore, all libraries were quantified using the highly sensitive fluorescent dye-based Qubit ds DNA Assay Kit (Qiagen). Whole RNA library preparation was performed using the TruSeq stranded mRNA HT technology (Illumina) according to the manufacturer's protocol. All single mRNA libraries were pooled into one sequencing library pool with an equal amount of DNA per sample. After quantification, the final sequencing library was diluted to 2.25 nm followed by denaturation with NaOH. Short RNA sequencing was performed using the NEBNext small RNA Library Prep Kit for Illumina using 500 ng total RNA, according to the manufacturer's instructions. Purified short RNA sequencing libraries were quantified using the fluorescent dye-based Qubit ds DNA HS Assay Kit (Thermo Fisher Scientific) on a 2100 Bioanalyzer. After quantification, 2 nm of the sequencing library were used for further analysis. Both mRNA and short RNA libraries were clustered and sequenced using the Illumina NovaSeq600 and performed using single reads with a length of 75 bp. Technical quality parameters were evaluated using the SAV software and revealed that more than 85% of bases passed the Q30 bases threshold for all samples.
Long and Short RNA Read Processing and Differential Gene Expression Analysis: Derived mRNA FastQ files were aligned to the human reference genome (GRCh38.p13, release 33) provided by GENCODE (www.gencodegenes.org) using STAR (v2.7) software with the default parameters. [99] Short RNA sequencing libraries were prepared using the NEBNext Small RNA Library Prep Set for Illumina. To further process small RNA reads adapter sequences (AGATCGGAAGAGCACACGTCT-GAACTCCAGTCAC) were trimmed using Flexbar v3 [100] and aligned to the human reference genome (Parameters can be found in Table S18, Supporting Information). Raw counts for mRNA and short RNA were obtained sorted .bam files using HTSeq (v 0.11.1, only considering unique counts for RNA and short RNA sequencing). Raw mRNA and miRNA counts were determined by annotation to the GENCODE comprehensive gene annotation file (GRCh38.p13_chr_patch_hapl_scaf.annotations.gtf, www.genecod-egenes.org) and the hsa.gff file provided by miRbase (release 22.1, www.mirbase.org), respectively.
Since non-coding RNAs can be highly amigous reads, the CLC Genomics Workbench was further used to align ncRNAs against miRBase release v.21 and the Homo_sapiens. GRCh38.ncRNA database counting multimappers too. This data is available within the NOCICEPTRA tool (multimapper miRNA section) and as tables in the GitHub repository. To avoid overestimation using multimappers for miRNA expression and www.advancedsciencenews.com www.advancedscience.com correlation analysis only uniquely mapped miRNA reads were used in this publication.
Differential Gene Expression: Detection of DEGs and miRNAs (DEmiR) was performed using the Bioconductor R-package DESeq2 (release 3.10) and RStudio. [101] A first model was deployed merging all samples from each timepoint and each cell line, and evaluating differential gene and miRNA expression using the likelihood ratio test (LRT) with timepoint as independent parameter, to detect significant changes at any timepoint during the differentiation. A second model was deployed determining the number of DEGs and DEmiRs for each cell line at any given timepoint. Quantification of DEGs and DEmiRs, as well as log2-fold changes were estimated using the normal shrinkage function, and the FDR-corrected significance threshold was set to FDR < 0.05 for DEGs and FDR < 0.1 for DEmiRs. Normalized counts were further variance stabilized using the vst function.
To identify samples with an increased similarity of gene and miRNA expression, principal component analysis (PCA) as well as hierarchical clustering were applied on log-transformed counts. To determine the number of DEGs for each gene biotype, the Bioconductor R-package BioMart was used [102] to assign the biotype to each gene.
Weighted Gene Co-Expression Network Analysis: A weighted gene coexpression network analysis (WGCNA) was used to detect modules of genes that were highly correlated. [103] To be as unbiased and unsupervised as possible, normalized read counts obtained from the DESeq2 analysis for all mRNAs and miRNAs were filtered for genes and miRNAs that had at least 5 counts in 3 samples using Python's "Numpy" and "Pandas" package. A consensus dataset for genes and miRNAs was constructed separately using an inner joint on shared filtered genes and miRNAs for all three cell lines and quality as well as consistency of the dataset were quantified using the checkSets function as well as the goodSamplesGenesMS function integrated in the WGCNA package. Sample trees for each of the three cell lines were constructed using the hclust function. The blockwiseCon-sensusModules function (softthresholding power = 16 (mRNAs)/7 (miR-NAs), maxBlockSize = 37 000, network-Type = "signed," minModule-Size = 30, mergeCutHeight = 0.15, minKMEtoStay = 0.2 (mRNA)/0.0 (miRNA)) was used to construct modules of significantly correlated genes and miRNAs. Modules of gene expression were defined as a cluster of densely interconnected genes based on co-expression. The co-expression matrix was determined by Pearson correlation between two genes and transformed into an adjacency matrix, using a soft threshold determined according to the scale-free topology criterion using the following formula: where cor(x i ,x j ) is the correlation of any possible gene pair and b ist the soft-threshold power to generate a weighted network adjacency matrix. This matrix was used to determine the topological overlap metric, which was then further used to cluster genes and miRNAs according to their expression trajectories.
The multiSetMEs function was used to obtain metrics about the correlation and the preservation of the consensus modules across the three cell lines and the inter-module correlation within each cell line. Hub genes were identified extracting the kModule Eigengene for each gene in each module. The grey module represented genes that could not be assigned to a specific module and was therefore excluded from further analysis.
Supercluster Analysis and Singular Value Decomposition: Identification of Eigengene modules generated a high number of modules (49 modules → mRNA, 17 modules → miRNA), of which subsets exhibited similar gene-expression patterns. To collect modules with correlated gene expressions, hierarchical clustering using the "average" agglomerative clustering algorithm (mRNA distance: 3.2, miRNA distance 3) was used on modulewide averaged z-scored standardized variance stabilized (vst) gene expression counts, which were further averaged across the three cell lines for genes and miRNAs. Furthermore, a second sanity check approach was employed using singular value decomposition on z-scored standardized vst count values averaged across the cell lines, according to Bennett et al. [104] using the python framework "scikit-learn" by decomposing the gene and miRNA expression matrices.
Gene Pathway Enrichment and Disease Enrichments: Every module and each supercluster underwent gene pathway enrichment analysis as well as disease enrichment analysis (DOSE, DisGeNet, www.disgenet. org). g:Profiler (https://biit.cs.ut.ee/gprofiler/gost, Ensembl 100, Ensembl Genomes 47) was used for the 4 annotation spaces Gene ontology (GO) biological process (GO:BP), GO cellular components (GO:CC), GO molecular function (GO:MF) and KEGG by reducing the number of significantly enriched pathways to best per parent by a custom written script, as described previously. [77,105] To further reduce the complexity of gene pathway enrichments, only the top 5 GO:BP, GO:MF, GO:CC, and KEGG pathways were visualized using the negative log10(p-value).
Reversed miRNA Targeting Prediction: To evaluate potential miRNA targets implicated in the differentiation of iPSC-derived nociceptors, each mRNA::miRNA time trajectory (z-score standardized vst counts) was correlated using Pearson correlation. Additionally, the databases StarBase, [108] miRTarBase, [109] and TarBase [110] were used to determine the validation status, and DIANAs microT-CDS algorithm [111] was used to predict a potential interaction. Moreover, miRNAs were ranked in percentiles based on their baseMean expression during the differentiation, to take the amount of available miRNA molecules for targeting into consideration. All the described parameters were used to determine an overall target score, calculated with the following formula: Target score = ((p + v) * r * rank) * 100 (2) where p is the microT-CDS v5 score considered only if score > 0.9, v is a boolean variable where 1 was validated and 0 was not validated yet, r is the Pearson correlation coefficient between the miRNA and the gene and rank is the percentile rank of the miRNA based on its expression. [77] miRNA Interaction Networks: mRNA::miRNA interaction network analysis, for miRNAs known to play a prominent role in the regulation of iPSC-derived nociceptor differentiation, was performed for each expressed miRNA (BaseMean > 5 counts in at least 3 samples).
Only gmRNA::miRNA interactions that were highly negatively correlated (r < −0.7) were taken into consideration. Those genes that were at least validated or predicted (DIANA microT-score > 0.9) with an overall target score > 50 if not differently stated were used for further stage enrichment analysis. Target gene PPI was identified using StringDB with the confidence score set to 0.7 if not stated otherwise and a chord plot with genes as nodes, the arc color as supercluster affiliation and the edges as PPI was constructed using "Holoviews v1.13.2," "Matplotlib v3.4.2" and "Bokeh v2.0.1." To identify potential supercluster stage enrichments the number of targeted genes for each supercluster was determined and compared to the global number of possible targeted genes per differentiation stage. A contingency table was constructed and a 2 test performed. To evaluate which modules were potentially enriched, standardized Pearson Residuals were determined and a Pearson Residual cutoff > 2 was used to distinguish highly enriched differentiation/module stages. To identify the regulation of functional pathways by single miRNAs g:Profiler analysis of the targeted genes was performed.
For ion channel specific mRNA::miRNA interactions a spring-forced network was constructed using the target-score as edges weight. Networks were drawn with the python module "NetworkX" and refined with Cytoscape and the NetworkAnalyzer (https://cytoscape.org). www.advancedsciencenews.com www.advancedscience.com Hub Genes::Hub miRNAs Interactions: Genes considered hub genes were characterized by the kME of above > 0.8 and a p-value < 0.05. Hubgene networks for each module were constructed using the top 30 hubgenes ranked based on the kME and the p-value for all three cell lines and the edges weight was defined as the correlation between two genes calculated using "scikit-learn v0.24" stats.pearsonr() function.
Only edges with a correlation > 0.8 were considered and networks were constructed using NetworkX using the kamada kawai layout. Hub miR-NAs were miRNAs that target hub-genes. To determine hub genes that share the same miRNA interaction space, an adjacency matrix using the target score (>70) for each miRNA::mRNA interaction was constructed and dimensional reduction along the gene axis was performed using the python module unifold manifold approximation and projection "UMAP" (n_neighbor = 20, min_dist = 0.1, metric = "manhattan"). To cluster genes potentially targeted by the same set of miRNA, kMeans clustering (n_cluster = 8) was performed using "scikit-learn v0.24" kmeans() function. g:Profiler enrichment analysis was further performed to determine k-Means cluster enriched for neuronal pathways.

Whole-Cell Patch Clamp Analysis:
Single cell voltage-clamp as well as current-clamp recordings were performed using the whole-cell configuration of the patch clamp technique. Cells were kept in extracellular solution (ECS) containing (in millimolar); NaCl (150), KCl (5), CaCl 2 (2), MgCl 2 (1), HEPES (10), Glucose (10), and the pH was set to 7.3 with NaOH. Borosilicate glass pipettes (Science Products) were pulled with a horizontal puller (P-97, Sutter Instrument Company) with an average pipette resistance of 5 MOhm after filling with intracellular solution (ICS) (in millimolar): K-Gluconate (98), KCl (59), CaCl 2 (0.5), MgCl 2 (2), EGTA (5), HEPES (10), MgATP (2), NaGTP (0.2), and pH adjusted to 7.3 with KOH. Serial resistance (R-series) was compensated in all recordings using the automatic compensation function provided by the HEKA Patchmaster (>60%, 10 μs). Recorded neurons were held at −60 mV for voltage-clamp recordings and 10 mV depolarizing voltage pulses applied from −60 to 40 mV to assess peak (minimal current value per step) as well as sustained current amplitudes (mean value per step). The current-clamp mode was used to retrieve the resting membrane potential at 0 pA. The minimum current to induce an action potential (Rheobase) was obtained using increasing 5 pA depolarizing pulses. Action potential parameters analyzed were duration (from threshold to repolarization), repolarization, threshold onset, and half-width using the integrated AP fitting function of the HEKA Fitmaster. All recordings were performed in AD3 and AD2 at D36 platted on Matrigel coated IBIDI dishes maintained in N2/B27 Media supplemented with NGF, BDNF, and GDNF.
Statistics: Statistical analysis was performed using the Python packages sklearn, statsmodel, numpy, pandas, networkX, and requests as well as the Bioconductor R packages DOSE and WGCNA. Differential gene expression analysis was performed using a hierarchical model approach implemented in the Bioconductor package DESeq2 with the threshold of significance after FDR correction set to FDR < 0.05 (mRNA) and FDR < 0.1 (miRNA).
WGCNA consensus module affiliation was calculated by means of the WGCNA package distributed by Bioconductor, using kME and metakME and the module associated p-value for a significant affiliation was set to < 0.05. 2 test was performed following contingency table analysis with the p-value threshold set < 0.05. Pearson residuals were calculated using Python's "statsmodel v0.12.2," and Pearson residuals above and below abs(2) were considered as indicator for enrichments/depletions. Correction for multiple comparisons was performed where appropriate using the non-negative Benjamini-Hochberg FDR correction.
Key Resource Table: A key resource table as well as software used for the project is deposited in Table S21, Supporting Information. www.advancedsciencenews.com www.advancedscience.com

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.