NOCICEPTRA2.0 - A comprehensive ncRNA atlas of human native and iPSC-derived sensory neurons

Summary Non-coding RNAs (ncRNAs) are pivotal in gene regulation during development and disease. MicroRNAs have been extensively studied in neurogenesis. However, limited knowledge exists about the developmental signatures of other ncRNA species in sensory neuron differentiation, and human dorsal root ganglia (DRG) ncRNA expression remains undocumented. To address this gap, we generated a comprehensive atlas of small ncRNA species during iPSC-derived sensory neuron differentiation. Utilizing iPSC-derived sensory neurons and human DRG RNA sequencing, we unveiled signatures describing developmental processes. Our analysis identified ncRNAs associated with various sensory neuron stages. Striking similarities in ncRNA expression signatures between human DRG and iPSC-derived neurons support the latter as a model to bridge the translational gap between preclinical findings and human disorders. In summary, our research sheds light on the role of ncRNA species in human nociceptors, and NOCICEPTRA2.0 offers a comprehensive ncRNA database for sensory neurons that researchers can use to explore ncRNA regulators in nociceptors thoroughly.


INTRODUCTION
Small non-coding RNA (ncRNA) species are increasingly emerging as important regulatory hubs in health-related biological processes and the pathogenesis of a multitude of diseases including pain disorders.These short ncRNA species include microRNAs (miRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNAs (snRNAs), tRNAs, rRNAs, PIWI-interacting RNAs (piRNAs), and small interfering RNAs, ranging from smaller molecules such as mature miRNAs ($22 nucleotides) to full-length sno-and tRNAs ($70-120 nucleotides). 15][6] In addition, fragmented mature ncRNAs such as tRNAs and snoRNAs can be incorporated into the Argonaute complex thereby mimicking the function of miRNAs. 7,8For example, upregulation of tRNA fragments (tRFs) inducing neuronal necrosis is detectable in post stroke patient blood samples. 9,10ramatic changes in the expression of coding genes, ncRNAs, pseudogenes, and splice isoforms are seen during the transition from pluripotent stem cells to early differentiating neurons and several long ncRNAs (lncRNAs) undergo dramatic expression changes suggesting important roles for ncRNAs in neurogenesis. 11,12Distinct ncRNA expression signatures are retrieved from different cell types indicating that ncRNA signatures can represent and resolve a particular phenotype. 4,6,13However, cell-specific ncRNA signatures are still in their infancy and temporal developmental dependencies of cell-type specific markers derived from single-cell and bulk RNA sequencing studies are underinvestigated or even missing for ncRNA species other than miRNAs.ncRNA transcriptome analyses are usually more complex compared to mRNA sequencing approaches due to relevant caveats in the alignment and counting process of ncRNAs, mostly correlated with shorter nucleotide length and the tendency of multiple occurrence and annotations within the genome.This specifically applies to shorter variant ncRNA species such as tRFs, snoRNAs, piRNAs, and miRNAs. 146][17] Thus, our understanding of ncRNA signatures and their roles throughout human neuronal development remains limited, and neither temporal single-cell transcriptomic data nor bulk-ncRNA transcriptomic datasets specific for human neuronal development are available.
Nociceptors are primary afferent sensory neurons serving the transduction of noxious stimuli leading to the perception of pain.Although rodent models are currently the gold standard for functional sensory neuron research, single-cell RNA studies identified substantial differences between human and rodent sensory neurons, indicating sex-and species-specific functionalities with potential relevance for the forward and reverse translation of findings from rodent models toward human pain disorders. 18,19Recent advances in the reprogramming of induced pluripotent stem cells (iPSCs) allow to develop human iPSC-derived nociceptors which offer improved opportunities for modeling human pathophysiologies in depth. 20,21evelopmental trajectories of mRNAs and miRNAs characterize five different stages of nociceptor differentiation from human iPSC 22 and likely play important roles in developmental processes via direct or indirect regulatory pathways.Tissue and cell-specific expression of ncRNAs is associated with cell type functionalities, phenotypes, and morphologies; however, mechanistic insight into their functionality remains sparse. 4Despite their critical importance for health and disease, developmental signatures of ncRNAs are mostly explored for miR-NAs 22,23 whereas unbiased knowledge on other ncRNA species relevant for human nociceptor development is currently not available.Therefore, we for the first time provide a comprehensive atlas of ncRNA signatures, which potentially contribute to human iPSC-derived nociceptor development and maturation, with a specific focus on individual ncRNA expression trajectories, their importance in specific developmental stages, and ncRNA-mRNA connectivity.

