- Research article
Reliable reference miRNAs for quantitative gene expression analysis of stress responses in Caenorhabditis elegans
BMC Genomicsvolume 15, Article number: 222 (2014)
Quantitative real-time PCR (qPCR) has become the “gold standard” for measuring expression levels of individual miRNAs. However, little is known about the validity of reference miRNAs, the improper use of which can result in misleading interpretation of data.
Here we undertook a systematic approach to identify highly stable miRNAs in different stress conditions such as low oxygen (hypoxia), UV-stress and high temperature (heat-stress) in the nematode Caenorhabditis elegans. We conducted genome-wide RNA-seq for small RNAs and selected abundant miRNAs with minimal variation of expression between the different conditions. We further validated the stable expression of a selection of those constitutively expressed candidates in the different stress conditions by SYBR Green qPCR. The selected miRNA candidates were analyzed for stability by applying the widely used geNorm logarithm. With this approach, we were able to successfully identify suitable reference miRNAs for each stress condition. Interestingly, we also found that 3 miRNAs, namely mir-2-5p, mir-46-3p and mir-47-3p, are stable in all the above-mentioned conditions suggesting that they might have general functions independent of stress.
Our analysis offers a comprehensive list of stably expressed miRNAs in different stress conditions that can be confidently used as reference miRNAs for qPCR analysis in C. elegans.
The level at which molecules are expressed in a biological sample, under different conditions, often holds information about their mechanistic role(s). Quantitative real-time PCR (qPCR) has been established as one of the standard methods for semi-quantitative measurements of the abundance of various transcripts in virtually any biological sample. However, qPCR has been the source of many misleading reports due to being prone to experimental errors . One of the most commonly encountered mistakes is the choice of unsuitable reference genes, which can lead to erroneous conclusions . Reference genes are transcripts against which qPCR results are normalized. Currently, a number of normalization strategies for qPCR results exist, however, the use of reference genes is the most accepted one . Reference genes have to exhibit stable expression between the conditions of interest in order to serve as suitable normalizers.
Caenorhabditis elegans has been widely used as an experimental model for a range or studies, spanning from development to drug discovery. Gene expression analysis in this model organism constitutes a powerful tool to discover new roles for different types of molecules, such as miRNAs. miRNAs are small RNA molecules (~22 nt) that mostly impact on gene expression post-transcriptionally  and have been implicated in different stress responses [5, 6]. qPCR has extensively been used to assess expression levels of both protein coding genes and miRNAs in C. elegans. A previous study reported the validation and selection of appropriate reference genes for qPCR for protein coding genes , however, there is no report of stable miRNAs that can be used as reference genes in qPCR analysis of C. elegans miRNAs.
miRNA data normalization has traditionally relied on expression of other small RNAs, such as small nuclear (snRNAs) or small nucleolar (snoRNAs) RNAs (e.g. U8, U6 or RNU44). Using such transcripts is not ideal because these small RNAs are of different size and function, and are produced and/or processed by different pathways to miRNAs. Moreover, the use of only one or few popular reference genes in each and every qPCR experiment without experimental validation of the stable expression of these genes between the conditions used, can lead to heavily distorted results and misinterpretations.
Here we undertook a systematic approach to identify stable miRNAs across a range of conditions (Figure 1). We chose to focus on three different, widely-studied, stress conditions in C. elegans, namely low oxygen (hypoxia), UV-stress and high temperature (heat-stress). We provide a selection of stably expressed and reliable reference genes for normalizing qPCR results obtained from stress studies in C. elegans.
A battery of miRNAs is stably expressed throughout different stress conditions
In an effort to identify miRNAs that are stably expressed in different stress conditions, we conducted genome-wide small RNA-seq on samples treated with hypoxia, UV-irradiation or high temperature, three widely studied stress conditions in C. elegans (Figure 1 and Additional file 1: Table S1). We then filtered the results for candidates that showed minimal variation of expression between unstressed conditions and each stress condition. We arbitrarily set a threshold of accepted expression variation to 7%, based on normalized number of reads for each miRNA in each sample. Furthermore, we discarded miRNAs with less than 500 reads as lowly expressed (Additional file 2: Table S2). We found 17 stably expressed miRNAs for UV-treated samples, 16 for the high temperature-treated samples, 32 for the hypoxia-treated samples and 28 for the hypoxia-treated samples, including samples of hif-1 mutant animals (Figure 2, Additional file 3: Table S3). The hif-1 mutant is a highly relevant strain for studying hypoxic responses and including hif-1 samples in the study aimed to reveal hif-1-independent, stably expressed miRNAs. Notably, we found mir-2-5p, mir-46-3p and mir-47-3p to be stably expressed in all conditions (Figures 2 and 3 and Additional file 4: Figure S1, Additional file 5: Figure S2, Additional file 6: Figure S3). Thus, the high-throughput RNA-seq study provided a number of miRNAs exhibiting stable expression in each of the three stresses when compared to control conditions and three potential ‘general stress’ miRNA reference genes. These highly stable miRNAs may be involved in the regulation of basic cellular functions that are not influenced by any of the stresses used, thus their expression remains unchanged.
Validation of stable expression of selected miRNAs by an inexpensive, selective and sensitive qPCR method and subsequent geNorm analysis
The use of three different reference genes has recently been shown to be sufficient for accurate normalization of results obtained by qPCR analysis for protein-coding genes in C. elegans[7, 8], as was previously suggested . For this reason we decided to validate the stable expression of three miRNAs in each stress condition, choosing candidates that span over a wide range of abundance (low, medium and high levels), as well as to validate the three miRNAs that are stable between all the conditions. We chose to use an optimized method for miRNA qPCR that is based on a single reverse transcription reaction for all miRNAs followed by a real-time qPCR reaction with a pair of DNA based miRNA specific primers . This technique is highly sensitive and can recognize between miRNAs that differ by as little as one nucleotide . Moreover, this technique is fairly inexpensive. We believe that the identified candidates would also be appropriate reference genes to be used with other qPCR chemistries, such as TaqMan or SYBR Green with LNA primers.
miRNAs selected for validation of stable expression for each sample group
UV: For the UV-treated samples, we chose mir-57-5p as highly expressed (215390,119 reads), mir-55-3p as intermediately expressed (8240,595 reads) and mir-357-3p as lowly expressed (1289,386 reads) to validate by qPCR (Figure 4).
Heat-stress: For the high temperature-treated samples, we selected mir-64-5p as highly expressed (163324,299 reads), mir-86-5p as intermediately expressed (6457,365 reads) and mir-794-5p as lowly expressed (636,433 reads) to validate by qPCR (Figure 4).
Hypoxia N2: For the hypoxia-treated samples, we selected mir-52-5p as highly expressed (19899,848 reads), mir-241-3p as intermediately expressed (1086,805 reads) and mir-254 as lowly expressed (585,136 reads) to validate by qPCR (Figure 4).
Hypoxia N2 + hif-1: For the hypoxia-treated samples including the hif-1 mutant, we chose mir-58-3p as highly expressed (375488,779 reads), mir-230-3p as intermediately expressed (2681,376 reads) and mir-357-3p as lowly expressed (1287,386 reads) to validate by qPCR (Figure 4).
All conditions: mir-47-3p (21043,167 reads), mir-46-3p (22186,815 reads) and mir-2-5p (5977,944 reads) were stably expressed in all the conditions and therefore chosen to validate by qPCR (Figure 4).
The detected Cq (Quantification cycle) values for all but one (mir-57-5p) of the above miRNAs correlated with the read counts obtained by RNA-seq. mir-57-5p is found to be more abundant than mir-55-3p in UV-treated samples when tested by qPCR, probably due to RNA-seq inaccuracy (Figure 4).
geNorm analysis results
The relative expression levels of all the selected miRNAs were calculated based on the Cq values obtained, with the primer efficiencies taken into account. Subsequently, the geNorm logarithm was used to calculate the average gene expression stability (M value) of each of the miRNAs. The more stable the expression of a gene is, the lower is its M value. All the miRNAs tested in our study showed an M value less than 0.5 (excluding mir-58-3p - M value equals 0.573), which is the suggested cut-off of geNorm corresponding to high expression stability. The coefficience of variation for all of the candidates was < 0.2. Based on M values we conclude that the selected miRNAs are stably expressed between the conditions tested and can confidently be used as reference genes in qPCR analysis in C. elegans, when the particular stress conditions are under study. The order for each group of miRNAs based on M value is shown in Figure 3.
Next, we wanted to determine the optimal number of reference genes to be used in qPCR studies of stress responses, as assessed by the pairwise variation value (Vn/n + 1 value). The V2/3 value for all but ‘Hypoxia N2 + hif-1′ miRNA groups was <0.15, which indicates that the use of the two best miRNAs, as previously evaluated by M value, is sufficient for accurate normalization (Figure 3). In ‘Hypoxia N2 + hif-1′ group the V2/3 value was >0.15 (0.189) and therefore qBase suggests the use of all 3 miRNAs (Figure 3).
Primer efficiency and specificity evaluation
All the primer pairs used for the qPCR reactions (Table 1) were highly efficient as calculated by the use of sequential dilutions of a representative cDNA sample (Table 2). Moreover, the specificity of the primers was assessed by standard melting curve analysis (Additional file 7: Figure S4). All primer pairs gave a single melting peak indicating their specificity in amplifying only one product. ‘No reverse transcription’ controls as well as ‘no template controls’ showed no amplification suggesting no DNA contamination in the samples and no contamination in the reagents used for qPCR experiment, respectively. Moreover, the RNA used was of excellent quality as assessed by Nanodrop and Bioanalyzer, a crucial prerequisite for reliable gene expression estimation (Additional file 8: Figure S5).
In recent years, significant effort has been directed towards re-evaluating the way that qPCR experiments are conducted , with much focus on the identification of suitable reference genes in different experimental conditions. In this study, we used genome-wide RNA-seq for small RNAs as an unbiased approach to discover stably expressed miRNAs in a range of different stress conditions. This approach allowed us to identify a number of miRNAs with stable expression levels ranging from 500 to many thousands of reads (Additional file 2: Table S2). We further evaluated a subset of these miRNAs by applying the commonly used geNorm logarithm on qPCR measurements. All the miRNAs selected for further evaluation showed high expression stability. This strong consistency between RNA-seq results and qPCR analysis suggests that the stable miRNAs found by RNA-seq could all be potentially used as reference genes in qPCR analysis. In fact, a previous sequencing study has shown that two of the miRNAs identified in our study, mir-46-3p and mir-47-3p, are remarkably stable when expression levels are compared between young adults and aged animals . However, we strongly advise the prospective user to perform qPCR validation in each experimental setup and test the stability of chosen miRNAs using geNorm or other similar analysis tool (e.g. NormFinder , BestKeeper ), before confidently using them as reference genes. Nevertheless, we have validated the stable expression of a sufficient number of miRNAs for each condition, which the user can confidently use in qPCR analysis.
The method of choice for performing the qPCR was described by Balcells et al. . This method has key advantages; it is simple (single RT-reaction) and inexpensive (regular DNA primers) with high sensitivity and specificity. Although we used this method to validate stable miRNAs, we believe that the identified candidates would be appropriate reference genes to be used with other qPCR chemistries, such as TaqMan or qPCR with LNA primers.
Our study provides a comprehensive list of miRNAs that can confidently be used for normalization of qPCR data in C. elegans stress biology. Moreover, it demonstrates an unbiased approach to identify and validate appropriate reference miRNAs based on their stable expression in genome-wide RNA seq data.
Worm maintenance and worm sample preparation
Mixed populations of worms were grown on E. coli (strain OP50) seeded NGM (Nematode-Growth-Medium) agar plates for 2–3 days at 20°C as previously described . Worms were washed off the plates and bleached to obtain eggs. These eggs were kept in sterile M9 buffer (3 g KH2PO4, 6 g Na2HPO4, 5 g NaCl, 1 ml 1 M MgSO4, H2O to 1 liter) overnight to hatch and developmentally arrest at the L1 stage. Starvation staged L1 worm populations were grown on OP50 seeded NGM plates until the L4 stage. For the hypoxia-treated samples, mid-L4 animals were placed in a hypoxic chamber with 0.5% O2 for 4 h and harvested immediately. For the heat-stress-treated samples, mid-L4 animals were heated at 35°C for 4 h and harvested immediately. For the UV-treated samples, mid-L4 animals were irradiated with 120 J/m2 UV-C and harvested 4 h post-treatment. Samples were harvested by washing with M9. One set of samples was used for the small RNA-seq and three biologically independent sets of samples for the qPCR validation.
Worm pellets were kept in Trizol and lysed by freeze-cracking. A round of chloroform phase separation was followed by an addition of 1.5 volumes of 100% ethanol. Consecutive miRNeasy column purification was performed as per supplier instructions starting from step nine in the miRNeasy guidebook (Qiagen) . RNA concentration and purity was assessed by Nanodrop (Thermo Scientific) and followed by RNA integrity assessment using Bioanalyzer RNA Nano Chips (Agilent Technologies). The RIN (RNA Integrity Number) of the total RNA samples ranged from 10.0-7.1 with 10 being the highest possible value (Additional file 8: Figure S5).
Genome-wide small RNA-seq
The preparation of the RNA samples for small RNA-seq and the analysis of sequencing data was performed as described in Friedlander et al. and is detailed below.
Small RNA library preparation for sequencing
Small RNA library preparation was obtained from 1 μg of total RNA per sample. 5′-monophosphate dependent small RNA cloning and library preparation was then performed with the TruSeq Small RNA Sample Preparation Kit (Illumina, USA). Cluster generation as well as sequencing of the prepared libraries was performed on the Illumina cluster station and GAIIx (Illumina, USA).
miRNA fold-change analysis
Preprocessing of raw sequencing data was performed with a custom PERL script to group reads with different barcodes into separate fasta files. These fasta files were processed with the mapper.pl module from miRDeep2  to clip sequencing adapters, collapse reads and map them to the C. elegans genome (WS205) with the following command: mapper.pl config.txt -c -d -i -j -k TGGAATTCTCGGGTGCCAAGG -l 18 -m -p c_elegans.WS205.dna -s reads.fa -t reads_vs_genome.arf -v. Raw read counts as output by the quantifier module from the miRDeep2 package were normalized to the total number of miRNA reads per barcode then multiplied by the mean value of all miRNA reads of all barcode files. A pseudocount of 10 was added and the normalized values were log2 transformed. A regression line with slope 1 was calculated in R for all pairwise sample comparisons and the deviation from each data point to this line (residual) was used to estimate the fold change between two samples.
qPCR DNA primers were designed according to previously published instructions  using a custom-made primer design software that can be downloaded from Aalborg University server:
“http://vbn.aau.dk/en/publications/mirprimerdesign3(bcd08649-3883-4016-b493-52dceedd04f8).html.” (Courtesy of Peter Kamp Busk, University of Aalborg, Denmark). A list of the primers used can be found in the Table 1.
cDNA synthesis and quantitative real-time PCR
cDNA synthesis was essentially performed as previously described . Briefly, single tube poly (A) synthesis and reverse transcription was performed with 100 ng of RNA in a final volume of 10 μl combined with 1 μl MuLV reverse transcriptase (New England Biolabs), 1 μl of poly (A) polymerase (New England Biolabs), 1 μl 10x poly (A) polymerase buffer, 1 μl 0.1 mM ATP, 1 μl 10 μM RT-primer (5′-CAGGTCCAGTTTTTTTTTTTTTTTVN where V is A, C and G and N is A, C, G and T). The reaction was incubated at 42°C for 1 hour followed by enzyme inactivation at 95°C for 5 minutes. A panel of three technical cDNA replicates for each of the three biological samples per condition was constructed. The cDNA was diluted 5× before used in qPCR reactions.
miR-specific qPCR was performed as previously described in . Briefly, 0.25 μl of 10 μM miRNA specific forward and 0.25 μl of 10 μM miRNA specific reverse primer, 5 μl of 2× LightCycler 480 SYBRGreen I Mastermix (Roche) and 3.5 μl of RNAse free H2O and 1 μl of cDNA of each samples was combined in 10 μl reaction and subjected to qPCR performed on LightCycler 480 (Roche).
Sequential dilutions of a representative cDNA sample were made in order to calculate the primer pair efficiency (E) using the formula E% = (10(1/slope)-1) × 100. The Cq values obtained from the qPCR runs were applied into the geNorm algorithm  as part of the qBase software package  to calculate the expression stability value (M value) and the pairwise variation Vn/n + 1. The M value is defined as the average pairwise variation of a single candidate reference gene to all other tested candidates. The commonly accepted cut off of 0.15 for the pairwise variation was used to assess the minimum number of reference genes required for accurate normalization .
Baker M: qPCR: quicker and easier but don’t be sloppy. Nat Methods. 2011, 8 (3): 207-212. 10.1038/nmeth0311-207.
Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, Mueller R, Nolan T, Pfaffl MW, Shipley GL, Vandesompele J, Wittwer CT: The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009, 55 (4): 611-622. 10.1373/clinchem.2008.112797.
Huggett J, Dheda K, Bustin S, Zumla A: Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 2005, 6 (4): 279-284. 10.1038/sj.gene.6364190.
Pasquinelli AE: NON-CODING RNA MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship. Nature Reviews Genetics. 2012, 13 (4): 271-282.
Pocock R: Invited review: decoding the microRNA response to hypoxia. Pflugers Arch. 2011, 461 (3): 307-315. 10.1007/s00424-010-0910-5.
Kagias K, Nehammer C, Pocock R: Neuronal responses to physiological stress. Front Genet. 2012, 3: 222-
Hoogewijs D, Houthoofd K, Matthijssens F, Vandesompele J, Vanfleteren JR: Selection and validation of a set of reliable reference genes for quantitative sod gene expression analysis in C. elegans. BMC Mol Biol. 2008, 9: 9-10.1186/1471-2199-9-9.
Zhang Y, Chen D, Smith MA, Zhang B, Pan X: Selection of reliable reference genes in Caenorhabditis elegans for analysis of nanotoxicity. Plos One. 2012, 7 (3): e31849-10.1371/journal.pone.0031849.
Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3 (7): RESEARCH0034-
Balcells I, Cirera S, Busk PK: Specific and sensitive quantitative RT-PCR of miRNAs with DNA primers. BMC Biotechnol. 2011, 11: 70-10.1186/1472-6750-11-70.
de Lencastre A, Pincus Z, Zhou K, Kato M, Lee SS, Slack FJ: MicroRNAs both promote and antagonize longevity in C. elegans. Curr Biol. 2010, 20 (24): 2159-2168. 10.1016/j.cub.2010.11.015.
Andersen CL, Jensen JL, Orntoft TF: Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64 (15): 5245-5250. 10.1158/0008-5472.CAN-04-0496.
Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP: Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper–excel-based tool using pair-wise correlations. Biotechnol Lett. 2004, 26 (6): 509-515.
Brenner S: The genetics of Caenorhabditis elegans. Genetics. 1974, 77 (1): 71-94.
Cirera S: Highly efficient method for isolation of total RNA from adipose tissue. BMC research notes. 2013, 6: 472-10.1186/1756-0500-6-472.
Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N: miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012, 40 (1): 37-52. 10.1093/nar/gkr688.
Hellemans J, Mortier G, De Paepe A, Speleman F, Vandesompele J: qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007, 8 (2): R19-10.1186/gb-2007-8-2-r19.
We thank Nikolaus Rajewsky and Sebastian Mackowiak for performing the RNA-seq experiment, Susanna Cirera for guidance regarding RNA isolation and qPCR setup as well as Peter Kamp Busk for assistance with qPCR primer design. This work was supported by a grant from the European Research Council (ERC Starting Grant number 260807 to R.P.).
The authors declare no competing financial or non-financial interests.
Conceived, designed and performed the experiments: AP and KK. Analyzed the data: KK and AP. Wrote the paper: KK, AP and RP. All authors read and approved the final manuscript.
Konstantinos Kagias, Agnieszka Podolska contributed equally to this work.