In total, 16 male, 8–12 weeks old C57Bl/6J mice (Janvier Labs, France) were used for the experiments. Animals were kept under standard pathogen free (SPF) conditions, at 24°C on a 12 h light/dark cycle and had free access to autoclaved pelleted food and water. All mice were treated in accordance with the Ethics Guidelines of Animal Care (Medical University of Innsbruck), the European Communities Council Directive of 22nd September 2010 on the protection of animals used for scientific purposes (2010/63/EU) and all procedures were approved by the Austrian National Animal Experiment Ethics Committee of the Austrian Bundesministerium für Wissenschaft und Forschung (permit number BMWF-66.011/0054-WF/V/3b/2015).

Mice were divided in three groups, non-treated (n = 5), sham (n = 6), and SNI (n = 5). The SNI model was adopted from Decosterd and Woolf (2000). Briefly, mice were anesthetized with a mixture of xylazine (10 mg/kg, AniMedica, Germany) and ketamine (100 mg/kg, Graeub, Switzerland). The skin of the lateral surface of the left thigh was incised and the sciatic nerve was exposed. For the SNI procedure, the common peroneal and the tibial nerves were ligated with 4-0 vicryl (Sh-1 plus, Ethicon, Austria) and a portion of approximately 1–2 mm length was excised, leaving the sural nerve intact. Mice subjected to sham surgery had their sciatic nerve exposed but not lesioned. After surgery, the skin was sutured using 4-0 vicryl and mice were placed on a heat block adjusted to 37°C until recovery. Non-treated mice were not subjected to any treatment.

Non-treated, sham and SNI mice 7 days post-surgery (only for sham and SNI mice) were deeply anesthetized with isoflurane (Forane, Abbott, United Kingdom) and euthanized by decapitation. Lumbar ipsilateral L3-L5 DRG and dhSC as well as contralateral mPFC were microdissected, snap-frozen in liquid nitrogen, and stored at –80°C until use.

RNA extraction was performed using peqGOLD TriFast reagent (Peqlab Biotechnologie, Germany), according to the protocol provided by the manufacturer [chloroform (C2432) and absolute ethanol (107017) were obtained from Merck, United States]. The air-dried RNA pellets were diluted in nuclease free water (R0582, ThermoFisher Scientific, United States). RNA quantity was estimated using Nanodrop 2000 (ThermoFisher Scientific) and RNA integrity [RNA integrity number (RIN)] was assessed using 2100 Bioanalyzer (Agilent Technologies, United States) in the deep-sequencing core facility of the Medical University of Innsbruck.