Time dependency of small RNAs throughout human iPSC-derived sensory neuron differentiation
In this study, we report on trajectories and hub ncRNAs to provide the first ncRNA atlas throughout the development of human iPSC-derived nociceptors (iDNs).Since our recent analysis 22 of mRNA and miRNA interactions excludes other relevant ncRNA species, we now extended the analysis toward tRNA, snoRNA, rRNA, lncRNA, piRNA, and long intergenic ncRNA (lincRNA) expression patterns during iDN development from three different iPSC clones.To retrieve counts of different ncRNA species, trimmed small RNA sequencing reads were counted and aligned using the recently published eXCerpt pipeline against the human genome (hg38), miRBase.org,gtRNAdb and piRNAdb database in the eXCerpt pipeline's default order.General raw count frequency analysis across samples revealed that the majority of ncRNA reads were annotated as miRNAs, followed by tRNAs, miscellaneous RNAs (miscRNAs), and snoRNAs and further mapped to other annotations such as intron retentions and protein-coding genes (Figure S2A; [24][25][26][27] ).In addition, distributional frequencies per biotype were also mostly conserved across samples and clones within a time point (Figure S2A).Aggregation on time point revealed that miRNAs exhibited low distributional changes in the frequency of counts throughout differentiation, where only day 26 showed a $10% reduction in the relative number of counts compared to all other time points (Figure 1A).
In contrast, tRNAs and snoRNAs exhibited pronounced time point-dependent frequencies, and tRNAs accounted for 10.93% and 14.04% of raw counts in the iPSC stage and early differentiation phases.Lower frequencies were identified at D09 (5.68%) and D16 (7.38%), and a small increase in tRNA count frequency was observed in more mature sensory neurons at D26 and D36 (Figure 1A).While snoRNAs showed a general decline in frequency throughout sensory neuron development, with 8.68% of counts annotated at D0 and only 5.97% of counts annotated at D36, miscRNA showed anti-correlated patterns with low frequencies in the iPSC stage (2.69%) and high number of counts at D26 (16.34%) and D36 (11.49%, Figure 1A).The majority of annotated tRNAs, rRNAs, snoRNAs, mt-tRNAs, and miRNAs maintained at least low expression levels throughout development (Figure 1B).Differential gene expression analysis using the model $clone*time point (reduced model: $clone) revealed that the majority of expressed ncRNAs were differentially expressed (8,000 aligned ncRNAs; 3,000 with a base mean above 20) throughout development at any of the time points (Figure 1B, right panel).To evaluate the most significantly changed ncRNAs, differentially expressed ncRNAs were ranked based on the negative log10 adjusted p value.In pluripotent iPSCs (D00), the lincRNAs CTD-244A18.1 and LINC00678, the tRNA-Thr-AGt-3, and SNORD90 were enriched.The piRNAs hsa_piR_018165 and hsa_piR_016659 were highly enriched at D16 (Figure S1A).Interestingly, several snoRNA also exhibited abundant expression at D16 such as SNORA80B and SNORA73B with currently unknown functions for neuron development (Figure S1A).As most regulated ncRNAs, tRNA-Glu-CTC-1 (upregulated at D26-D36), tRNA-Gly-GCC-1 (upregulated at D5), and hsa_piR_016658 (upregulated at D16) were highly abundant and significantly expressed throughout development (Figure S1B).
To further evaluate the general time dependency of ncRNAs during development, principal component analysis (PCA) revealed a high degree of variance (PC1, 66% variance ratio) explained by the different time points, which implicates pronounced ncRNA changes throughout the development of iDNs (Figure 1C).In contrast, much smaller clonal variances were detected, indicating that developmental changes were much more drastic than clonal differences (Figure 1C).This supports the idea that conserved biological ncRNA expression signatures may serve specific developmental stages during neuron differentiation and maturation.
In addition, human dorsal root ganglion (hDRG) samples from seven organ donors were analyzed, and distributional analysis revealed inter-sample variance of ncRNA biotypes, where miRNAs (71.24%) and tRNAs (18.44%) were the predominant expressed biotypes, followed by snoRNAs (4.32%) and snRNAs (4.09%) (Figure 1D).The pronounced distributional variance between the samples might be correlated to donor metadata information, however, currently, the sample size is too small to infer confounding donor covariates (Figure S2B).
We used correlation analysis of normalized, variance stabilized counts to infer how similar iDN samples were to hDRG samples at different time points.iDN ncRNA signatures became increasingly correlated with hDRG samples throughout the differentiation process (Day 00 R = $0.69,Day 36 R = $0.84; Figure 1F).Although inter-sample variance of hDRG samples was relatively high, correlation analysis with iDN samples revealed conserved similarity scores across all hDRG samples indicating similar intra-sample ncRNA expression signatures between hDRG samples (Figures S2C and S2D).hDRG and late-stage iDNs shared several top expressed miRNAs such as the let-7 family which were highest expressed in hDRG, and miR-182/183, which were the top expressed miRNAs in iDNs (Figures S3A and S3B).Trajectories of the let-7 family showed a strong upregulation towards more mature iDNs (>Day26, Figure S2F) supporting a critical role of the let-7 family in maturation.miR-21, miR-100, miR-30d, and miR-30a were consistently found among the list of top-20 expressed miRNAs in hDRG and all developmental iDN stages (Figures S3A and S3B).ncRNA signature differences were likely explained by the fact that hDRG samples not only contained sensory neurons but also other cell types such as Schwann cells, satellite glia cells, and immune cells whose ncRNA profiles differ considerably from neuronal ones.

