Integrative miRNA-mRNA profiling of human epidermis: unique signature of SCN9A painful neuropathy

Abstract

Personalized management of neuropathic pain is an unmet clinical need due to heterogeneity of the underlying aetiologies, incompletely understood pathophysiological mechanisms and limited efficacy of existing treatments. Recent studies on microRNA in pain preclinical models have begun to yield insights into pain-related mechanisms, identifying nociception-related species differences and pinpointing potential drug candidates. With the aim of bridging the translational gap towards the clinic, we generated a human pain-related integrative miRNA and mRNA molecular profile of the epidermis, the tissue hosting small nerve fibres, in a deeply phenotyped cohort of patients with sodium channel-related painful neuropathy not responding to currently available therapies. We identified four miRNAs strongly discriminating patients from healthy individuals, confirming their effect on differentially expressed gene targets driving peripheral sensory transduction, transmission, modulation and post-transcriptional modifications, with strong effects on gene targets including NEDD4. We identified a complex epidermal miRNA-mRNA network based on tissue-specific experimental data suggesting a cross-talk between epidermal cells and axons in neuropathy pain. Using immunofluorescence assay and confocal microscopy, we observed that Nav1.7 signal intensity in keratinocytes strongly inversely correlated with NEDD4 expression that was downregulated by miR-30 family, suggesting post-transcriptional fine tuning of pain-related protein expression. Our targeted molecular profiling advances the understanding of specific neuropathic pain fine signatures and may accelerate process towards personalized medicine in patients with neuropathic pain.

Keywords: SCN9A; miRNA; neuropathic pain; skin biopsy; small fibre neuropathy.

Figure 1

Skin biopsy profiling. Workflow of integrative miRNA and mRNA profiling that allowed miRNA–mRNA pathway network construction. (A) The epidermis was dissected from skin biopsy samples and total RNA was isolated. (B) The miRNA profiling was performed using Taqman Array Human cards containing 754 miRNAs. (C) The bioinformatic pipeline allowed the identification of four significantly downregulated miRNAs. (D) Next, from the same RNA sample, mRNA profiling was performed using custom-designed microfluidic TaqMan™ Array Card with 93 pain-related genes. (E) This allowed identification of significantly dysregulated genes in the epidermis of 11 NavNP patients compared to 7 HC. (F) In silico analysis using publicly available databases (TargetScan, miRTarBase and Tarbase) was used to identify miR-30 family, miR-181a-2-3p and miR-203a-3p putative gene targets matched with differentially expressed genes (DEGs). (G) Six DEGs were found to be miRNA-target genes, (H) while six other DEGs did not emerge as direct gene targets of miR-30 family, miR-181a-2-3p and miR-203a-3p. (I) KEGG and REACTOME pathways and biological processes were employed for enrichment analysis to identify pain-related epidermal terms. (J) Correlation analysis was used to identify functional relationships between miRNA and target mRNA expression, and was applied to each miRNA and its putative target that was significantly dysregulated in NavNP patients. (K) The functional network was reconstructed integrating the obtained data. The workflow scheme was designed using BioRender.com.

Figure 2

Downregulated miRNA in epidermis of painful SCN9A-related neuropathy patients (NavNP). Microfluidic analysis of miRNA profiling in total RNA extracted from the epidermis of 11 NavNP patients and 7 healthy controls (HC) demonstrated a significant reduction of miR-30d-5p (P-value 3.23 × 10⁻⁴, FC −5.83), miR-30a-5p [P-value 4.40 × 10⁻⁴, fold change (FC) −4.95], miR-203a-3p (P-value 4.40 × 10⁻⁴, FC −3.64), miR-181a-2-3p (P-value 4.40 × 10⁻⁴, FC −2.21) expression in NavNP patients compared to HC. Bar graph indicates the 2^−(ΔΔCq). The comparisons are made applying Wilcoxon rank sum test and corrected for Bonferroni multiple test. The fold change in expression was calculated as 2^−(ΔΔCq). We provide the fold change reduction in expression in the NavNP group compared to HC applying the negative inverse of 2^−(ΔΔCq).

Figure 3

miRNA target prediction and gene expression validation. (A) Unsupervised heat map of the ΔCq values, the difference between the tested gene and geometric mean as normalization factor, of significantly dysregulated miRNA-targets in NavNP patients (black) compared to HC (grey), highlighting the two separated phenotypic classes. Heat map colours correspond to mRNA expression as indicated in the colour key: red (lower expression) and blue (higher expression). (B) Scatter correlation plots between miRNAs and their gene targets. Pearson coefficient and P-value are shown in the graphs. NavNP patients are coloured in red and HC in black showing a different expression distribution between phenotypic groups.

Figure 4

miRNA–mRNA pathway functional network. Experimental miRNA and mRNA profiling data together with in silico terms enrichment are used to reconstruct functional biological network in the epidermis of NavNP patients. Four downregulated miRNAs (green hexagon) are paired with their gene targets (white ellipse with purple edge). Significant correlations between miRNA–mRNA pairs are represented with red and light blue dotted lines. All genes were then mapped within pain-related terms, grouped into four main classes, defined according to the scientific literature, as follows: (i) transduction (yellow rectangle); (ii) transmission (light orange rectangle); (iii) modulation (light blue rectangle); and (iv) post-transcriptional modification (cyan rectangle). Dysregulated genes that are not miRNA targets are represented with white ellipse with grey edge and are paired with pain-related pathways. The network represents the complex interactions occurring in the epidermis, related to pain signalling.

Figure 5

miR-30 family regulates Nav1.7 signalling in keratinocytes. (A) Representative confocal microscope image of epidermis. Nav1.7 (green) in keratinocytes of NavNP patients and HC. Scale bar = 20 µm. (B) Boxplot of mean Nav1.7 immunofluorescence intensity for two studied groups, NavNP patients and HC, respectively. NavNP patients show significantly higher Nav1.7 signal intensity than HC (P = 8.798 × 10⁻⁵). ***P < 0.001 according to Wilcoxon rank sum test. (C) Scatter correlation plot between mean Nav1.7 immunofluorescence intensity and NEDD4 expression values in NavNP patients and HC samples. (D and E) Scatter correlation plot between mean Nav1.7 immunofluorescence intensity and (D) miR-30a-5p and miR-30d-5p (E) expression values in NavNP patients and HC samples. Spearman coefficient and P-value are shown in the graph.