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.
Keywords: development, differentiation, human, induced pluripotent stem
cells, pain, RNA‐seq resource, temporal transcriptome analysis
__________________________________________________________________
A resource is provided investigating temporal short and long RNA
dynamics in induced pluripotent stem cell‐derived sensory neuron
differentiation, with a special focus on miRNA::mRNA interactions
regulating synapse development, neurite outgrowth, and ion channels.
The containerized docker tool NOCICEPTRA allows for exploration of
temporal trajectories, Kyoto Encyclopedia of Genes and Genomes pathway
onsets, and more.
graphic file with name ADVS-8-2102354-g008.jpg
1. Background
Nociceptors are specialized primary afferent neurons that transduce and
transmit painful stimuli to alert the brain of bodily 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).^[ [36]^1 , [37]^2 , [38]^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.^[ [39]^4 , [40]^5 , [41]^6 , [42]^7 ,
[43]^8 ^] Notable differences between mice and men challenge the
applicability of rodent model systems for the investigation of human
pain pathogenesis.^[ [44]^9 , [45]^10 ^] Human‐cell‐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,^[ [46]^11 , [47]^12 ,
[48]^13 , [49]^14 , [50]^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.^[
[51]^16 , [52]^17 , [53]^18 , [54]^19 , [55]^20 , [56]^21 , [57]^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.^[ [58]^23 ,
[59]^24 ^] Protein‐coding genes and long non‐coding RNAs (lncRNAs) have
been the focus of recent gene expression studies during cell
differentiation.^[ [60]^25 , [61]^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.^[ [62]^27 , [63]^28 , [64]^29
, [65]^30 , [66]^31 , [67]^32 , [68]^33 , [69]^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 time‐course specific gene
expression signatures of human iPSC‐derived 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.
2. Results
2.1. 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 [70] 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 [71]1B,[72]C; for analysis of individual cell lines
see Figure [73]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 [74]S1,
Supporting Information).
Figure 1.
Figure 1
[75]Open in a new tab
Differential gene expression analysis, time trajectories, and pathway
enrichments. A) Protocol characteristics and timepoints for RNA
harvesting based on Chambers et al. (2014). B) Principal component
analysis and explained variance analysis of the 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^[ [76]^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.
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 [77]1D,[78]E, Table
[79]1 , Figures [80]S1B,C, [81]S5–[82]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 [83]1E,[84]F; Table [85]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 [86]1F; Table [87]S3, Supporting Information).
Table 1.
Protein‐protein interaction network of each module
Module Nodes PPI enrichment p‐value
Bisque 4 58 3.60E‐03
Black 584 2.23E‐11
Blue >2000 n.d
Brown >2000 n.d
Brown4 75 5.92E‐02
Cyan 231 2.45E‐08
Darkgreen 75 1.00E‐16
Darkgrey 137 2.04E‐06
Darkmagenta 124 1.00E‐16
Darkolivegreen 114 1.21E‐01
Darkorange 121 3.47E‐01
Darkorange2 65 3.52E‐01
Darkred 179 3.20E‐02
Darkslateblue 41 7.13E‐01
Darkturqoise 145 3.75E‐09
Floralwhite 76 1.63E‐05
Green 813 1.00E‐16
Greenyellow 324 1.00E‐16
Grey60 217 1.45E‐07
Ivory 82 2.98E‐04
Lightcyan 232 1.23E‐05
Lightcyan1 78 1.00E‐16
Lightgreen 195 3.60E‐05
Lightsteelblue1 88 2.94E‐05
Lightyellow 184 5.80E‐04
Magenta 384 1.00E‐16
Mediumpurple3 88 2.58E‐02
Midnightblue 235 1.54E‐01
Orange 154 3.82E‐05
Orangered4 76 4.21E‐14
Paleturqoise 129 1.00E‐16
Pink 351 7.24E‐05
Plum1 95 1.94E‐01
Plum2 44 2.70E‐01
Purple 307 7.37E‐08
Red 785 1.00E‐16
Royalblue 170 3.26E‐09
Saddlebrown 139 3.14E‐03
Salmon 234 1.00E‐16
Sienna3 85 1.00E‐16
Skyblue 125 1.16E‐07
Skyblue3 107 1.64E‐02
Steelblue 133 1.75E‐01
Tan 258 1.00E‐16
Turqoise >2000 n.d
Violet 133 1.00E‐16
White 152 1.09E‐07
Yellow >2000 n.d.
Yellowgreen 92 1.01E‐06
[88]Open in a new tab
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 [89]1E,[90]F). Analysis of
important pluripotency, neural crest, and nociceptor markers revealed a
general time‐dependent 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 NEUROG1 was abundantly expressed
(≈200 TPM versus ≈4 TPM at D9, respectively, Figure [91]1H). This
finding supports successful reprogramming into small‐diameter sensory
neurons.^[ [92]^35 , [93]^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 [94]1F,[95]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 [96]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 [97]S2C and [98]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 [99]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 [100]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.^[ [101]^37 ,
[102]^38 , [103]^39 , [104]^40 , [105]^41 , [106]^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.^[ [107]^42 , [108]^43 , [109]^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 [110]1F,[111]H, Tables [112]S3,[113]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 [114]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 [115]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 expression (Figure [116]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
[117]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 [118]S1D,E,
[119]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.^[ [120]^26 ,
[121]^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 [122]S1G, Supporting
Information developmental marker). Other murine nociceptor markers,
such as SYT13, GAL, TH, and PRDM2, peaked around days 9–16 (Figure
[123]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 [124]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,^[ [125]^46 , [126]^47 ^]
Figures [127]S1J,K, Supporting Information).
Given that many of the key genes are well known checkpoints in murine
sensory neuron differentiation and human sensory neurons,^[ [128]^26 ^]
our analysis indicates a largely conserved sensory neuron developmental
program from mouse to human.
2.2. 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 [129] 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 [130]2B–[131]D). Enrichment analysis indicative of the
functional relevance for each respective stage showed that
neuron‐related pathways were induced as early as in the 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 [132]2E, Table [133]S5, Supporting Information) and
several neuronal related diseases, such as cognitive disorders,
neurodegenerative disease, and nervous system disease (Table [134]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 [135]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 [136]1F, Tables [137]S3,[138]S7, Supporting Information).
Figure 2.
Figure 2
[139]Open in a new tab
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).
2.3. 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 intra‐timepoint clustering although with
slightly higher variances (Figure [140] 3A,[141]B; for analysis of
individual cell lines see Figure [142]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 [143]S8E, [144]S9, Supporting
Information). SVD resulted in five main differentiation stages based on
miRNA expression trajectories that overall resembled the gene
expression trajectories (Figure [145]3C–[146]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 [147]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 [148]S8F, Supporting Information).
Figure 3.
Figure 3
[149]Open in a new tab
Differential miRNA expression and trajectory‐based differentiation
stages. A) Principal component analysis along the sample axis of miRNA
expression during iDN development projected onto the first two
principal components. Explained variance for each principal component
(PC1, PC2) 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
iPSC‐derived 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.
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 [150]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 [151]3F,[152]G).
2.4. 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 [153] 4A). Cluster analysis revealed four
miRNA::mRNA interaction clusters associated with the nociceptor stage
(Figure [154]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
[155]S8, Supporting Information).
Figure 4.
Figure 4
[156]Open in a new tab
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.
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 [157]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 [158]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 synapse‐related pathways (mainly in the
red and yellow consensus gene modules, Figure [159]4D, Table [160]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) (Figure [161]4D,[162]E, Table [163]S9, Supporting Information). In
contrast, miRNAs of the miR‐302 family, hsa‐miR‐25‐3p, and
hsa‐miR‐363‐3p targeting genes within the yellow module, were more
implicated in the regulation of chemical synapse transmission (SCN2A,
SHTN1, and PRICKLE2, cluster 1) and postsynapse (GRIK, RBFOX1, COL12A1,
NRXN1 and SATB1, cluster 2, Figure [164]4D,[165]F).
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 [166]4D,[167]F). Enrichment analysis of all high confidence
targets revealed an important function in nervous system and
postsynapse development (Table [168]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
high‐confidence targets of hsa‐miR‐25‐3p (Figure [169]4G) highlighting
this miRNA as an interesting novel target for neuron‐related disorders.
2.5. 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,^[ [170]^48 ^] but also
to allow for nerve regeneration following peripheral nerve injury.^[
[171]^49 ^] Of the seven consensus gene modules related to neurite
outgrowth (Figure [172]1E), three (royalblue, violet, and floralwhite)
were strongly enriched for neurite outgrowth, cell migration, and
cytoskeleton regulation (Figure [173] 5A,[174]B, Table [175]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 [176]5A,[177]B, Table [178]S4, Supporting Information).
Figure 5.
Figure 5
[179]Open in a new tab
Neurite outgrowth related miRNA::gene interactions. A,B) Top 30 hub
gene networks (edges represents gene::gene correlations) based on the
ranked kME for the royalblue, violet, and floralwhite modules as well
as the averaged standardized gene expression trajectories (vst(counts))
for all three cell lines and modules highly correlates with neurite
outgrowth. C) Number of top hub genes regulated by miRNAs predominantly
targeting k‐Means cluster 1 and 2 for the royalblue and violet module;
clusters were named according the G:Profiler enrichments. D) miRNAs
with the highest number of targets for the royalblue and violet module
were determined for each k‐Means cluster and hierarchical clustering
along the z‐scored cluster axis performed. E) Chord network analysis of
hsa‐miR‐103a‐3p, hsa‐miR‐31‐5p, and hsa‐miR‐155‐3p targets with edges
representing StringDB interactions and arc colors representing the main
differentiation stage affiliation of the genes.
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 [180]5C,[181]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 high‐confidence 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 [182]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 [183]5C–[184]E, Table [185]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 [186]5E, Table
[187]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 [188]5D,[189]E), which regulate membrane trafficking.^[
[190]^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^[ [191]^51 ^] (Figure [192]5D,[193]E,
Table [194]S14, Supporting Information).
2.6. 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 all currently
assigned ion channels (Figure [195] 6A). Six consensus gene modules
were significantly associated with ion channel activity, including the
brown and cyan modules (Figure [196]6B; brown: top hub genes: ESRP1,
KCNK6; enrichments: gated channel activity, integral component of the
plasma membrane; cyan: top hub genes: [197]AL365203.2, CYB561D1,
SLCO2B1, NGFR, BIN1; enrichments: plasma membrane, ion channel
regulator activity, Table [198]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).^[ [199]^52 ^]
Figure 6.
Figure 6
[200]Open in a new tab
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.
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 [201]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 [202]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 function in mammalian neuronal tissues and in the pathogenesis
of neuropathic pain^[ [203]^53 , [204]^54 , [205]^55 ^] peaked already
at D9 at the transition from pluripotency to neuronal differentiation
(Figure [206]1F, Figure [207]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.^[ [208]^56 ^]
2.7. 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 [209]6C) and later during
differentiation (D16‐D36; Figure [210]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 [211]6D,[212]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 [213]6D,[214]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^[
[215]^57 , [216]^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 [217]6E).
2.8. 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 miRNAs 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) ([218]www.genome.jp/kegg) and
disease pathways ([219]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
([220]www.biit.cs.ut.ee/gprofiler) and StringDB PPI analysis
([221]http://www.string‐db.org/, Figure [222] 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 [223]7A) and single genes or gene sets miRNA target
spectra determined.
Figure 7.
Figure 7
[224]Open in a new tab
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.
Since a number of severe human pain disorders are strongly associated
with gain or loss‐of‐function mutations of ion channels in neurons,^[
[225]^59 , [226]^60 , [227]^61 , [228]^62 , [229]^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 [230]7B, Table [231]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 [232]7B, contingency table statistics > 0.05; Table
[233]S19/[234]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 [235]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 [236]7B, Table [237]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, [238]www.cuilab.cn) revealed two groups of miRNA which were
associated with disease entities (Figure [239]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 [240]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 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,^[ [241]^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 [242]S10A, Supporting Information). Second, WNT
signaling is crucial for corticogenesis/cortex formation, and WNT7B and
WNT8A expression increases during cortical specification.^[ [243]^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
[244]S10B,C, Supporting Information).^[ [245]^64 , [246]^65 ^] Both
complex NOCICEPTRA 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.
3. 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.^[ [247]^23 , [248]^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.
3.1. 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.^[ [249]^16 , [250]^24 , [251]^67 ^] 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.^[ [252]^3 , [253]^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 (early‐differentiation—neural progenitor
stage), SOX10 (neural progenitor stage), and CLCN3/SCN9A (nociceptor
stage), as well as gene pathway enrichments for each stage.^[ [254]^68
, [255]^69 , [256]^70 ^] Characteristic trajectories of mRNAs and
miRNAs emerged for the three relevant features of maturating
nociceptors: neurite growth, synapse, and ion channel expression.
3.2. 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.^[ [257]^71 , [258]^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.^[
[259]^73 ^] In particular, the miR‐17≈92 cluster, was highly expressed
and associated with the early differentiation of iPSC‐derived
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.^[ [260]^74 ,
[261]^75 , [262]^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.
3.3. 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,^[ [263]^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.^[ [264]^78 , [265]^79 ^] Interestingly, miR‐20b/miR‐106
directly suppress NEUROG2,^[ [266]^80 ^] thereby indirectly favoring
NEUROG1 expression and nociceptor differentiation.
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
RABs—brain‐associated Ras small G‐proteins, which are well associated
with neurite outgrowth.^[ [267]^81 , [268]^82 , [269]^83 ^] The
antiproliferative miRNAs hsa‐miR‐31‐5p but also hsa‐miR‐103a‐3p were
highly abundant early in differentiation.^[ [270]^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.^[ [271]^85 ,
[272]^86 ^]
3.4. 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 voltage‐gated 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.^[ [273]^87 , [274]^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.^[ [275]^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.^[ [276]^90 , [277]^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,^[ [278]^57 ,
[279]^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.
3.5. Differentiation of the Synaptic Release Machinery
Although native primary sensory afferents do not form autapses,^[
[280]^92 ^] iDNs express synaptic proteins.^[ [281]^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.^[ [282]^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.^[ [283]^94 , [284]^95 ,
[285]^96 , [286]^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.^[ [287]^98 ^] To further
investigate timepoint specific mRNA and miRNA regulation, we developed
the open access NOCICEPTRA tool to extract expression patterns and time
trajectories of miRNAs and their target gene space for mechanistic
studies.
3.6. 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 NOCICEPTRA 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 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 miRNAs 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.
4. 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.
5. Experimental Section
Induced Pluripotent Stem Cell (iPSC) Lines
Three different iPSC lines, adult clone 2 (AD2; SFC‐AD2‐1, male), adult
clone 3 (AD3; SFC‐AD3‐1, female), and clone 840 (SFC‐840‐03‐01,
female), derived from healthy donors with no known history of neuronal
diseases as described,^[ [288]^18 ^] were used throughout the
experiments for reprogramming, and paired long and short RNA
co‐sequencing (Table [289]S16, Supporting Information). Cells were
stored in liquid nitrogen until thawing. After thawing, cells were
passaged twice and maintained in mTeSR1 medium at 37 °C and 5% CO[2.]
iPSC‐Derived Sensory Neuron (iDN) Differentiation
AD2, AD3, and 840 iPSCs were cultured on Matrigel (Corning) coated
6‐well plates and maintained in mTeSR1 medium (STEMCELL Technologies).
Nociceptor differentiation was performed as previously published.^[
[290]^18 , [291]^19 ^] In brief, ≈80 × 10^4 cells per 6‐well plate were
seeded as single‐cells at passage 2 after thawing. Cells were
maintained overnight in mTeSR1 medium containing the Rho‐associated
kinase (ROCK) inhibitor Y‐27632 (Sigma Aldrich) until cells reached
60–70% of confluency. Subsequently, dual‐SMAD inhibition was initiated
by inhibiting BMP4 and TGF‐β in knockout serum replacement (KSR;
Knockout DMEM [Gibco], 15% KO serum replacement [Gibco], 1%
non‐essential amino acids (PAA), 2 mm Glutamax [Gibco], 1×
Penicillin/Streptomycin [Sigma Aldrich], and 100 µm β‐mercaptoethanol
[Gibco] medium supplemented with 100 nm of LDN‐1931189 (STEMGENT) and
10 µm SB431542 (selleckchem.com) for 5 days. Additionally, 3 µm
CHIR99021 (GSK3beta inhibitor, selleckchem.com), 10 µm DAPT
(gamma‐secretase inhibitor, selleckchem.com) and 10 µm SU5402 (VEGFR‐2,
FGFR‐1, and PDGFRB inhibitor, selleckchem.com) were added from day 2
(D2) to D12. From D4 on N2/B27 (Neurobasal medium/GLUTAMAX/2% B27 w/o
vitamine A [Gibco]/1% N2 [Gibco] was added to the KSR medium every
second day in increments of 25%. On D12, 80 × 10^4 cells were passaged
onto Matrigel coated 6‐well plates according to the supplier's
recommendations and maintained in N2/B27 medium supplemented with 25 ng
mL^−1 nerve growth factor (hNGF, PreproTech) 25 ng mL^−1 brain derived
neurotrophic factor (hBDNF, PreproTech) and 10 ng mL^−1 glia derived
neurotrophic factor (hGDNF, PreproTech). Y‐27632 was added for 1 day
following passaging, to promote survival of neuronal‐like progenitor
cells. Cytosine‐β‐D‐arabinofuranoside (4 µm, Sigma Aldrich) was
administered on D14 for ≈24 h to reduce the amount of
non‐differentiated proliferating cells. Cultures were maintained until
D36 and medium was changed every other day.
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 [292]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
([293]www.gencodegenes.org) using STAR (v2.7) software with the default
parameters.^[ [294]^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
(AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC) were trimmed using Flexbar v3^[
[295]^100 ^] and aligned to the human reference genome (Parameters can
be found in Table [296]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,
[297]www.genecod‐egenes.org) and the hsa.gff file provided by miRbase
(release 22.1, [298]www.mirbase.org), respectively.
Transcripts per million (TPMs) were generated from counts using the CLC
Genomics Workbench 12.0.3 (parameter: Missmatch Cost‐2,
Insertion/Deletion Cost‐3, Length‐fraction‐0.8,
Similarity‐fraction‐0.8), as a second unbiased pipeline aligning
against (GRCH38.p7, NCBI).
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 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.^[ [299]^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^[ [300]^102 ^] to assign the biotype to each gene.
Weighted Gene Co‐Expression Network Analysis
A weighted gene co‐expression network analysis (WGCNA) was used to
detect modules of genes that were highly correlated.^[ [301]^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
blockwiseConsensusModules function (softthresholding power = 16
(mRNAs)/7 (miRNAs), maxBlockSize = 37 000, network‐Type = “signed,”
minModuleSize = 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:
[MATH: aij=corxi,xjb :MATH]
(1)
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 module‐wide 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.^[ [302]^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,
[303]www.disgenet.org). g:Profiler
([304]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.^[ [305]^77 , [306]^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).
The Bioconductor R‐package DOSE was used to determine potential disease
enrichments for each supercluster/module (parameters: p‐value cutoff:
0.01, pAdjustMethod: “BH,” minGSize = 10, maxGSSize = 10 000,
qvalueCutoff = 0.2) using DisGeNET database ([307]www.disgenet.org) and
the human microRNA disease database (HMDD v3) was used to determine
miRNA specific disease enrichments,^[ [308]^106 , [309]^107 ^] using
contingency analysis provided by the python framework “statsmodel
v0.12.2” and the function contingency table. Enrichments and
contributions of each individual differentiation stage were determined
by means of standardized Pearson Residuals using the python package
“statsmodel v0.12.2” and a Pearson residual above > abs(+‐2) were
considered as putative stage specific enrichment or depletion.
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,^[ [310]^108 ^] miRTarBase,^[
[311]^109 ^] and TarBase^[ [312]^110 ^] were used to determine the
validation status, and DIANAs microT‐CDS algorithm^[ [313]^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:
[MATH:
Targetscore=p+v∗r∗rank∗100 :MATH]
(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.^[
[314]^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 ([315]https://cytoscape.org).
Hub Genes::Hub miRNAs Interactions
Genes considered hub genes were characterized by the kME of above > 0.8
and a p‐value < 0.05. Hub‐gene networks for each module were
constructed using the top 30 hub‐genes 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 miRNAs
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.
Indirect Immune Fluorescence Microscopy
AD3 (SbAd‐03‐01) and 840 (SFC840‐03‐01) iPSC clones differentiating
into iDNs were fixed on D0, D16, D26, D36 with 4% paraformaldehyde in
PBS at room temperature, for 10 min. Fixed iDNs were blocked with 5%
normal goat serum (NGS) in PBS supplemented with 0.2% BSA and 0.2%
Triton X‐100, for 30 min. Primary antibodies were applied at 4 °C
overnight in a humidified chamber, and detected by
fluorochrome‐conjugated secondary antibodies (Alexa goat anti‐rabbit
A594 (#A32740), goat anti‐mouse A594 (#A32742), goat anti‐rabbit A488
(#A32731), goat anti‐mouse A488 (#A32723), and goat anti‐rat (#A11007,
all 1:4000; Thermo Fisher Scientific). Primary antibodies used were
anti‐TUJ1 (1:600, mouse monoclonal, R&D Systems, #MAB1195);
anti‐Substance P (SP) (1:1000, mouse monoclonal Abcam #SP‐DE4‐21),
anti‐CGRP #149598 (1:500, rabbit polyclonal, Cell Signalling),
anti‐Nav1.7 (SCN9A) #ASC‐008 (1:1000,rabbit polyclonal, Alomone Labs),
anti‐PIEZO2 #APC‐090 (1:500, rabbit polyclonal, Alomone Labs),
anti‐panTRK #ab76291,(1:500, rabbit monoclonal, Abcam), anti‐NANOG
#ab62734 (1:500, mouse monoclonal, Abcam), anti‐Ki‐67,
#ab156956,(1:250, rat monoclonal, Abcam), anti‐TMEM119, #ab209064,
(1:500, rabbit monoclonal, Abcam), anti‐PERIPHERIN, #ab4666, (1:500,
rabbit polyclonal, Abcam), anti‐GFAP #04‐1062 (1:250, rabbit
monoclonal, Merck Millipore). Nuclei were counterstained with DAPI
(4′,6‐diamidino‐2‐phenylindol) 1:10000 (Thermo Fisher Scientific).
Images were recorded using an Axioimager 2 Microscope with cooled CCD
camera (SPOT Imaging Solutions) for glass coverslips and an inverted
Axiovert 200 m setup for Ibidi dishes (both Carl Zeiss Light
Microscopy). Average fluorescence intensities were quantified using
Metaview software (version 7.8.0.0, Molecular Devices, LLC, San Jose,
CA95134) with the line scan plug‐in to quantify fluorescent intensities
along a defined line cutting individual neuronal structures. Cell
counts positive and negative for primary antibodies were correlated to
DAPI positive nuclei and compared between the individual
differentiation days. In general 2 dishes/cellline/condition (Day 16‐D
36) were stained and analyzed in at least 2 individual
RT‐qPCR Analysis
Reverse transcription and qPCR reactions were performed on the same
samples that were used for sequencing according to the protocol
provided by the supplier (Thermofisher scientific). In short, each
reverse transcription reaction contained 10 ng of total RNA, 1X reverse
transcription buffer, 5.5 mm MgCl2 (GeneAmp 10X PCR Buffer II &
MgCl[2], #N8080130), 1 mm dNTPs, RNase inhibitor (#N8080119), 50 units
of MultiScribe Reverse Transcriptase (#4311235) and 1X RT specific
primers (Table [316] 2 ) volume of adjusted to 15 µL with nuclease free
water (#R0582). The RT program was set to 30 min at 16 °C, 30 min at
42 °C, 5 min at 85 °C, followed by a holding step at 4 °C. After
reverse transcription, reactions were stored overnight at −20 °C. Each
qPCR reaction contained 1.33 µL of the RT product, 1X TaqMan Universal
Master Mix II, no UNG (#44440049), 1X of the appropriate assay
(Table [317]2) and nuclease free water up to a final volume of 20 µL.
Reactions for each sample were prepared as technical duplicates,
alongside reverse transcription non‐template controls, loaded on
MicroAmp Fast Optical 96‐well reaction plates (#4346906), and placed in
the QuantStudio 6 Pro Real‐Time PCR system (all reagents for RT‐qPCR
were purchased from ThermoFisher Scientific). The PCR cycle protocol
was: 10 min at 95 °C, 40 two‐step cycles of 15 s at 95 °C and 1 min at
60 °C. In order to account for potential methodological variabilities
all RT (day 1) and qPCR reactions (day 2) were processed in parallel.
Threshold was set at 0.1 and results were obtained from the QuantStudio
Design and Analysis software. The expression levels of the miRNA of
interest were quantified and expressed as relative expression to the
respective expression of the reference gene (U6) using the 2^−∆Cq
method.
Table 2.
miRNA TaqMan Assay IDs
miRNA Assay ID
hsa‐miR‐25‐3p 000403
hsa‐miR‐103a‐3p 000439
hsa‐miR‐146a‐5p 000468
hsa‐miR‐302a‐5p 002381
hsa‐miR‐335‐5p 000546
Reference gene
U6 001973
[318]Open in a new tab
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 [319]S21, Supporting Information.
Conflict of Interest
The authors declare no conflict of interest.
Author Contributions
M.Z., K.K., C.L.S., T.K., and G.K. performed the experiments; M.Z. and
K.K. developed the NOCICEPTRA tool; M.Z., K.K., and M.K. conceptualized
the study; M.Z., K.K., M.Z.C., and M.K. wrote the manuscript.
Supporting information
Supporting Information
[320]Click here for additional data file.^ (5.2MB, pdf)
Supplemental Tables
[321]Click here for additional data file.^ (223.1KB, xlsx)
Acknowledgements