ncRNA modules with preserved trajectories
To classify ncRNAs into highly interconnect modules representing co-expressed units, we filtered ncRNAs (baseMean >20 as a metric of abundance and to remove spurious correlations) and used weighted gene correlation network analysis (WGCNA).This revealed 14 modules of highly correlated ncRNAs, and allowed us to identify modules with conserved expression trajectories found in all three iPSC clones such as the black, yellow, green, and blue module.Conversely, we also detected modules associated with expression changes of ncRNAs found only in one of the iPSC clones indicating clonal effects (e.g., magenta and red module) (Figures 2A and 2B).Module trajectories of correlated ncRNAs were obtained using Eigengene vector calculations (WGCNA R package) and associated with all stages of iDN development as defined previously 22 (Figure 2A).Especially, the black and yellow modules were associated with the iPSC and early differentiation stages and hub ncRNAs were classified as tRNAs such as tRNA-Thr-AGT-3 and tRNA-Val-TAC-4 (both black module), tRNA-MET-CAT-chr6 and tRNA-Lys-TTT-7 (both yellow module), as well as snoRNAs such as SNORD88B (yellow module), and lincRNAs such as RP3-467D16.3 (black module) and MIR17HG (yellow module).In contrast, hub snoRNAs, such as SNORD90, SNOR24, SNORD116, and SNORD32A, associated with the green and blue modules, were highly upregulated during later stages (D16-D36).Of these, SNORD116 variants are known as important regulators of neuronal differentiation (Figure 2B). 28Interestingly, the majority of ncRNAs associated with the brown module were annotated as Y RNAs (fragments).These Y RNAs such as Y_RNA.439and Y_RNA.700 were highly enriched from D05-D16 implicating putative functions in neural crest cells and progenitor cells developing into a sensory neuron fate (Figure 2B).(B) Top 20 ncRNA hub gene network determined using determining the Pearson correlation coefficient between all module members using variance stabilized counts.Networks were drawn using the NetworkX python module using the correlation between ncRNA nodes as weight and the Netgraph module with the spring layout and the correlation score as edges weight.

