Reverse transcription quantitative real-time polymerase chain reaction reference genes in the spared nerve injury model of neuropathic pain: validation and literature search.

BACKGROUND: The reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) is a widely used, highly sensitive laboratory technique to rapidly and easily detect, identify and quantify gene expression. Reliable RT-qPCR data necessitates accurate normalization with validated control genes (reference genes) whose expression is constant in all studied conditions. This stability has to be demonstrated.We performed a literature search for studies using quantitative or semi-quantitative PCR in the rat spared nerve injury (SNI) model of neuropathic pain to verify whether any reference genes had previously been validated. We then analyzed the stability over time of 7 commonly used reference genes in the nervous system - specifically in the spinal cord dorsal horn and the dorsal root ganglion (DRG). These were: Actin beta (Actb), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), ribosomal proteins 18S (18S), L13a (RPL13a) and L29 (RPL29), hypoxanthine phosphoribosyltransferase 1 (HPRT1) and hydroxymethylbilane synthase (HMBS). We compared the candidate genes and established a stability ranking using the geNorm algorithm. Finally, we assessed the number of reference genes necessary for accurate normalization in this neuropathic pain model.
RESULTS: We found GAPDH, HMBS, Actb, HPRT1 and 18S cited as reference genes in literature on studies using the SNI model. Only HPRT1 and 18S had been once previously demonstrated as stable in RT-qPCR arrays. All the genes tested in this study, using the geNorm algorithm, presented gene stability values (M-value) acceptable enough for them to qualify as potential reference genes in both DRG and spinal cord. Using the coefficient of variation, 18S failed the 50% cut-off with a value of 61% in the DRG. The two most stable genes in the dorsal horn were RPL29 and RPL13a; in the DRG they were HPRT1 and Actb. Using a 0.15 cut-off for pairwise variations we found that any pair of stable reference gene was sufficient for the normalization process.
CONCLUSIONS: In the rat SNI model, we validated and ranked Actb, RPL29, RPL13a, HMBS, GAPDH, HPRT1 and 18S as good reference genes in the spinal cord. In the DRG, 18S did not fulfill stability criteria. The combination of any two stable reference genes was sufficient to provide an accurate normalization.

Background

While 20%-30% of the general populace are affected by chronic pain, a significant proportion, about 7%, experiences neuropathic pain characteristics [1,2]. Defined by the International Association for the Study of Pain (IASP) as a pain “caused by a lesion or disease of the somatosensory nervous system”, neuropathic pain represents a daily challenge to medical practice and a worldwide public health problem since it is often refractory to treatment. Several animal models of neuropathic pain are used to study its different pathophysiological mechanisms and hopefully uncover new targets for treatment. Among them, the spared nerve injury (SNI) model [3] represents an easy-to-perform and robust model of peripheral nerve injury. The axotomy of the tibial and peroneal branches of the sciatic nerve, sparing the sural nerve, allows a straightforward testing of the hind paw’s sural nerve territory for signs of hyperalgesia. After a peripheral nerve injury in rodents, numerous alterations in the somatosensory nervous system can be observed, which together lead to a hyperexcitability of the system [4]. Besides changes affecting the primary nociceptors whose cell bodies are in the dorsal root ganglion (DRG) or the secondary neurons in the dorsal horn projecting into the brain, glial cells also react dramatically to a peripheral nerve injury like SNI. A common approach used to elucidate such changes in molecular machinery is to explore modifications in the expression of relevant genes [5,6]. These variations can be detected and quantified in a sensitive, specific way using a reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR), assuming that an accurate normalization has been performed with reference genes that have proved stable in all biological replicates and experimental conditions [7-10]. The recently published guidelines on the “minimum information for publication of quantitative real-time PCR experiments” (MIQE) strongly recommended the validation of reference genes used in different cell types and conditions [7]. These guidelines outline the importance of RT-qPCR accuracy and recommend following detailed procedures to produce the most robust and reproducible RT-qPCR data.