We evaluated the stability of all ncRNAs, suggested as controls and available by Thermo Fisher Scientific. The expression of small nucleolar (snoRNA) and small nuclear (snRNA) RNAs, as potential reference genes, was quantified using Taqman MicroRNA Control Assays, which through a two-step protocol ensure high specificity and sensitivity and do not require prior DNase treatment of the RNA template, since the stem-loop transcription amplifies only the mature sequence (Mestdagh et al., 2008). Reverse transcription and qPCR reactions were prepared according to the protocol provided by the supplier. Briefly, each reverse transcription reaction contained 10 ng of total RNA, 1X reverse transcription buffer, 5.5 mM MgCl₂ (GeneAmp^® 10X PCR Buffer II and MgCl₂, #N8080130), 1 mM dNTPs, RNase inhibitor (#N8080119), 50 units of MultiScribe^TM Reverse Transcriptase (#4311235), and 1X RT specific primers (Table 1), in a final volume of 15 μL, adjusted with nuclease free water (#R0582). The RT program was: 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 at –20°C until next day. 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 1), 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 7500 Fast RT-PCR system (all reagents for RT-qPCR were 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) per tissue assessed were processed on the same day. The MIQE checklist for authors, reviewers, and editors can be found in Supplementary Table S1.

Results were extracted from the 7500 Software v2.3 (ThermoFisher Scientific), by manually setting the threshold at 0.1 (single threshold method) for all amplicons and keeping the automatic baseline (baseline start cycle: 3, baseline end cycle: 15). Additionally, raw amplification data were imported to LinRegPCR program (Ramakers et al., 2003; Ruijter et al., 2009, 2015; Tuomi et al., 2010), which performs a baseline correction of the amplification data, determines a window of linearity, and through a linear regression analysis determines the qPCR efficiency and the Cq value per reaction, in order to estimate the mean PCR efficiency per primer set or assay [52]. The efficiencies for each amplicon were calculated using LinRegPCR (v2014) for hydrolysis probes with the following parameters: window of linearity: four points; exclude no plateau samples; include efficiency outlier samples; log-linear phase criterion: strictly continuous log-linear phase. Mean Cq values for each sample and ncRNA were calculated as the average of the technical duplicates. The stability of all evaluated ncRNAs was estimated based on the Cq values obtained from single threshold settings and LinRegPCR program adjustments, by using four different statistical approaches that are commonly used for stability assessment of reference genes: BestKeeper (Pfaffl et al., 2004), the comparative delta-Cq method (Silver et al., 2006), geNorm (Vandesompele et al., 2002), and NormFinder (Andersen et al., 2004).

For BestKeeper (Pfaffl et al., 2004), the appropriate reference genes need to exhibit low variability and comparable levels of expression. Variability was defined by the standard deviation (SD) of mean Cq values and low SD values corresponded to stable ncRNAs. Candidates with SD > 1 were considered inappropriate for normalization. Similar levels of expression were explored using Pearson’s linear correlation, by comparing each ncRNA with the BestKeeper index (BKI), which was the geometric mean of the mean Cq values of all evaluated ncRNAs. The closer the correlation was to 1, the more stably an ncRNA was considered to be expressed. A mean ranking score for each ncRNA was calculated by averaging the scores obtained by SD and correlation with the BKI rankings.

In the comparative delta-Cq method (Silver et al., 2006), all pairs of ncRNAs were compared with each other, according to their Cq differences (ΔCq). Variability was determined by calculating the average SD between the pairwise comparisons of all ncRNAs and a low average SD indicated increased stability.

GeNorm (Vandesompele et al., 2002) calculated the average expression stability M, which represents the average pairwise variation of each evaluated ncRNA compared to all the rest. Low M values suggested that the ncRNAs expression was stable. Subsequently, geNorm performed repeated stepwise exclusions (by repeatedly omitting the candidate with the worst M value) and identified the most stable pair of ncRNAs. Additionally, geNorm identified the optimal number of reference genes, by calculating the pairwise variation coefficient (V value).

NormFinder (Andersen et al., 2004), besides calculating the variation of expression, estimated the variation between groups. Specifically, for each ncRNA the stability value ρ was estimated as well as the intra- and inter-group variations. The lowest ρ values represented more stably expressed candidates. Furthermore, a stably expressed ncRNA needed to exhibit intra- and inter-group variations as close as possible to zero.

Each ncRNA was subsequently ranked according to its performance in each of the four different methods (with 1 indicating the most stable ncRNA). The overall stability of each candidate was estimated by calculating the mean of the rankings provided by each statistical approach.

The expression levels of miR-21a-5p were quantified in the DRG of non-treated, sham, and SNI mice and expressed as fold changes relative to the respective expression in the non-treated mice using the 2^–ΔΔCq method. ncRNAs were assessed for their suitability for the quantification of miR-21a-5p in DRG by using as a reference gene(s): (1) three (sno420, sno429, and sno202) and two (sno420 and sno429) most stable ncRNAs from the overall ranking, (2) two most stable pair of ncRNAs as suggested by geNorm (sno202 and sno420) and Normfinder (sno234 and sno429), and (3) each ncRNA as a single normalizer.

For statistical data analyses, GraphPad Prism 7.00 and IBM SPSS Statistics 21 were used. Statistical tests used were specified in the text or in their respective figures or tables. Statistical analysis for miR-21a-5p expression was performed using one-way analysis of variance (ANOVA) with the group of animals (non-treated, sham, and SNI) as the independent factor, followed by Tukey’s multiple comparison test, when appropriate, in order to identify differences between groups. The level of statistical significance was predefined at p < 0.05. Graphs were plotted in GraphPad Prism 7.00 and illustrated in CorelDRAW X7.

RNA quantity and integrity varied in between samples and tissues (Supplementary Table S2). For DRG the median total RNA concentration was 102.3 ng/μL with an interquartile range (IQR) of 88.85–122.6 ng/μL, for dhSC the median RNA concentration was 127.8 ng/μL with an IQR of: 106.2–176.2 ng/μL, and for mPFC the median RNA concentration was 68.35 ng/μL with an IQR of: 58.43–88.95 ng/μL. Analysis of integrity provided median values of RIN = 6.85 (IQR: 6.45–7.55) in DRG, RIN = 7.55 (IQR: 6.825–8.0) in dhSC, and RIN = 7.35 (IQR: 7.1–7.75) in mPFC.

Mean qPCR efficiencies for all evaluated ncRNAs were calculated by the LinRegPCR program (Supplementary Table S3). R² values for all amplicons, as calculated by LinRegPCR, were > 0.998. Correlation coefficients for RIN values and amplicons’ efficiencies were not statistically significant (Supplementary Table S4). Furthermore, no statistically significant difference was observed for individual assay performance on different experimental days (Supplementary Table S5A) or in different experimental groups (Supplementary Table S5B), suggesting that all assays performed similarly in all experimental groups, experimental days, and tissues.

After setting the common threshold at 0.1 and after importing the raw amplification data into LinRegPCR, we assessed the raw Cq values (Figures 1A–C). All ncRNAs were successfully amplified and both settings provided similar results in DRG, dhSC, and mPFC. In all three tissues, the most abundantly expressed ncRNA was sno202, followed by U6, whereas the least abundant ncRNA was sno412 (Figures 1A–C and Supplementary Table S6). We further analyzed the SDs of the technical duplicates per sample and amplicon, since reproducibility is critically important for RT-qPCR quantification (Figures 1D–F and Supplementary Table S6). SDs for technical replicates ≤0.167 are required to detect a 2-fold change of expression in 99.7% of cases (Applied Biosystems, 2016). Since sno55 exhibited SDs > 0.167 in the technical duplicates in all three tissues, it was excluded from further analysis. Likewise, sno412 was excluded from assessment in the mPFC. In DRG and dhSC, U6 also showed SDs > 0.167 according to LinRegPCR settings, however, was kept for further analysis, since the common threshold setting suggested appropriate variability in technical replicates.

Data extracted using both the common threshold and LinRegPCR adjustments, provided similar results in BestKeeper evaluation (Figure 2 and Supplementary Table S7). In DRG (Figure 2A and Supplementary Table S7A), U6 showed the lowest SD (SD = 0.336 for the common threshold and 0.341 for LinRegPCR) across non-treated, sham, and SNI mice, whereas sno251 showed the highest SD (SD = 0.602 for the common threshold and 0.604 for LinRegPCR). All SDs were below the proposed cut-off value of < 1 for the exclusion of a candidate. Correlation analysis between each ncRNA and BKI revealed r > 0.92 for sno202, sno420, and sno429, whereas sno135, sno292, and U6 had r values <0.67. In the dhSC (Figure 2B and Supplementary Table S7B), sno429 exhibited the lowest SD (0.138 and 0.113, according to the common threshold and LinRegPCR, respectively) and sno234 the highest SD (0.333 using the common threshold and 0.325 according to LinRegPCR). r values for most ncRNAs were surprisingly low and p values were not significant for sno135, sno142, sno251, sno412, and sno429, indicating no correlation with the BKI. Out of the ncRNAs, which significantly correlated with the BKI, the best correlation was determined for sno292 (r = 0.738), followed by U6 (r = 0.732) for the common threshold, whereas for LinRegPCR processed data, the best correlation was observed for U6 (r = 0.777), followed by sno292 (r = 0.735). In the mPFC (Figure 2C and Supplementary Table S7C), sno202 showed the lowest SD (0.494 for the common threshold and 0.498 for LinRegPCR) and sno135 the highest SD (0.900 for the common threshold and 0.905 for LinRegPCR). Furthermore, in the data normalized with the common threshold, sno202 exhibited the highest correlation (r = 0.918) with the BKI, whereas in the LinRegPCR processed data, the best correlation was observed for sno429 (r = 0.923).

The comparative delta-Cq method provided similar results for both common threshold and LinRegPCR processed data (see Figure 3 and Supplementary Table S8 for details). In DRG, the most stable ncRNA was sno429 (with a ΔCq SD = 0.398 in the common threshold and 0.388 in LinRegPCR), whereas sno135 was the least stable, with ΔCq SD > 0.66 (Figure 3A and Supplementary Table S8A). In dhSC, sno429 showed the lowest ΔCq SD (0.322 in the common threshold and 0.332 in LinRegPCR) and sno234 the highest (0.529 in the common threshold and 0.524 in LinRegPCR; Figure 3B and Supplementary Table S8B). Overall, in the mPFC all ncRNAs had higher ΔCq SDs (Figure 3C and Supplementary Table S8C). Among them, sno202 displayed the highest stability (with a ΔCq SD = 0.521 in the common threshold and equal to 0.522 in LinRegPCR processed data) and sno135 the lowest (ΔCq SD = 0.947 in the common threshold and 0.956 in LinRegPCR data).

According to the developers’ recommendations, the M values should not exceed 1.5 and all ncRNAs in all tissues analyzed with both common threshold and LinRegPCR settings met this requirement. Initial M values (with all ncRNAs included in the analysis), as computed by geNorm are reported in Supplementary Table S9. After step-wise exclusion of the worst ncRNA and pairwise variation for the determination of the optimal number of ncRNAs for normalization, the average expression stability of the remaining candidates was determined (Figure 4). In DRG, sno135 was identified as the least stable ncRNA; whereas the sno202/sno420 pair was the most stable (Figure 4A). In the dhSC, sno234 showed the lowest and the sno251/sno429 pair the highest stability (Figure 4B). In the mPFC, sno135 was determined as the least stable ncRNA and the sno142/sno202 pair as the most stable (Figure 4C). Pairwise variation analysis showed that in all cases the V values were below 0.15, suggesting that the best pair of ncRNAs would already be sufficient for appropriate normalization (Figures 4D–F).

Stability values ρ, inter- and intra-group variability for all ncRNAs, were evaluated (Figure 5; detailed data are provided in Supplementary Table S10). Both threshold settings provided similar ranking results: In DRG, the lowest ρ value was observed for sno420 when using the common threshold (ρ = 0.121) and sno429 when using LinRegPCR (ρ = 0.090; Figures 5A,B). The highest ρ value was observed for sno292 (ρ = 0.306 for the common threshold and ρ = 0.265 with LinRegPCR). The sno234/sno429 pair was identified as the best pair combination (ρ = 0.069 and ρ = 0.048, respectively for the common threshold and LinRegPCR settings; Table 2). In the dhSC, the most stable ncRNA was sno429 (ρ = 0.071 using the common threshold and ρ = 0.077 using LinRegPCR; Figures 5C,D) and the best pair combination sno429/U6 (ρ = 0.062 for both settings; Table 2), whereas sno234 exhibited the worst performance (ρ = 0.158 for both settings). In the mPFC, the setting of a common threshold revealed that sno142 was the most stable ncRNA (ρ = 0.113), whereas sno420 showed the best performance when using LinRegPCR (ρ = 0.088; Figures 5E,F). The most stable pair was sno202/sno420 (ρ = 0.078) for the common threshold and sno142/sno420 (ρ = 0.070) for LinRegPCR (Table 2). Both the common threshold and LinRegPCR settings identified sno135 as the least stable ncRNA (ρ = 0.374 for the common threshold and ρ = 0.317 for LinRegPCR).

Each potential reference gene was assigned a rank, according to its performance (1 was assigned to the most stable ncRNA) in each analysis. For BestKeeper analysis, ncRNAs were ranked based on both their SD and correlation with the BKI. Finally, we calculated the arithmetic mean and compiled an overall ranking list (Figure 6 and Table 3). Since all four methods employed for assessing the stability of the ncRNAs, and their potential as reference genes, were taken into account as of equal importance, we did not use the geometric mean as it indicates the tendency of a set of values and dampens the effect of extreme values. Nevertheless, the rating was also calculated using the geometric mean and only sno251 in the dhSC for the common threshold exhibited a slightly different ranking (one rank higher), compared to the ranking obtained by using the arithmetic mean (Supplementary Table S11). Both the common threshold and LinRegPCR settings revealed the same three ncRNAs as the most stable and the same two candidates as the least stable, which, however, were different between the three neuronal tissues. Specifically, in DRGs, sno420, sno429, and sno202 ranked as the top three stable ncRNAs, whereas the lowest stability was observed for sno135, sno292, and sno251 (Figure 6A). In the dhSC, sno429, sno202, and U6 had the best stability, whereas sno234 and sno135 were the most variably expressed ncRNAs (Figure 6B). In the mPFC, the best three candidates were sno202, sno420, and sno142, whereas the least stable ncRNAs were sno135, sno234, and U6 (Figure 6C). The overall ranking for DRG and mPFC resembled the ranking provided by geNorm analysis. Furthermore, the overall ranking did not substantially differ between the common threshold and LinRegPCR datasets.

The ranked ncRNAs were assessed for their suitability to quantify miR-21a-5p, which is upregulated in DRG after peripheral nerve injury (Strickland et al., 2011; Zhou et al., 2015; Hori et al., 2016; Chang et al., 2017; Karl et al., 2017; Simeoli et al., 2017; Zhang Z. J. et al., 2018). miR-21a-5p was found upregulated in the SNI group when three (sno420/sno429/sno202; F₂,₉ = 10.72, p = 0.0042; non-treated vs. SNI p = 0.0058, sham vs. SNI p = 0.0110; Figure 7A, left panel) or two (sno420/sno429; F₂,₉ = 12.05, p = 0.0028; non-treated vs. SNI p = 0.0035, sham vs. SNI p = 0.0097; Figure 7A, right panel) most stable ncRNAs from the overall ranking list were used as reference genes. The most stable pair of genes as computed by geNorm (sno202/sno420; F₂,₉ = 9.781, p = 0.0055; non-treated vs. SNI p = 0.0086, sham vs. SNI p = 0.0121; Figure 7B) and NormFinder (sno234/sno429; F₂,₉ = 24.34, p = 0.0002; non-treated vs. SNI p = 0.0003, sham vs. SNI p = 0.0009; Figure 7C) provided similar results. However, when U6 was used as a normalizer, miR-21a-5p expression was different only between non-treated vs. SNI (F₂,₉ = 5.72, p = 0.0249; non-treated vs. SNI p = 0.0321; Figure 7D), whereas the use of sno135 completely abolished the significance (F₂,₉ = 1.048, p = 0.3897; Figure 7E). The use of single reference genes revealed that the five most stable ncRNAs, namely, sno420 (F₂,₉ = 11.06, p = 0.0038; non-treated vs. SNI p = 0.0052, sham vs. SNI p = 0.0104; Supplementary Figure S1A), sno429 (F₂,₉ = 11.21, p = 0.0036; non-treated vs. SNI p = 0.0039, sham vs. SNI p = 0.0153; Supplementary Figure S1B), sno202 (F₂,₉ = 8.319, p = 0.0090; non-treated vs. SNI p = 0.0160, sham vs. SNI p = 0.0158; Supplementary Figure S1C), sno234 (F₂,₉ = 33.09, p < 0.0001; non-treated vs. SNI p = 0.0001, sham vs. SNI p = 0.0002; Supplementary Figure S1D), and sno412 (F₂,₉ = 17.07, p = 0.0009; non-treated vs. SNI p = 0.0023, sham vs. SNI p = 0.0014; Supplementary Figure S1E) could also detect the upregulation of miR-21a-5p in the SNI mice in comparison to non-treated and sham mice. Normalization to sno142 showed a significant difference only between non-treated vs. SNI mice (F₂,₉ = 8.652, p = 0.0080; non-treated vs. SNI p = 0.0071; Supplementary Figure S1F), whereas sno251 (F₂,₉ = 16.63, p = 0.0010; non-treated vs. SNI p = 0.0010, sham vs. SNI p = 0.0058; Supplementary Figure S1G) and sno292 (F₂,₉ = 48.52, p < 0.0001; non-treated vs. SNI and sham vs. SNI p < 0.0001; Supplementary Figure S1H) exaggerated the upregulation of miR-21a-5p, indicating that these two ncRNAs likely were regulated by the treatment.