Differential temporal distribution of tRNA fragments during sensory neuron development
tRNAs represent a well-known ncRNA species serving RNA translation into protein; however, novel regulatory functions for tRNA fragments (tRFs) have recently been reported. 29,302][33][34] tRNAs and tRFs exhibited pronounced time-dependent expression patterns during iDN differentiation supporting the idea that tRNA regulatory pathways contribute to neurodevelopmental processes (Figure 1A).A more in-depth analysis of count frequency for amino acid specific tRNAs revealed that tRNA Gly and tRNA Glu were the most frequent tRNA fragments throughout differentiation with inversely correlated frequency patterns: while tRNA Gly was more abundant early in development, tRNA Glu predominated in transition stages such as D09 and was highly abundant in late iDN stages (Figure 3A).Most counts (average read length $25 nucleotides long) emerging from the tRAX pipeline represented tRFs rather than full-length tRNAs, which may partially result from potential limitations of the underlying sequencing technology (Figure S4A). 35However, the percentages of tRNA Lys (7.94%) and tRNA Asp (6.94%) and tRNA His fragment counts were strongly increased at D09 in comparison to all other days, indicating a putative shift not only in individual tRNA species but also tRNA fragment types, and this finding supports specific regulatory roles of tRFs to control gene sets required for pluripotency vs. neuron differentiation and maturation.
To explore the dependencies between neuronal development and tRNA expression signatures, two approaches were used: first, only tRFs determined by eXcerpt were used for analysis (Figures 3B-3D); second, tRFs were realigned using tRNA-specific tools such as MintMap2 to reclassify tRNA fragments based on mapping occurrence (Figure 3E).eXcerpt normalized counts of annotated tRNAs followed by PCA revealed pronounced time dependency for tRNA fragments (Figure 3B).To retrieve information about co-expressed tRFs, all expressed tRNA sequences were subject to hierarchical clustering which revealed six stages characterized by tRF expression during sensory neuron differentiation (Figure 3C).By ranking of each cluster based on the FDR-corrected p value and plotting of trajectories, we found top-regulated tRFs for each cluster such as tRNA-Lys-TTT-9-1 (iPSC stage, early differentiation, cluster 1), tRNA-SeC-TCA-2 (nociceptor stage, cluster 4), tRNA-Gly-GCC-2 (nociceptor stage, cluster 5), and tRNA-Gly-CCC-5-1 (early differentiation phase/nociceptor stage, cluster 6, Figure 3D).In addition, MintMap2 was used to classify the exact fragment type of the exclusive tRF spaces, and distributional analysis throughout iDN differentiation was performed. 36Most reads classified as 5 0 tRFs (5 0 halves, 5 0 tRF, as supported by the tRAX analysis) followed by 3 0 tRFs (3 0 halves, 3 0 tRFs) throughout differentiation, but also showed proportional changes between 5 0 halves and 5 0 tRFs.While 5 0 tRFs were enriched in iPSCs, 5 0 halves became more abundant throughout differentiation, which might indicate differential cleavage of tRNA fragments throughout development.In addition, internal tRFs showed the highest compositional fraction at D09 (Figure 3E).To further evaluate putative roles of the most significant tRNA fragments per cluster, we used the DIANA MR-microT software and queried the aligned sequences derived from the MintMap2 alignment for putative target mRNA spaces.The tRNA GlyCCC was likely targeting genes implicated in Axon Guidance (KEG:04360) and GluCTC appeared to regulate cell projections and synapse-related pathways (Table S5).Interestingly, both tRNA-GluCTC as well as tRNA-GlyCCC significantly decreased at Day 09, which is considered the starting point of extensive axon growth.This might support the idea that these two tRNA fragments may play a significant role in axonogenesis/synaptogenesis (Figure S5A).In contrast, ThrAGT and ThrTGT were likely involved in RNA transcription and metabolic process regulation (Table S5).
We also evaluated tRNAs likely mimicking miRNAs and queried tRNA-targeted genes for each individual tRNA fragment to retrieve miRNA interactions using microTCDS v.5.GlyCTC and GlyCCC mimicked the functionality of the miR-548 family and CycACA and AlaAGC shared target spaces with miR-3163, miR-4282, and miR-4536-5p (Figure S5B).In addition, overlap analysis of hDRG and iDNs revealed that a significant portion of the top-20 expressed tRNA fragments were found both in hDRG and iDNs such as members of the tRNA-GlyGCC, tRNA-GluCTC, tRNA-AspGTC, and tRNA-LysCTT family (Figures S3A and S3B).
Finally, we evaluated expression changes of mitochondrial (mt-) tRNAs throughout differentiation and found that all mt-tRNAs were differentially expressed with increasing expression throughout sensory neuron development, which putatively correlated with an increased expression and translation of mitochondrial proteins to secure the increasing energy requirements of maturating neurons (Figure S4). 37,38fferential temporal distribution of snoRNAs in sensory neuron development snoRNAs are well accepted for their role in ribosomal protein regulation but expression signatures and regulatory roles in human neuronal development are mostly unknown. 3,39PCA on all expressed snoRNAs revealed high temporal dependencies of snoRNAs indicating a strong regulatory effect ofnsnoRNA expression throughout differentiation (Figure 4A).Six clusters emerged with differential temporal snoRNA expression signatures associated with different developmental stages, with cluster 1 and 2 representing the late stages of development, stage 3 the progenitor stage/immature neuron stage, clusters 4 and 5 the early developmental stages, and cluster 6 which represented snoRNAs downregulated in iPSCs (Figure 4B).SNORD93 (baseMean: 1618 normalized counts), SNORA69.1 (baseMean: 37 normalized counts), and SNORD88B (baseMean: 74 normalized counts) were the most significant differentially expressed snoRNAs associated with the iPSC or early development stages while SNORA80B (base-mean: 25 normalized counts), SNORA73B (base-mean: 731 normalized counts), and SNORD94 (base-mean: 281 normalized counts) were significantly associated with D16.Their roles and the functional consequences of this regulation are currently unknown.SNORD71 (baseMean: 8436), SNORD115 (baseMean: 40 normalized counts), and SNORD116 (baseMean:1629 normalized counts) were upregulated in more mature neuron-like cells (Figures 4B and 4C).1][42] The precise function and regulated pathways are currently unknown for the majority of snoRNAs such as for abundantly expressed SNORD71, SNORD89, and SNORD90, which both exhibited a gradual increase in expression throughout differentiation (Figure 4C) indicating that specific regulatory roles may apply to other snoRNAs as well.This is further supported by our finding that the majority of the top 20 expressed snoRNAs were conserved (E) Sequencing reads were further aligned using MintMap2, a specific pipeline to count tRNA reads and determine location of tRNA fragments, while also revealing specific sequence of ncRNA contigs as well as the specific fragment type.Fragment types per time point were aggregated using the sum function within samples and the average function across samples from the same time point.
between hDRG and iDNs such as SNORD30, SNORD104, SNORD48, and SNORD69.Since snoRNAs modifying RNAs via methylation (C/D) (C/D Box snoRNAs) are distinct from snoRNAs which introduce pseudouridination (H/ACA) (H/ACA snoRNAs), we investigated the distribution of the two snoRNA classes throughout development (Figure 4D).Significantly regulated C/D box snoRNAs (cluster 1, 2) were enriched during late differentiation (Figure 4D), while expression H/ACA snoRNAs was higher in the early stages of iDN development (cluster 5).Since functional consequences of snoRNAs are difficult to depict, we focused further on the expression of host genes (gene that is harboring the ncRNA) associated with snoRNAs and correlated host gene with snoRNA expression.Significantly correlated and anti-correlated host genes with R > abs[0.7]were subject to ontology enrichment analysis.This revealed detrimental changes in ribosomal regulation throughout sensory neuron development for both positively as well as negatively correlated snoRNA host genes as indicated by the expression patterns of ribosomal host genes (Figures S4B and S4C; Tables S3/S4 positive correlated/negative correlated).