Algorithms like geNorm or NormFinder were developed to help researchers assess candidate reference genes. GeNorm has been well described [9] and widely accepted as a very useful tool for the normalization of RT-qPCR data [11]. It measures the stability of potential reference genes by comparing their expression levels against one another. The expression ratio is supposed to be the same in all biological and technical replicates. This pairwise comparison permits the assessment and ranking of candidate reference genes and additionally demonstrates how many are necessary for an accurate normalization [9,11]. To our knowledge, a validation of reference genes for the rat after spared nerve injury had never been published. We therefore performed a literature search for all studies using quantitative or semi-quantitative PCR in the SNI model and listed all the reference genes used. We then tested a total of 7 genes – including commonly used reference genes in the nervous system and those used in the SNI model – at different time points after nerve lesion in the spinal cord dorsal horn and in the DRG. These genes were Actin beta (Actb), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), ribosomal proteins 18S (18S), L13a (RPL13a) and L29 (RPL29), hypoxanthine phosphoribosyltransferase 1 (HPRT1) and hydroxymethylbilane synthase (HMBS) [12-15]. Selected genes were related to various metabolic pathways to avoid any co-regulation [9].

Methods

Animals

We carried out our experiments on Sprague–Dawley rats (Charles River, L’Abresle, France) weighing 250–300 grams. They were housed under 12 h day/night photoperiodic conditions, at a constant temperature and with free access to water and food. All procedures were approved by the Canton of Vaud’s Animal Experimentation Committee and were in accordance with the Swiss Federal Law on Animal Welfare and the IASP’s guidelines [16].

Surgery

Rats were anesthetized with isoflurane 1.5%-2.5% (Abott, Baar, Switzerland) and surgical procedures were performed as previously described [3,17]. Briefly, the left sciatic nerve was exposed at the mid-thigh level distal to the trifurcation, tibial and common peroneal branches were tightly ligated with 5.0 silk and axotomized, leaving the sural branch intact. Muscle and skin were closed in two layers and animals were allowed to recover. Three rats were sacrificed at 2, 4, 7, 10, 14 and 21 days after surgery.

Literature search

We performed a literature search in the PubMed database using the keywords “spared nerve injury” and selected every study using semi-quantitative PCR or RT-qPCR in the rat SNI model (13 April 2013). We only selected articles in English. Reference genes used were noted, as were indications of their validation in the models.

RNA isolation and reverse transcription

Animals were transcardially perfused with NaCl 0.9% in order to expel the blood from collected tissue. Ipsilateral L4 and L5 DRG and spinal cord dorsal horn were rapidly dissected and immediately stored in RNAlater® (Qiagen, Switzerland) for stabilization at 4°C overnight, and then frozen at −80°C until processed. Samples were first homogenized with a POLYTRON® homogenizer and then total RNA was extracted using the RNeasy Plus Mini Kit (Qiagen). To ensure complete removal of gDNA, we treated the dorsal horn samples with a DNase (RNAqueous®-4PCR Kit, Ambion) and then the inactivation reagent recommended in the manufacturer’s instructions.

Nucleic acid purity was assessed measuring the A260/A280 ratio by spectrophotometer (NanoDrop). Total RNA integrity (RIN) and quantity were determined with an Agilent 2100 Bioanalyser (Agilent Technologies). RNA used had to have both a A260/A280 ratio between 1.8 and 2.1, as well as a RIN ≥8 for dorsal horn RNA extractions and ≥6.9 for DRG samples. The reverse transcription was achieved using Omniscript RT Kit (Qiagen), according to manufacturer’s instructions. All reactions took place for 1h at 37°C in a final volume of 40 μl containing 1 μg total RNA in the presence of 20 units RNase inhibitor (RNasin® Ribonuclease Inhibitor, Promega), 1 μg random hexamers (Microsynth), 4 μl 10× Buffer RT, 0.5 mM dNTPs and 8 units Omniscript reverse transcriptase (Omniscript RT Kit, Qiagen).