Temporal regulation of rRNAs, snRNAs, and piRNAs
We also evaluated the temporal trajectories of rRNAs, snRNAs, and piRNAs using PCA.snRNAs and rRNAs showed only minor time-dependent changes throughout nociceptor development (Figures 5A and 5B).Nonetheless, clustering analysis revealed five clusters associated with snRNA and three clusters associated with rRNA expression (Figure S6).Clonally conserved RNU5E-6p and RNU7-3P ncRNAs were highly expressed early in development (associated with cluster 1) and decreased with increasing differentiation and maturation indicative of a regulatory role in suppressing pluripotency genes.RNU2-3P, RNU5E-1, and RNU7-19P expression increased from Day 09 to Day 36 (cluster 1, 3, Figure S4B).In contrast, piRNA fragments with a usual length of 21-24 nucleotides exhibited pronounced temporal dependency and clustering based on time point associated with sensory neuron development (Figure 5C).Small silencing piRNAs are implicated in the silencing process of transposable elements in germline cells thereby ensuring germline DNA integrity. 2Although the cumulative number of piRNA reads was low, hierarchical clustering revealed five well-distinguishable clusters of piRNA expression signatures (Figures 1B, 5D, and 5E).Cluster 1 and 2 piRNAs with the most significant members hsa_piR_011968 and hsa_piR_019628 exhibited abundant expression in the iPSC stage and early during development.Cluster 3 piRNAs (hsa_piR_018165) were upregulated at D09-D16 and cluster 4-5 piRNAs (hsa_piR_00170 and hsa_piR_002158) showed abundant expression at late stages (D26-D36, Figures 5D and 5E).Interestingly, the overall percentage of piRNA counts was considerably lower in hDRGs compared to iPSC (0.2% in hDRGs and $1% in iDNs).However, the highest expressed piRNAs showed a strong overlap between hDRG and iDNs, and the highest expressed piR-016658, piR-017033, and piR-018780 were conserved (Figures S3A and S3B).Based on the specific expression patterns and their trajectories, regulatory roles in neurodevelopment can be anticipated also for piRNAs.

lncRNAs
We further investigated lncRNAs with and without poly-A tails by integrating both the mRNA with the small ncRNA sequencing dataset.lincR-NAs derived from small RNA sequencing and lncRNAs from long RNA sequencing exhibited pronounced temporal dependencies during iPSC differentiation into sensory neurons (Figures 6A and 6B).Clustering analysis revealed 5 lincRNA groups and 7 lncRNA groups with distinct temporal signatures (Figures 6C and S4A).lncRNAs B3GALT5-AS1 (ENsG00000184809) and LINC01356 (ENSG00000215866) were highly enriched in the iPSC stage.In contrast, the lncRNAs LINC00205 (ENSG00000223768), and FTX (ENSG00000230590) exhibited signatures of nociceptor stage-related gene expression, 43 support an important role of FTX in the survival and possibly function of mature human iPSC-derived nociceptors.In addition, SOX1-OT (ENSG00000224243) and RENO1 (ENSG00000287431) were highly abundant in cluster 2 ($Day05-Day09).5][46] Based on these reports and our current findings, functional relevance is anticipated for SOX1-OT and RENO1 trajectories for iDN fate decision and differentiation.Since overall count numbers for lincRNAs were low, we merged lincRNA and lncRNA clusters for further downstream analysis.In general, lncRNAs/lincR-NAs act as sponges for miRNAs or other RNA fragments in neurons. 47Therefore, we specifically searched for enriched miRNAs that putatively targeted lncRNAs per cluster.lncRNAs associated with pluripotency were enriched for miRNAs of the let-7b family, which is increasingly expressed throughout development and essential for neural stem cell differentiation and cell fate commitment. 22,48,49lncRNAs abundant at D00-D05 were enriched for hsa-miR-16-5p, hsa-miR-15a|b-5p, and hsa-miR-434, while at the mature nociceptor stage (D36) enriched lncRNAs targeted hsa-miR-574-5p, hsa-miR-485-5p, hsa-miR-369-3p, and hsa-miR-143-3p.Especially lncRNAs enriched at D00 and D05 showed highly anti-correlated trajectories with their enriched miRNAs, suggesting putative regulatory controls of these miRNAs in pluripotent cells by interactions with the enriched lncRNAs or vice versa (Figure 6D).
Finally, we compared iDN lncRNA expression signature to hDRG signature by means of a correlation analysis between all stages of iDN development and hDRG samples prepared from lumbar L1-5 and thoracic T12 DRG.This revealed that maturating iDNs became more like hDRG neurons in terms of lncRNA expression signatures.Interestingly, iDNs were strongly correlated to neurons located in lumbar segment L2 and L3 hDRG whereas weaker correlation was identified for iPSCs (D0, Figure 6E).

DISCUSSION
While we previously determined miRNA trajectories and hub-miRNAs, this work extends the landscape of ncRNA regulation by revealing that $40%-50% of the retrieved sequences are not classified as miRNAs but as other non-coding RNA species such as tRNA fragments, snoRNAs, or miscRNAs in human iPSC-derived sensory neurons.Our distributional analysis of all currently known ncRNA species revealed conserved clonal temporal changes in count distribution at different time points throughout development for tRNAs, snoRNAs, piRNAs, lncRNAs, mis-cRNAs, and lincRNAs.This putatively indicates that ncRNAs in general and not only miRNAs are functionally important and need to be considered as phenotypic cellular signatures of pluripotent cells and differentiating human nociceptors development in greater detail.
Several ncRNA species including miRNAs, snoRNAs, and lncRNAs are emerging as important cell and tissue-specific hub regulators of processes involving protein-coding genes and other RNA products. 2,7,9,50We carried out ncRNA transcriptomics from native human DRG to provide insight in the distribution and expression signatures of different ncRNA biotypes and understand similarities and dissimilarities between native hDRG and iDNs.The combined analytical approach of iDN and hDRG putatively helps to identify ncRNAs that are relevant for human sensory neuron differentiation and maturation.This is highly relevant since in-depth analysis of ncRNAs is not possible through scRNA sequencing technologies yet.We developed an atlas for ncRNA expression of nociceptors differentiated from human iPSCs and revealed that additional ncRNA biotypes next to miRNAs can contribute significantly to the developmental process in human stem-cell derived neurons. 5,51To enable researchers to re-use our data and investigate ncRNAs in depth, our significantly extended NOCICEPTRA2.0tool will be provided locally at GitHub (https://github.com/MUIphysiologie/NOCICEPTRA2.0) or online via Streamlit cloud (https://nociceptra. streamlit.app)as an open-source web framework connected to the snowflake cloud, to enable easy accessibility to analyzed data.NOCICEP-TRA2.0 is the first comprehensive ncRNA atlas to explore all currently known RNA species and their trajectories in detail, together with WGCNA modules for mRNA/miRNA and ncRNAs as well as indications of ncRNAs differentially expressed throughout development.
Five to seven developmental mRNA::miRNA clusters characterize specific stages of human nociceptor development, 22 and these were also conserved for the majority of ncRNA biotypes.The majority of ncRNAs showed time-dependent expression trajectories throughout development, suggesting that specific ncRNAnome signatures are necessary to drive human iPSC-derived nociceptor development.Interestingly, few studies report detrimental changes in snoRNAs and tRNAs in differentiating mouse embryonic stem cells. 52,53Our data confirmed increasing expression following differentiation initiation until D16 for Snora53, Snora54, and highly expressed tRNA-Gly-GCC early in development supporting important roles for Snora53/54 throughout the entire development process for neural progenitor cells whereas tRNA-Gly-GCC appeared to be involved in the general initiation process of differentiating neurons.
The majority of significantly regulated ncRNAs of different biotypes have not been reported or functionally associated with neuron development.This includes tRNA fragments such as tRNA-Thr-AGT-3 and tRNA-Ala-AGC-20, the snoRNAs SNORD71 and SNORD90, the piRNAs hsa_piR_019420 and hsa_piR_018165, and a whole group of miscRNA as most significantly regulated ncRNAs.tRNAs most likely represented tRNA fragments (cleaved tRNAs) while complete tRNA transcripts were rare probably due to the sequencing approach. 54While the majority (D) lncRNA and lincRNAs clusters were merged into one cluster holding both groups, based on trajectories and assigned lincRNAs and lncRNAs and subjected to the lncSEA algorithm to determine potential enrichment of miRNAs per cluster.Enriched miRNAs per cluster were determined based on the number of lncRNAs interacting with the miRNA (dot-size).(E) Correlation analysis using variance stabilized counts of lncRNAs at all time points in iDN development and in hDRG samples revealed a significant increase of similarity in lncRNA expression throughout iDN development.
of tRNAs were annotated as 5 0 and 3 0 tRFs, differential cleavage throughout the development may account for shifting from tRFs to longer tRNA halves, indicating alterations in tRNA silencing and target spaces or a compensatory response to stress throughout development as tiRNA (tRNA halves) are upregulated in response to stress. 9,55In addition to the differential composition of fragment types throughout differentiation, we also showed that tRNA-isodecoder distribution shifted from glycine to glutamine throughout differentiation, which might indicate an increased usage of glutamine for translational processes in later stage iDNs.
Since it is well described that interactions of ncRNA fragments with mRNAs as well as target spaces are many-to-many relationships, we clustered all ncRNA biotypes into transcriptional units exhibiting the same expression signatures. 56Modules highly correlated with late-stage, more mature nociceptors consisting of multiple ncRNA biotypes such as SNORD90, which is enriched in brain tissues, orchestrating the development into iDNs. 4In future, it would be interesting to integrate these modules with temporal ATAC sequencing data to explore putative transcription factor binding sites correlated with the expression trajectories of ncRNAs. 57,58n addition to the short ncRNAs, lncRNAs play important roles in neuron development by regulating cell fate determination, progenitor pool maintenance and proliferation, as well as neuronal survival. 46,59Non-surprisingly, lncRNAs exhibited the highest temporal to clonal variance ratio throughout iDN development, indicative of conserved roles for lncRNAs not only in neurodevelopmental processes of the brain [59][60][61][62] but also in human iDN and nociceptor development.In addition, lncRNA expression signatures of iDNs became more similar to hDRG lncRNA signatures with increasing differentiation suggesting that the latter predominantly originated from sensory neurons.This bolsters the suitability of the iDN model for in vitro human sensory neuron research.Moreover, overexpression of those ncRNAs that are highly expressed in hDRGs but lowly in earlier stages of differentiation of iDNs such as the let-7 family holds great potential as promising strategy to further enhance maturity and similarity of human iDNs and native nociceptors.
Evaluating the roles and contributions of specific ncRNA biotypes throughout iDN development revealed differential importance of ncRNA biotypes for the developmental process.While tRNAs, snoRNAs, piRNAs, and lncRNAs exhibited pronounced temporal development-specific dependencies, which putatively indicated critical contributions to development processes, snRNAs and rRNAs were only to a minor extent associated with the differentiation process in our analysis.Although this might be confounded by the overall shallow depth of the ncRNA sequencing approach, which only allowed to investigate more abundantly expressed ncRNAs, it could also support the idea towards a more general role of snRNA and rRNA regulatory signatures conserved throughout all developmental stages.Changes in lncRNA, tRF, and miRNA expression patterns with increasing maturation may also reflect increased electrical activity of more mature neuron and immediate changes in expression after activity, during which fluctuations in mRNA abundance subject to ncRNA regulation are observed. 63n summary, developmental ncRNA expression signatures of iDNs were approximating signatures of native hDRG throughout development, and top expressed and regulated ncRNAs were surprisingly similar as indicated by our analysis.This strengthens the usage of iDNs for sensory neuron research and provides further opportunities to improve the maturity of iDNs.Overall, NOCICEPTRA2.0provides the first comprehensive and online available, searchable atlas of currently known RNA species in developing and maturating iPSC-derived human nociceptors which become increasingly important as human model systems for preclinical pain research.

Limitations of the study
While this study marks significant progress in unraveling putative functional roles and expression trajectories of ncRNAs in human sensory neuron development, several limitations still warrant consideration.Firstly, human DRG inherently comprise a mixture of various cell types, posing challenges for precise correlative analyses between human DRGs and iPSC-derived sensory neurons.Additionally, the degree of maturity achieved by iPSC-derived sensory neurons remains inadequately characterized.Moreover, our study encompasses data from only three distinct iPSC-derived clones, which may not comprehensively capture the full spectrum of heterogeneity among iPSC clones.To this end, we also used just a single donor for total RNA sequencing, which might not properly depict the population variance and might introduce bias in the analysis.However, the aim was rather to retrieve an idea which DRGs (spinal location) might be the most correlated to the iPSC-derived sensory neurons.In addition, we added also multiple hDRG samples for small RNA sequencing to compare the ncRNA species in depth and reduce bias.
Furthermore, the annotation of ncRNA fragments remains a challenge, as substantial overlap exists among ncRNA species of different family sequences, potentially impeding accurate annotation.While more recent analysis pipelines aim to mitigate annotation errors, certain ncRNA fragments, particularly tRNA fragments and piRNAs, often exhibit minimal sequence differences, making precise mapping a complex task.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  S1 and S2) was done through a collaboration with the Southwest Transplant Alliance and all protocols were approved by the Institutional Review Board of University of Texas at Dallas (MR 15-237 -The Human Dorsal Root Ganglion Transcriptome).RNA was isolated from lumbar and thoracal human DRG samples and total RNA was sequenced using the Illumina TrueSeq sequencing kit for this study.In addition, small RNA sequencing was performed using the NEBNext Small RNA Library Prep Set (The Donor Table at Table S2).