Selection of reference genes and primer design

Reference genes were chosen from those used previously in the SNI model. Other potentially suitable reference genes were selected among those used in published literature on the nervous system (Table  1). Whenever possible, primers fulfilled the following recommended criteria: amplicon length of 60 bp-150 bp, location of primers on two different exons, primer sequence length of 18 bp-25 bp, melting temperature of 60°C +/−1°C and GC content of 40%-60% [7,18-20]. Primer specificity was checked in silico (Primer-BLAST Tool from http://www.ncbi.nlm.nih.gov/tools/primer-blast/) [21]. All oligonucleotides were supplied unmodified and desalted (Microsynth AG, Switzerland).

Table 1

Specification of used reference genes

Gene symbol	Sequence (5'->3')	Position	Product length	Intron inclusion	RTqPCR efficiency in DH/DRG	Ref.
Full name
Accession number
18S	Fw: GGCTCATTAAATCAGTTATGGTTCCT	94-119	147	no	78%/70%	[22]
18S ribosomal RNA	Rev: GTTGGTTTTGATCTGATAAATGCACG	240-215
V01270
RPL29	Fw: ACAGAAATGGCATCAAGAAACCC	96-118	105	yes	81%/77%	[15]
Ribosomal protein L29	Rev: TCTTGTTGTGCTTCTTGGCAAA	200-179
NM_017150
RPL13a	Fw: TCTCCGAAAGCGGATGAACAC	185-205	145	yes	88%/73%	[14]
Ribosomal protein L13A	Rev: CAACACCTTGAGGCGTTCCA	329-310
NM_173340
Actb	Fw: GGAGATTACTGCCCTGGCTCCTA	1023-1045	150	yes	85%/63%	[14]
Actin, beta	Rev: GACTCATCGTACTCCTGCTTGCTG	1172-1149
NM_031144
HPRT1	Fw: GCATCTAAGAGGTTTCCCCAGT	1133-1154	76	no	70%/74%	1
Hypoxanthine phosphoribosyl-transferase 1	Rev: GCATTTAAAAGGAACGGTTGAC	1208-1187
NM_012583
HMBS	Fw: GAGACCATGCAGGCCACCAT	1007-1026	97	yes	78%/72%	[15]
Hydroxymethyl-bilane synthase	Rev: TTGGAATGTTCCGGGCAGTG	1084-1103
NM_013168
GAPDH	Fw: CCCCCAATGTATCCGTTGTG	780-799	118	yes	89%/87%	[23]
Glyceraldehyde-3-phosphate	Rev: TAGCCCAGGATGCCCTTTAGT	877-897
dehydrogenase
NM_017008

Fw forward, Rev reverse primer sequence, DH dorsal horn, DRG dorsal root ganglion.

1 from our laboratory database.

Quantitative real-time PCR (qPCR)

We performed qPCR on an iQ5 Cycler (Bio-Rad) with SYBR Green I. The reactions were carried out in 96-well plates (Thermo-Fast® 96 Semi-Skirted PCR Plate, Thermo Scientific, Switzerland) each with a total volume of 20 μl. Each well contained 5 μl of a 100-fold dilution of cDNA, 10 μl of iQ™ Sybr® Green Supermix (2× qPCR mix contains dNTPs, 50 U/ml iTaq DNA polymerase, 6 mM MgCl2, SYBR Green I, enhancers, stabilizers, 20 nM fluorescein) (Bio-Rad, Switzerland), 2 μl of each primer 1–3 μM and 1 μl water. We optimized qPCR conditions on the iQ5 thermal gradient cycler and by testing different concentrations of primers and templates. The qPCR program began with an initial three-minute denaturation step at 95°C to activate the hot-start iTaq™ DNA polymerase. This was followed by 45 cycles of 10 s at 95°C for denaturation and 45 s at 60°C for annealing and extension.

We confirmed the amplification of specific qPCR products by performing a melting-curve step at the end of each run. Serial dilution curves for each primer allowed us to calculate qPCR efficiencies. The 100-fold diluted cDNA utilized for all the amplifications was within the linear dynamic range of the calibration curve – between 10 and 1000-fold dilution. Across all the assays, none of the quantification cycle (Cq) values was higher than 30. No-template and no-reverse transcription controls were run to determine any contamination or the generation of primer dimers. All amplifications were run in triplicate, and any doubtful curves were excluded. To minimize technical variation between samples through different runs we preferred the sample maximization method [11], i.e. a run contained all the samples for one gene of interest respective to one reference gene.

Data analysis

Raw data was collected and computed by iQ5 BioRad software with an automated analysis of the baseline and threshold of each run, and then exported in Excel files for further analyses. We rescaled all Cq values for each gene to the lowest Cq value as an internal control, converted these rescaled Cq logarithmically into linear, relative quantities taking into account the gene specific amplification efficiency [relative quantity = (1+efficiency) ^ (Cqinternal control-Cqsample)]. Finally, we calculated arithmetical means from the replicates [11,24]. To explore the stability of each candidate gene we calculated M-values from the geNorm algorithm and coefficients of variation (CV) [9,11]. The M-value corresponds to the average pairwise variation between a gene and the other candidate genes. The more this value tends to zero, the more stably the gene is expressed in comparison to the others. For a heterogeneous sample, a candidate gene with an M-value below 1 can be considered as a reliable reference gene [11]. The CV for a gene was obtained by calculating the ratio of the standard deviation (σ) to the mean (μ) of relative quantities (CV = σ / μ). For a heterogeneous panel, a value below 50% is proposed as satisfactory [11].

To rank the genes that satisfied the above criteria, M-values were recalculated after stepwise elimination of the worst candidate. To assess the number of reference genes needed for an accurate normalization, the effect of inclusion of a supplementary gene to those already considered stable was determined by calculating the pairwise variation of the added candidate gene to the others (standard deviation of logarithmically transformed expression ratios). If this pairwise variation goes below 0.15, no supplementary gene is needed [9].

Results

Reference genes used in the SNI model

We found 484 articles in the PubMed database using “spared nerve injury” as keywords. Twenty-six articles were retrieved with semi-quantitative PCR or RT-qPCR performed in the rat DRG, spinal cord or brain. Despite its rigorous validation never having been mentioned, GAPDH was the most commonly used reference gene (16 times) whether in the DRG or in the dorsal horn. This was followed by Actb and 18S (3 times), HPRT1 (twice) and HMBS (once). HPRT1 and 18S were the only two reference genes with validation in the rat SNI model (Table  2).

Table 2

List of reference genes in the rat spared nerve injury model

Article	PCR method	Sample	Reference genes	Validation
Kanda et al., 2013 [25]	PCR	L4/L5 SC	GAPDH	no
Zhou et al., 2013 [26]	RT-qPCR	L4/L5 DRG	GAPDH	no
Kashimoto et al., 2013 [27]	RT-qPCR	L4/L5 SC	GAPDH	no
Shankarappa et al., 2012 [28]	RT2 Profiler PCR array	L4/L5 DRG	HPRT1	validated1
Inquimbert et al., 2012 [29]	PCR	L4/L5 DH	GAPDH	no
Zapata et al., 2012 [30]	RT-qPCR	RVM	18S and HMBS	no
Samad et al., 2013 [31]	PCR	L4 DRG	GAPDH	no
Kobayashi et al., 2012 [32]	RT-qPCR	L4/L5 SC	GAPDH	no
Liu et al., 2012 [33]	RT-qPCR	L4–L6 SC	Actb	no
Tochiki et al., 2012 [34]	PCR	L4-L6 DH	HPRT1	no
Okubo et al., 2012 [35]	RT-qPCR	L4-L5 SC	GAPDH	no
Del Rey et al., 2011 [36]	PCR	hippocampus	not mentioned	no
Kühlein et al., 2011 [37]	RT-qPCR and PCR	DH and DRG	Actb and 18S	no
Yamanaka et al., 2011 [38]	RT-qPCR	L4, 5 DRG	GAPDH	no
de Novellis et al., 2011 [39]	RT-qPCR	mPFC	Actb	no
Vega-Avelaira et al., 2009 [40]	PCR	L4/L5 DRG	GAPDH	no
Costigan et al., 2009 [41]	PCR	L4/L5 DH	GAPDH	no
Okubo et al., 2010 [42]	RT-qPCR	L4-L5 SC	GAPDH	no
Staaf et al., 2009 [43]	Taqman Low Density Arrays	L4 DRG	18S	validated2
Moss et al., 2008 [44]	PCR	L4/5 DRG	GAPDH	no
Berta et al., 2008 [45]	PCR	L4 and L5 DRGs	GAPDH	no
Millecamps et al., 2007 [46]	PCR	mPFC	GAPDH	no
Moss et al., 2007 [47]	PCR	L4/L5 dorsal horn	GAPDH	no
Apkarian et al., 2006 [48]	RT-qPCR	brainstem, thalamus, and prefrontal cortex	GAPDH	no
Pertin et al., 2005 [49]	PCR	L4/L5 DRG	GAPDH	no
Takahashi et al., 2003 [50]	RT-qPCR	L5 DRG	GAPDH	no

Literature search in the PubMed database of studies using semi-quantitative PCR or RT-qPCR in the rat SNI model.

PCR semi-quantitative polymerase chain reaction, RT-qPCR reverse transcription quantitative real-time polymerase chain reaction, SC spinal cord, DRG dorsal root ganglion, DH dorsal horn, RVM rostral ventromedial medulla, PL-IL prelimbic and infralimbic, mPFC medial prefrontal cortex.

1 HPRT1 validated among 5 reference genes.

2 See Discussion.

Validation of reference genes in spinal cord dorsal horn over time after SNI

Rat spinal cord dorsal horn samples were analyzed in an overall panel without differentiating any timepoints, i.e. throughout a temporal follow-up with naive animals, 2nd, 4th, 7th, 10th, 14th and 21st postoperative day animals. The pairwise comparison of all the potential reference genes (Actb, 18S, GAPDH, RPL13a, RPL29, HPRT1 and HMBS) calculated using geNorm resulted in M-values below 0.5 for all but the 18S gene, which showed 0.52 (Figure  1). All genes were below the M-value cut-off of 1.0 proposed for a heterogeneous panel. CV ranged from 20% for HMBS to 34% for GAPDH and RPL13a, i.e. below the 50% cut-off for stability in a heterogeneous panel.

Figure 1

Analysis of reference gene stability. Assessment of reference gene stability, in the ipsilateral dorsal horn (A) and L4-L5 dorsal root ganglia (B) in the rat spared nerve injury model using M-value (left) from the geNorm algorithm and coefficient of variation (CV) (right). Candidate genes are classified from left to right by increasing order of M-value, from most to least stable. The cut-off for stability is 1 for M-value and 50% for CV (dotted line).

Validation of reference genes in DRG after SNI in a time-dependent panel

As for the dorsal horn, DRG samples were analyzed in an overall panel without differentiating any timepoints. M-values were below 1.0 and CV below 50% for Actb, GAPDH, RPL13a, RPL29, HPRT1 and HMBS. The 18S reference gene had the worst M-value, at 0.8, but this was still acceptable within the selection. However, we rejected it from our gene sample because of its 61% CV (Figure  1).

Ranking of reference genes

We reanalyzed the 7 tested reference genes (Actb, GAPDH, RPL13a, RPL29, HPRT1, HMBS and 18S) to rank them. After each step, the gene with the highest M-value was eliminated until we obtained the two most stable genes from the list. No further discrimination was possible.

In the spinal cord dorsal horn, the 2 most stable genes remaining after the stepwise elimination were RPL29 and RPL13a, and then in increasing order GAPDH, Actb, HMBS, HPRT1 and 18S (Figure  2A).

Figure 2

Ranking of reference genes. Stepwise calculation of stability by excluding the least stable gene and averaging the remaining M-values in the spinal cord dorsal horn (A) and the dorsal root ganglia (B). Insert: list of reference genes ranked from the most to least stable.

In the DRG, the 2 most stable genes remaining after the stepwise elimination were HPRT1 and Actb, and then in increasing order RPL29, RPL13a, HMBS, GAPDH and 18S (Figure  2B). The latter is used because this ranking only takes into account the M-value and not the CV.

Minimal number of reference genes

In the dorsal horn samples, the calculation of pairwise variations resulted in values much lower than 0.15 along the entire time course and in all samples. Using the two best reference genes the pairwise variation in dorsal horn was 0.07 (Figure  3). The addition of a supplementary reference gene did not significantly lower the pairwise variation. Results for DRG were similar, with a low pairwise variation of 0.05 if only the two most stable reference genes were taken in consideration.

Figure 3

Determination of the optimal number of control genes for accurate normalization. Determination of the number of reference genes to consider for an accurate normalization in spinal cord dorsal horn (A) and dorsal root ganglia (B) by calculating pairwise variation (V) after stepwise elimination of the worst reference gene.

We performed the same analysis using the least stable reference genes on our list for dorsal horn and DRG, and their pairwise variation was still below 0.15.

Discussion

Here, we tested the stability of seven candidate reference genes according to the geNorm algorithm to ensure accurate RT-qPCR normalization. M-values below 1.0 together with CV lower than 50% for Actb, GAPDH, RPL13a, RPL29, HPRT1, and HMBS identified these candidate genes as good controls in the rat SNI model in the spinal cord dorsal horn and DRG. The 18S gene could not be validated in DRG because of a CV above 50%.

After SNI, gene expression changes in both sites can be enormous compared to the baseline and also over time [51]. Because of this heterogeneity, meticulous RT-qPCR procedures relying on MIQE recommendations are essential to obtain consistent data. The most frequently chosen reference gene for RT-qPCR normalization after SNI in the rat is GAPDH, but its expression level had, to our knowledge, never been confirmed as being stable in this model (Table  2). Here, we confirm that GAPDH is a good candidate gene, even if it is not the most stable in our series.

The only 2 reference genes with previously mentioned validation in the rat SNI model are HPRT1 and 18S. HPRT1 was used by Shankarappa et al. in DRG [28] at 2 timepoints after injury and chosen from among 5 different reference genes on an array. After showing that 18S had a very low variation of Cq values in all data sets, from both SNI and naive groups, Staaf et al. used it to normalize and compare the expression of the transient receptor potential (TRP) family of genes in DRG after SNI. They mentioned neither ranking between the different reference genes, nor the criteria used to choose 18S from among their list of them [43]. These demonstrations are part of the workflow of RT-qPCR arrays where several reference genes are proposed, and their analysis is usually fully integrated in the procedure.

GeNorm is a popular algorithm for assessing different reference genes from a given panel [9,11]. The calculation of M-values permits measurement of an average variation in the expression ratios of one potential reference gene with the others and consequently informs us whether the tested gene is constantly expressed compared to them. The CV of relative quantities of a gene merely reflects the dispersion of its expression level. It has been previously demonstrated that a CV below 25% was typically observed in stable reference genes as long as the experimental conditions were homogeneous. In a heterogeneous panel, such as in the present study, a CV below 50% is acceptable [11]. As it measures the gene-specific variation, the CV in assessing the stability of a reference gene assumes a standardization of the procedures for all samples. Taken together, M-values and the CV are very useful tools for determining the most accurate reference genes because the information that they convey is different. All our candidate genes satisfied an M-value cut-off below 1.0 for a heterogeneous panel. 18S, despite its acceptable M-value, failed to satisfy the CV cut-off in DRG. Considering these results 18S could not be validated as a stable reference gene.

The more the M-value tends to zero, the more the candidate gene is stable compared to the other tested genes. Hence the stability of the different genes can be ranked. A first approximate classification was established after analyzing all the samples together (Figure  1). Actb, RPL13a and RPL29 were the 3 best genes in both tissues. A potentially more accurate classification can be performed if M-values are recalculated in a series of steps, after elimination of the worst scoring gene at each step. At the end, the two most reliable reference genes remain, but they cannot be further discriminated by definition of the M-value (Figure  2). Interestingly the ranking differed between dorsal horn and DRG. We again found RPL 29 and RPL13a to be the most stable genes in the spinal cord dorsal horn, however, in the DRG, HPRT1 and Actb were the best, followed immediately by RPL13a and RPL29. In the present study, differences in the ranking were very slight, and the former “rough” method would have been sufficient. In other circumstances, or with another panel of candidate genes, it could be useful to compare them further until a definitive ranking is obtained, so as to choose only the best control genes from amongst the candidates. M-value calculations and subsequent ranking were also performed independently for each different time point (data not shown) without being able to identify any time-dependent regulation. The primary goal of this study was to validate some suitable RT-qPCR reference genes at different time points, in different tissues, in a widely used neuropathic pain model.

According to MIQE guidelines, the use of a single reference gene is not enough to compensate for the intrinsic variation within that reference gene and normalization against several reference genes should be performed [7,9]. To define the minimal number of necessary reference genes, the pairwise variation permits an assessment of whether the removal of the worst reference gene may result in a greater variation (Figure  3). Here, in both DRG and dorsal horn, the two best candidate genes were sufficient to ensure an accurate normalization. This demonstration allows us to state that in our model, using only the two best reference genes would be sufficient: namely RPL 29 and RPL13a in the spinal cord dorsal horn and HPRT1 and Actb in the DRG. However, even the worst pair satisfied stability criteria, with a pairwise variation less than 0.15. In fact, all pairs of validated reference genes (not 18S because of its CV above 50%) provided a robust validation in both tissues and at all time points. GAPDH is often used, but according to MIQE guidelines it should not be used alone.

The criterion for choosing a reference gene is more than just its stability. 18S, for example, is disputed as a reference gene by some authors who argue that it by far exceeds the mRNA quantity of most genes of interest and it might be absent from pure mRNA extracts [9].

Conclusions

We validated 7 reference genes in the spinal cord dorsal horn: (Actin beta (Actb), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), ribosomal proteins L13a (RPL13a) and L29 (RPL29), hypoxanthine phosphoribosyltransferase 1 (HPRT1), hydroxymethylbilane synthase (HMBS) and 18S. In DRG, 18S did not fulfill all our criteria and could not be retained as a good reference gene, although the other 6 were. We were able to show that in both tissues, using any two of the proposed stable reference genes is sufficient to ensure an accurate validation. We were also able to confirm the stability of GAPDH, which is the most commonly used reference gene in the SNI model, but it should be used in combination with another control gene.

Availability of supporting data

The data sets supporting this article’s results are available in the Labarchives repository under http://dx.doi.org/10.6070/H4CZ3532.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

NP performed the surgeries, the tissue extraction, the RT-qPCR, the literature search, the statistical analysis and drafted the manuscript. ID participated in the study design and in the drafting of the manuscript. MRS conceived the study, participated in its design and coordination and supervised the drafting of the manuscript. All authors read and approved the final manuscript.