Human native DRG analysis
Fastq files were aligned against the human transcriptome provided by Gencode using Salmon.Obtained transcript per million reads were transformed into gene aggregated row counts using the tximport module.Libraries were normalized using DESeq2 and variance stabilized counts retrieved for lncRNAs and correlated against all the individual timepoints of iPSC-derived sensory neuron development, only considering lncRNAs.
In addition, ncRNA libraries of human organ donors DRG samples were generated using the Nextflex Small RNA-Seq kit (Table S2).In short ncRNA fastq files were quality controlled using FastQC and trimming of adapter as well as random 4-bp long nucleotides trimming from each end of the read were performed using cutadapt with parameters suggested by the Nextflex kit.Processed fastq files were then aligned and counted using the exceRpt pipeline against the hg38 reference genome and DESEq2 was used for count and library size normalization as well as generation of variance stabilized counts (vst).Finally, vst and normalized counts were compared with counts obtained from the iDN analysis.

RNA read processing and differential gene expression analysis
Derived mRNA FastQ files were pseudo aligned against the human transcriptome (Genecode.v39.transcripts.fa)provided by GENCODE (www.gencodegenes.org) using Salmon (v0.13.1) software with the following parameters (-l A, -p 16, -numBootstraps 30, -useVBOpt -seq-Bias -validateMappings).Small RNA sequencing libraries were prepared using the NEBNext Small RNA Library Prep Set for Illumina as previously described (Zeidler et al. 22 ).Small RNA reads were further processed and adapter sequences trimmed off (AGATCGGAAGAGCACA CGTCTGAACTCCAGTCAC) using Flexbar v3 and subsequently aligned against the human reference genome using the exceRpt pipeline and the Manatee Pipeline. 36,64,65Reads were aligned and counted using the exceRpt pipeline using the default parameter order aligning first to miRNAs, tRNAs, piRNAs and then against the genome and circRNAs.

Differential gene expression analysis mRNA
Salmon retrieved mRNA quant files were aggregated on gene level counts using the tximport (v1.20.09)package and the Genecode.v39.annotationfile.Differential gene expression was performed using the model $timepoint + cell line + timepoint:cell line, where time represents the different timepoints throughout development, cell line one of the three iPSC clones (AD3, AD2 and 840) and the interaction and compared this using the likelihood ratio test against a reduced model of $cell line.Log2Fold changes acquired by DESeq2 were shrinked using the

Figure 1 .
Figure 1.Basic exploration of ncRNA distribution Basic exploratory analysis of ncRNA reads aligned using the exCerpt pipeline.(B) Timepoint and ncRNA biotype aggregated raw count frequencies were calculated using the average of each biotype per timepoint.(B) Stacked-bar chart showing percentage of expressed vs. unexpressed genes (below 5 counts, left panel), and the percentage of significant genes from the subset of expressed genes as indicated by the DESeq2 analysis (right panel).(C) Principal component analysis of variance stabilized counts were used to categorize samples into cell line and time point-specific shapes and colors.(D) Distributional analysis of count frequencies aggregated based on ncRNA biotype for hDRG.(E) Set intersection analysis of the top 500 ncRNAs expressed in late-stage iDN (Day26-36) and hDRG.(F) Correlation analysis using Pearson correlation of variance stabilized counts between hDRG samples and iDN at the different time points revealed strong timedependent approximation between hDRG and iDNs ncRNA signatures.

Figure 2 .
Figure 2. WGCNA analysis of ncRNAs WGCNA analysis of ncRNA variance stabilized counts excluding miRNAs and only using ncRNAs with a base mean count above 20.(A) Eigengene vector trajectories of WGCNA-derived modules were clustered using soft-threshold correlation into classes exhibiting the same temporal signatures.(B)Top 20 ncRNA hub gene network determined using determining the Pearson correlation coefficient between all module members using variance stabilized counts.Networks were drawn using the NetworkX python module using the correlation between ncRNA nodes as weight and the Netgraph module with the spring layout and the correlation score as edges weight.

Figure 3 .
Figure 3. tRNA fragment distribution and tRNA-specific analysis of developmental trajectories (A) tRNA-specific raw counts were aggregated on time point and tRNA amino acids and count distribution was depicted as pie chart.(B) PCA of tRNA variance stabilized counts was clustered based on cell line and time points.(C) Hierarchical cluster analysis using the distance metric (1-correlation) as distance and the ward algorithm to determine clusters revealed 6 clusters of unique tRNA expression trajectories.Trajectories are represented as Z scored scaled variance stabilized counts and both, aggregated counts (boxplots) and unique counts are shown.(D) Top regulated tRNAs per cluster are depicted by means of the negative log10 p-adjusted value as well as the expression patterns.(E)Sequencing reads were further aligned using MintMap2, a specific pipeline to count tRNA reads and determine location of tRNA fragments, while also revealing specific sequence of ncRNA contigs as well as the specific fragment type.Fragment types per time point were aggregated using the sum function within samples and the average function across samples from the same time point.

Figure 4 .
Figure 4. snoRNA distribution throughout sensory neuron development snoRNA-specific analysis of development trajectories.(A) PCA of variance stabilized snoRNA counts revealed separation of time points.(B) Hierarchical clustering (using the distance 1-correlation) analysis grouped snoRNA into 6 cluster describing iDN development (C) Ranked top regulated members based on the p-adjusted FDR corrected p value of snoRNAs per cluster and expression signatures were represented as Z score variance stabilized counts using a heatmap.(D) snoRNA class affiliation was determined using the snoDB database.Number of snoRNAs members per group (C/D, H/ACA) were determined per time point summing counts per each sample individually and average across the hierarchical cluster.

Figure 5 .
Figure 5. snRNA, rRNA, and piRNA distribution throughout sensory neuron development Analysis of snRNA, rRNA, and piRNA trajectories derived from variance stabilized counts.(A) PCA of snRNA variance stabilized counts shows low temporal dependency as depicted by the distribution of samples and the variance explained.(B) PCA of rRNA exhibits low temporal dependencies based on the variance explained.(C) PCA of piRNA shows high temporal dependency and time point-dependent clustering.(D) Hierarchical clustering (distance: 1-correlation) analysis of Z scored variance stabilized piRNA counts revealed 5 distinguishable clusters.(E) Top ranked dysregulated piRNAs per cluster visualized by the negative log10 p value and expression signatures are represented by the Z scored standardized variance stabilized counts.

Figure 6 .
Figure 6.lincRNA and lncRNA distribution analysis Analysis of lncRNAs and lincRNAs from long and small RNA sequencing reads.(A) PCA from lncRNA variance stabilized counts obtained from the mRNA sequencing.(B) PCA from lincRNAs obtained from small RNA sequencing without poly-A tail.(C) Analysis of lncRNAs using hierarchical clustering revealed 7 cluster of lncRNA regulation throughout sensory neuron development.(D)lncRNA and lincRNAs clusters were merged into one cluster holding both groups, based on trajectories and assigned lincRNAs and lncRNAs and subjected to the lncSEA algorithm to determine potential enrichment of miRNAs per cluster.Enriched miRNAs per cluster were determined based on the number of lncRNAs interacting with the miRNA (dot-size).(E) Correlation analysis using variance stabilized counts of lncRNAs at all time points in iDN development and in hDRG samples revealed a significant increase of similarity in lncRNA expression throughout iDN development.

TABLE
EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS Dorsal Root Ganglia total RNA and small RNA sequencing Recovery of human DRG samples derived from organ donors (Tables d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS B Dorsal Root Ganglia total RNA and small RNA sequencing d METHOD DETAILS