The dilution effect and the importance of selecting the right internal control genes for RT-qPCR: a paradigmatic approach in fetal sheep
BMC Research Notes volume 8, Article number: 58 (2015)
The key to understanding changes in gene expression levels using reverse transcription real-time quantitative polymerase chain reaction (RT-qPCR) relies on the ability to rationalize the technique using internal control genes (ICGs). However, the use of ICGs has become increasingly problematic given that any genes, including housekeeping genes, thought to be stable across different tissue types, ages and treatment protocols, can be regulated at transcriptomic level. Our interest in prenatal glucocorticoid (GC) effects on fetal growth has resulted in our investigation of suitable ICGs relevant in this model. The usefulness of RNA18S, ACTB, HPRT1, RPLP0, PPIA and TUBB as ICGs was analyzed according to effects of early dexamethasone (DEX) treatment, gender, and gestational age by two approaches: (1) the classical approach where raw (i.e., not normalized) RT-qPCR data of tested ICGs were statistically analyzed and the best ICG selected based on absence of any significant effect; (2) used of published algorithms. For the latter the geNorm Visual Basic application was mainly used, but data were also analyzed by Normfinder and Bestkeeper. In order to account for confounding effects on the geNorm analysis due to co-regulation among ICGs tested, network analysis was performed using Ingenuity Pathway Analysis software. The expression of RNA18S, the most abundant transcript, and correlation of ICGs with RNA18S, total RNA, and liver-specific genes were also performed to assess potential dilution effect of raw RT-qPCR data. The effect of the two approaches used to select the best ICG(s) was compared by normalization of NR3C1 (glucocorticoid receptor) mRNA expression, as an example for a target gene.
Raw RT-qPCR data of all the tested ICGs was significantly reduced across gestation. TUBB was the only ICG that was affected by DEX treatment. Using approach (1) all tested ICGs would have been rejected because they would initially appear as not reliable for normalization. However, geNorm analysis (approach 2) of the ICGs indicated that the geometrical mean of PPIA, HPRT1, RNA18S and RPLPO can be considered a reliable approach for normalization of target genes in both control and DEX treated groups. Different subset of ICGs were tested for normalization of NR3C1 expression and, despite the overall pattern of the mean was not extremely different, the statistical analysis uncovered a significant influence of the use of different normalization approaches on the expression of the target gene. We observed a decrease of total RNA through gestation, a lower decrease in raw RT-qPCR data of the two rRNA measured compared to ICGs, and a positive correlation between raw RT-qPCR data of ICGs and total RNA. Based on the same amount of total RNA to performed RT-qPCR analysis, those data indicated that other mRNA might have had a large increase in expression and, as consequence, had artificially diluted the stably expressed genes, such as ICGs. This point was demonstrated by a significant negative correlation of raw RT-qPCR data between ICGs and liver-specific genes.
The study confirmed the necessity of assessing multiple ICGs using algorithms in order to obtain a reliable normalization of RT-qPCR data. Our data indicated that the use of the geometrical mean of PPIA, HPRT1, RNA18S and RPLPO can provide a reliable normalization for the proposed study. Furthermore, the dilution effect observed support the unreliability of the classical approach to test ICGs. Finally, the observed change in the composition of RNA species through time reveals the limitation of the use of ICGs to normalize RT-qPCR data, especially if absolute quantification is required.
Quantitative reverse transcription real time polymerase chain reaction (RT-qPCR) is the most sensitive and accurate method to measure mRNA gene expression. Accurate normalization is critical in accounting for quantity of input RNA, varying amounts of cDNA input, sample loss during handling, and activity and variation in the kinetics of reverse-transcription enzyme during reaction [1,2]. A reliable normalization may be even more important when the experiment aims to evaluate low-fold changes of a target gene .
Among several proposed method of normalization, including total RNA, genomic DNA, and spike-in of an external artificial reference , the use of internal control genes (ICGs) appear to be the most accurate . Validated ICG are essential for unbiased interpretation of RT-qPCR given the vast number of confounding factors that can affect the transcription of any gene . Classically, housekeeping genes (HKG) were used as ICGs because it was thought that they expression is stable across different tissue types, ages, and treatments; however, several studies revealed that also the expression of HKG is under active regulation . This was suggested by some studies where it was shown that expression of HKG may vary as a result of neoplastic growth, hypoxia or other experimental treatments [6-9]. Among classic HKG, beta-actin (ACTB) mRNA expression for example is increased by a maximum of 4-fold in fibroblasts during 8 hrs of serum stimulation . ACTB and cyclophilin A (PPIA) expression in human tumor cell lines also varied widely after hypoxia treatment . The expression of ribosomal protein, large, P0 (RPLP0) was significantly changed after hindlimb suspension intervention in rat tendon and muscle . With lipopolysaccharide stimulation, the expression level of hypoxanthine phosphoribosyltransferase 1 (HPRT1) and tubulin, beta (TUBB) increased in alveolar macrophages from chronic obstructive pulmonary disease patients . All the above clearly demonstrated that the use of classic HKG for RT-qPCR normalization without any assessment is highly unreliable. In addition, they also indicated that the ICGs do not need to be HKG. Therefore, there was a need of a new approach for testing the reliability of ICGs for each experiment. In order to evaluate the reliability of ICGs for RT-qPCR normalization several algorithms have been developed; those include geNorm Visual Basic (geNorm v3.5) , Bestkeeper , and NormFinder . Studies dealing with very different conditions, including tissue/organ development, have successfully used geNorm in selecting reliable ICGs [14-22], a method that has been statistically validated . Those studies have shown the importance of using the geometric mean of at least two stably expressed ICGs. In addition, some of those studies have also demonstrated the detrimental effect of improper selection of ICGs on RT-qPCR results of target genes of interest. The geNorm program assesses the stability of gene expression between candidate reference genes by performing a pair-wise comparison of the expression ratio. The fundamental rationale of the geNorm algorithm is that the higher the stability of expression ratio between two non-coregulated candidate genes across the samples, the higher the likelihood that those are stably expressed; thus, highly reliable ICGs. Therefore, it is critical prior to geNorm analysis to verify the absence of co-regulation among potential ICGs which otherwise would bias geNorm results .
Despite the critical importance of proper selection and evaluation of ICGs, these are often chosen based upon previously scientific publications rather than upon empirical data. Furthermore, the evaluation of ICGs reliability is performed by analyzing changes in the raw RT-qPCR in respect to the experimental conditions (e.g., gestational age, sex or the effects of a special treatment). This approach, however does not take the accumulation of errors and variations in each single sample during the analytical procedure or the dilution effect of stably expressed genes due to large increase in expression of very abundant transcripts into account .
Prenatal synthetic glucocorticoid (GC) treatment is commonly used in the management of women at risk of early preterm birth and has succeeded in reducing neonatal mortality and morbidity . Despite these benefits, GC treatment has also been associated with a decrease in birth weight and alterations in glucose homeostasis and hypothalamic–pituitary–adrenal (HPA) axis function and fetal programming [26-28]. Although the exact mechanisms underlying those responses are unknown, fetal hepatic development may play an important role . The liver, largest of the body’s organs, plays an important role in coordinating metabolic homeostasis, nutrient processing and detoxification . Many studies suggest the liver to be a target organ of GC treatment and partly responsible for fetal growth restriction [26,31,32] and gene expression profiling might help dissecting underlying processes. Early maternal GC treatment (dexamethasone, DEX) is used in suspected cases of congenital adrenal hyperplasia to protect female fetuses from virilization . We have shown that DEX treatment analogous to that used in human subjects in the first third of pregnancy, resulted in profound changes in fetal HPA axis development which persisted during postnatal [26,34].
In order to study the effect of early DEX treatment on liver gene expression using RT-qPCR it is essential to have reliable ICGs. In previous studies involving sheep under different experimental conditions, internal references such as 18S ribosomal RNA (RNA18S), ACTB, HPRT1, RPLP0, PPIA and TUBB have been widely used to investigate a range of target genes. All those ICGs were chosen in these studies based solely on previous publications in different experimental models [35-38] and not a thorough evaluation of reliability was performed. There are no previous data on the selection of reliable ICGs in fetal sheep liver which are independent of gestational age, sex and not regulated by GC treatment.
Within this framework we set out to identify reliable ICGs to be used for normalize RT-qPCR data from the experiment that investigate the effects of early DEX treatment on the fetal sheep liver transcriptome during gestation. For this we have tested the 6 ICGs previously used in sheep with the aim to test the expression stability, identify the most reliable, and uncover the number of ICGs that should be used for accurate RT-qPCR normalization. We evaluated the reliability of ICGs following two approaches: (1) we performed a statistical analysis of the raw RT-qPCR data to assess effect of DEX, gestation, sex, and their interaction (i.e., “the classical analysis”) and (2) by testing the stability and reliability of the normalization by using algorithms, particularly geNorm Visual Basic application (“geNorm analysis”). The effect on normalization based on the two approaches and different subsets of ICGs was assessed on the RT-qPCR data of glucocorticoid receptor (NR3C1).
Approach 1: Expression levels of ICGs, “the classical analysis”
Raw RT-qPCR data of all ICGs tested significantly (p < 0.05) decreased due to gestation between 50 and 140dG in all tested ICGs; however patterns were slightly different between ICGs (Figure 1). RPLPO, TUBB and PPIA significantly decreased between 50 and 140dG (Figure 1A, B, F), whereas ACTB levels stayed stable between 50 and 100dG (p > 0.05), but significantly decreased between 100 and 140dG (Figure 1C, D). HPRT1 and RNA18S levels stayed stable between 50 and 125dG and significantly decreased at 140dG (Figure 1E). Compared to control, TUBB raw RT-qPCR data was significantly lower at 125dG in female treated with DEX; raw RT-qPCR data of the other ICGs were not affected by DEX treatment (Figure 1). Gender differences were only found in RPLPO at 125dG, with significantly higher expression levels in females as compared to males (Figure 1A). In summary, all ICGs tested were affected by gestation and some by DEX and/or gender. Among the tested ICGs the raw RT-qPCR of HPRT1 and RNA18S were the least affected by the tested conditions; thus, can be considered the best choice for normalization of RT-qPCR of target genes among the 6 tested ICGs.
Approach 2: Determination of ICGs expression stability using algorithms
The raw RT-qPCR data of the 6 tested ICGs were analyzed with geNorm, Bestkeeper, and Normfinder. The rank from the most to the least reliable ICG was similar between the three algorithms (Table 1), with PPIA being one of the two most reliable ICGs among the ones tested. In geNorm the lowest M value indicates ICGs with the most stable expression. In the overall analyses, stepwise elimination of successive genes showed that PPIA and HPRT1 were the most stable ICGs across gestation followed by RNA18S and RPLPO (Figure 2). The determination of the optimal number of ICGs for normalization is performed by geNorm by calculating the pairwise variation (V-value) of adding the subsequent more reliable gene. A V-value below 0.15, cut-off reported by Vandesompele et al. , was obtained by adding the 5th more reliable gene (V4/5); however, the addition of the 4th more reliable gene (i.e., V3/4) had a V-value of 0.157, which is similar to V4/5 (Figure 2). Based on this observation and based also on practicality, we deemed that the use of 4 most reliable ICGs among the one tested, that is PPIA, HPRT1, RNA18S and RPLPO, can provide a trustworthy normalization factor.
geNorm is one of several algorithms available to evaluate ICGs. We have evaluated our potential ICGs also using Bestkeeper and Normfinder (Table 2). There was an overall agreement about the 3 algorithms to rank the best vs. the worst ICGs tested with some small differences, the major one being RNA18S as the most reliable ICG in Bestkeeper when considering the calculated raw RT-qPCR data instead of Ct values. The latter is the method originally used to develop Bestkeeper; thus, likely more reliable .
The geNorm algorithm is based on the assumption that the tested genes are completely independent; i.e., they are not co-regulated because they do not have up-stream common regulator(s). Therefore, it is essential to verify that the best ICGs are indeed independent. For this, potential co-regulation was analyzed using Ingenuity Pathway Analysis (IPA). The use of IPA revealed a co-regulation (or common up-stream regulator(s)) between RNA18S and ACTB via tumor protein p53 (TP53), methyl CpG binding protein 2 (MECP2), and methyl-CpG binding domain protein 2 (MBD2) [39-46]. Direct co-regulation has been observed between TUBB, RPLPO and ACTB via v-myc myelocytomatosis viral (MYCN) [47-49] and between TUBB and RPLPO via FBJ murine osteosarcoma viral (FOS) [49,50]. However, no direct co-regulation has been observed between HPRT1, PPIA, RNA18S and RPLPO (Figure 3), the ICGs uncovered by geNorm to be reliable, and the most reliable among the one tested.
Normalization of a target gene by using best ICG(s) uncovered in approach 1 and 2
We have compared the results of normalizing RT-qPCR data of glucocorticoid receptor (NR3C1), used as target gene, using the geometrical mean of the 4 best ICGs as indicated by approach 2 (Figure 4A), the two best ICGs (or the more “flat” ICGs HPRT and RNA18S) indicated by the approach 1 (Figure 4B and 4C) and ACTB, the ICG with the lowest average expression stability among the one tested but also one of the most used ICGs in literature (Figure 4D). The statistical differences observed between comparisons on the quantity of NR3C1 mRNA expression levels obviously differed between subsets of ICGs used for normalization. For example, NR3C1 mRNA expression significantly increased between 50 and 100dG and further increased between 100 and 140dG when normalized by the geometrical mean of HPRT1, PPIA, RNA18S and RPLPO (Figure 4A). However, when normalized only to RNA18S, the increase between 100 and 140dG was not significant (Figure 4B). Normalizing to ACTB resulted in a significant increase of NR3C1 mRNA expression between 50 and 125dG (Figure 4D). Normalizing to HPRT1, the ICG with the least time effect on the raw RT-qPCR data resulted in a significant reduction of the quantity of NR3C1 mRNA expression levels at 140dG in DEX compared to controls (Figure 4C), which was not significant when normalized to the other subsets of ICGs (Figure 4A, B, D).
Dilution effect and ICGs
The total RNA was significantly affected by gestation similar to mean raw RT-qPCR data of all ICGs tested (Figure 5A). For this, it is not surprising that we have observed a significant positive correlation between RNA concentration and raw RT-qPCR data for all ICGs tested with exception of HPRT1, where the P-value of correlation was 0.06 (Table 2). The raw RT-qPCR data between all ICGs (with exception of RNA18S) and insulin-like growth factor 1 (IGF1) and glucose 6 phosphatase (G6PC) were negatively correlated (Table 2). The raw RT-qPCR data of very abundant expressed ribosomal gene RNA18S had the least correlations with the liver-specific genes and total RNA, while raw RT-qPCR data of RNA28S did not have any correlation (Table 2). This was mostly due to the very limited effect of gestation on the raw RT-qPCR data (Figure 5B).
Reliability and accuracy of target gene expression levels in RT-qPCR greatly depends on the selection of ICGs; therefore a set of appropriate ICGs must be determined. An ideal ICG for our experiment would be stably expressed under all conditions, i.e., is independent of gestational age, gender or treatment. The previously reported ICGs used for RT-qPCR normalization in fetal sheep liver lacked empirical data demonstrating suitability where many studies have used a single ICG [32,35,52], and the stability of ICGs was not measured or reported.
Various strategies have been applied to normalize for the amount of starting material, differences in tissues, enzymatic efficiencies or overall transcriptional activity. The traditional way of analyzing the accuracy of ICGs is to pick out one or two widely reported ICGs and to test whether experimental conditions affect ICG genes expression levels. Our present study analyzed the effect of gender and DEX treatment at four different gestational time points in fetal liver in six commonly used ICGs in sheep. Significant gestation effect was observed in all measured ICGs. A similar association was found in bovine mammary tissue, in which ICG expression levels significantly changed during the whole lactation . Despite being all ICGs tested affected by gestation, RNA18S and HPRT1 were the least affected. Others have reported changes in expression of many target genes in fetal sheep after GC treatment [32,35,53]. However, in the present study, only raw RT-qPCR data of TUBB in females was affected by DEX treatment. Based on the above data and using the criteria of the classical approach (i.e., no effect in expression by the studied condition) we should have rejected all tested ICGs as being unreliable.
Although statistical effects of gestational age for all the ICGs tested were observed in our study, this did not exclude the possibility of a reasonable fit for these genes as ICGs, as indicated in other studies [1,19]. The absence of a statistical effect on potential ICGs is not an essential condition to consider them reliable for normalization and, oppositely, the use of such criteria to select proper ICGs is a serious limitation, especially in experiments dealing with temporal transcriptomics during differentiation [1,19,21]. In fact, the classical approach has two major limitations for selection of reliable ICGs: 1) relying exclusively on the lack of differences between means does not account for the sample-specific variation (i.e., the normalization is performed per each single sample and the use of the means does not account for that specific variation); 2) cannot account for the dilution effect (see below). Very often associated with such an approach there is also the use of a single ICG, which presents several limitations in itself .
In order to overcome such limitations, in the current study, we relying on algorithm, particularly geNorm, to evaluate the most stable ICGs and acquire accurate data on the optimal number of ICGs that should be used for normalization . Contrary to the use of the classical approach, geNorm analysis of pairwise variation among the 6 evaluated ICGs in the present study revealed that most of the ICGs tested can be considered reliable and that the use of the geometrical mean of 4 ICGs (PPIA, HPRT1, RNA18S and RPLPO) can provide a robust normalization.
The effect of an inappropriate set of ICGs for normalization of a target gene was demonstrated by normalizing the RT-qPCR data of a target gene (NR3C1) using the two ICG with the lowest time effect (i.e., the most ‘flat’ genes), i.e., RNA18S and HPRT1 or the 4 most stable ICGs indicated by geNorm analysis among the pool of ICG candidates chosen (geometrical mean of HPRT1, PPIA, RNA18S and RPLPO). Despite the fact that all the tested genes had a similar pattern (Figure 1), the analysis revealed obvious differences in the statistical results of the normalized target gene with the best subset of ICGs or the use of only one of the genes least affected by time (RNA18S and HPRT1). The difference was even more pronounced when the normalization was performed by the very commonly used ACTB (Figure 4). This comparison clearly confirmed the limitation of using the classical approach to select ICGs or the use of previously published ICGs without validation using specific algorithms.
The results from the use of geNorm (approach 2) contradict the results from the classical approach (approach 1). In approach 1 all the tested ICGs would have been rejected based on the fact that the raw RT-qPCR data were significantly affected by gestation and/or DEX treatment. All the tested ICGs had a decrease in raw RT-qPCR data through gestation (Figure 1). How is it possible that ICGs can be both deemed to be reliable and have the raw RT-qPCR data significantly affected by a condition (by definition an ICG should not be affected by any condition)? In a previous study in bovine mammary tissue from pregnancy to end of subsequent lactation it was observed a similar pattern of ICGs . In that study it was demonstrated that this pattern is an artifact and consequence of a dilution effect due to equal amount of total RNA is used for RT-qPCR between samples concomitant to a very large expression of few abundant mammary-specific genes with a consequent relative decrease in abundance of the stably expressed genes. In such a case the apparent decrease in raw expression of stably expressed genes is due to the decrease of proportional abundance of their mRNA among all RNAs in the sample and not to active transcription repression.
Based on the similarity between the previous study in bovine mammary tissue and the present study, we have assessed if a dilution effect could have been the cause of the observed reduction of raw RT-qPCR data during gestation (Figure 1; Additional file 1: Figure S1). For this, and similarly to the previous study , we have evaluated if a correlation existed between raw RT-qPCR data of ICGs and total RNA concentration. We have also measured the expression of RNA28S, the most abundant RNA, and correlated with ICGs. In addition, we have run a correlation with raw RT-qPCR data of the liver-specific genes IGF1  and G6PC . Contrary to the previous study , in the present study a significant positive correlation was observed among the raw RT-qPCR data of the tested ICGs and the amount of RNA per μg of tissue. As for the ICGs tested the abundance of RNA decreased through time (Figure 5A). In the same time period the measured amount of rRNA, which is known to make up the majority of total RNA, had a minor change through the differentiation, although tended to decrease. Most of the ICGs decreased more than 3-fold, the rRNAs decreased less than 2-fold, while the total RNA decreased 2-fold (Figure 5A and 5B). Despite the observed decrease through time and correlations, the data are uncoupled. In fact, the total RNA can be considered an absolute value while the values for ICGs and the rRNA are relative. This uncoupling (or lack of direct relationship) is due to the fact that the starting amount of total RNA was equal between all samples (see material and methods) despite the change in total RNA. Therefore, the apparent large decrease in raw RT-qPCR of ICGs and, even with lower magnitude, of the rRNA, strongly support that something else composing the total RNA must have increased dramatically during differentiation and was not measured in the present experiment.
The above observations allow us to conclude that at the beginning of pregnancy there was an overall higher RNA expression leading to higher total RNA concentrations compared to the end of pregnancy. The expression of RNA includes not only mRNA but all RNA species. A similar observed decrease in most of mRNA has been previously reported during the differentiation of human hepatocytes (HepaRG) in vitro where a decrease in expression of most of the genes was observed . Therefore, there was an active depression of transcription during differentiation, likely driven by epigenetic specialization . Based on the negative correlations of total RNA and ICGs with liver specific genes, we postulate that there was a large increase in liver specific genes, likely very abundant, such as albumin and fatty acid binding protein 1. This large increase of abundant liver-specific gene can have caused the dilution of other genes, especially the genes with a constant copy number/cells, such is the case for the tested ICGs. Increase in liver-specific abundant genes was also previous observed  and are known to be under active transcriptional control by specific liver-specific transcription regulators .
A dilution of stably expressed genes in the present experiment can explain the observed apparent, but likely misleading, decrease in raw RT-qPCR data of ICGs (Figure 1 and Additional file 1: Figure S1). Not accounting for such effect, would bias the final normalization of RT-qPCR data of target genes. In our case if a target gene that has an increase in mRNA copy/cell is subject to the dilution effect and it probably appear as having a lower or no change through time if it is normalized by a “flat” gene (i.e., not subjected to the dilution, but, for this has in reality an up-regulation in copy mRNA/cell that is inversely proportional to the dilution). The normalization of RT-qPCR data of that same target gene by using reference genes that proportionally follow the dilution can correct for that dilution revealing the true pattern of the target gene.
The above interpretation of the data has a caveat. If the total mRNA is changing through differentiation this means that there is an overall increase or decrease in overall transcription, as also suggested by previous data from others , resulting in an increased or decreased amount of mRNA/cell. This means that the proportion of the 3 types of RNA can change. This observation reveals a limitation on relying only on ICGs for RT-qPCR normalization. The normalization by ICGs is methodologically correct if a relative expression between specific mRNA is considered, but it is not appropriate if an absolute or real expression (i.e., amount of specific mRNA/cell) is considered. The caveat discussed above highlights a shortcoming of using expression of other mRNA to normalize target mRNAs and prompts to find even better approach for RT-qPCR normalization.
In the present analysis using geNorm we have uncovered that the geometrical mean of HPRT1, PPIA, RNA18S and RPLPO can be suitable for calculating the normalization factor to be used with RT-qPCR data in developing fetal sheep liver. The data also indicated that the use of one ICG, even among the list of the 4 more stable ICGs, to normalize a target gene can bias the final results. For this we recommend using the geometric mean of at least the four indicated ICGs as an accurate normalization factor. However, this can be considered reliable only in the present experimental setup. Different experiments, even using the same model, require a separate validation of ICGs, since those are experiment-specific.
In addition, the data clearly indicated that, in order to have a proper relative quantification of target genes, the normalization must be able to account for the dilution effect and also for the variation in total RNA and ratio between RNA species. In the present case the data supported a dilution effect that affected all genes but most visible in stably expressed genes; therefore, proper normalization needs to account for such dilution effect and reliable reference genes should have a decrease through gestation in measured raw RT-qPCR data.
Finally, the present data appear to highlights a limitation of using ICGs for normalization when a change in overall transcription during tissue differentiation is observed, i.e. increase/decrease of mRNA as was suggested by the data in the present experiment. Therefore, even if the proposed normalization can be considered reliable for relative gene expression, the quest for an absolute way to normalize RT-qPCR data is ever more necessary.
Animals and tissues
All procedures were approved by the Animal Experimentation Ethics Committee of the University of Western Australia and/or the Western Australian Department of Agriculture .
Pregnant Merino ewes (Ovis aries) with singleton pregnancies (total n =106 liver samples) of known gestational age were allocated randomly to receive maternal injections of DEX or saline (control). Maternal DEX (Mayne Pharma, Victoria, Australia) injections were given in a dose of 0.14 mg/kg ewe body weight suspended in 2 mL of saline, consisting of four intramuscular injections at 12 h intervals over 48 h between 40 and 41 days of gestation (dG). Control animals received saline injections of a comparable volume (2 mL saline/ewe) .
Tissue samples were collected at 49–51 (50), 101–103 (100), 125–127 (125), and 140–142 (140) dG. Fetal liver and other major fetal organs were removed, weighed, and collected for use in other studies (e.g., ). Right liver lobe (100, 125 and 140dG) and total liver (50dG) biopsies were snap frozen in liquid nitrogen before storage at −80°C.
Real time PCR protocol
Total liver RNA was extracted using the RNeasy Midi kit (QIAGEN, Clifton Hill, Victoria, Australia) and stored at-80°C until further use. RNA concentration in μg/μg tisse was determined by NanoDrop 1000 Spectrophotometer (average 260/280 = 2.04 ± 0.08; Thermo Fisher Scientific, Wilmington, U.S.A.) and integrity was analyzed in a selected number of samples with Agilent 2100 BioAnalyzer (RNA 6000 Nano Kit 5067–1511; RIN mean 8.2). Possible genomic DNA contamination was removed from each sample using a DNase treatment (Fermentas, Thermo fisher scientific, Catalog #EN0521), and then RNA sample (1 μg) was reverse transcribed in a 20 μl reaction mixture (iScript cDNA synthesis kit, BIO-RAD, Catalog #170-8890) according to the manufactures manual. DNase treatment and RT reactions were carried out in a Mastercycler (Eppendorf, Germany). For each run a no template control sample containing no RNA was reverse transcribed to provide a negative control for real time PCR applications. For RT-qPCR, primer pairs for sheep (Table 3) were either designed using Primer 6.0 (TUBB, PPIA, ACTB, RNA28S) or were previously reported (RPLP0 , RNA18S , HPRT1 , NR3C1 , IGF1 , G6PC ). All primer sequences have been tested and verified with fluorescent color band sequencing (PCR purification kit #83050, Seqlab, Sequence Laboratories Göttingen, Germany).
The RT-qPCR analysis was run in triplicates on an ABI 7500 Real Time PCR System (Applied Biosystems, Foster, USA). All measurements were run in a total volume of 10 μl including 5 μl power SYBR green master mix (Applied Biosystems, Foster, USA), 4 μl cDNA, 0.4 μl primer pair in a final concentration of 400 nM and 0.2 μl water. A non template control was included in each RT-qPCR run for each primer. After gradient PCR optimization with a 10 step annealing temperature gradient according to the manufactures protocol (Mastercycler gradient, Eppendorf Germany), all transcripts were amplified using the following cycling conditions: 95°C for 10 min for one cycle and 95°C for 20 sec, 60°C (RPLPO, TUBB, RNA18S, ACTB, HPRT1, PPIA, IGF1, G6PC and RNA28S) or 58.9°C (NR3C1) for 32 sec, and 72°C for 60 sec for 40 cycles (Table 3). Melting-curve analysis (95°C for 15 sec, 60°C for 60 sec, 95°C for 15 sec and 60°C for 15 sec for one cycle) demonstrated a single PCR product in all measurements. Efficiencies of all genes were calculated in the same sample with E(%) = 100x(101/slope -1)  with same ABI 7500 PCR program and deemed to be similar to each other (range from 95.3% to 100%; see Table 3).
Quantification and ICG stability evaluation
The geomean of threshold cycle (Ct) values from the 7500 system software SDS software version 1.4 (Applied Biosystems, Foster, USA) were exported to Microsoft Excel. The quantity of each ICG was calculated by Q = EfficiencyΔCt where ΔCt = Ct min-Ct sample for each gene. Statistical analyses were performed by using SPSS 20 statistical software (SPSS Inc., Chicago, USA). Data was tested for normal distribution and equal variance (Levene test, p > 0.05). Data that were not normally distributed (RPLPO, TUBB, ACTB, HPRT1, PPIA, NR3C1) were log transformed to achieve normality. To determine treatment, day of gestation, and gender effects as well an interaction between them on raw or normalized RT-qPCR data, a MANOVA with treatment, gender and day of gestation as factors, followed by a pairwise comparison (Holm’s Sidak) when main effects were p < 0.05 was performed. Main effects are indicated in the figure legend; post hoc p-values (Holm-Sidak) are indicated in figures. Data are presented as mean ± S.E.M. Statistical significance was accepted for values p < 0.05. The relationship between the average expression of ICGs and RNA concentration (ug/ug tissue) was analyzed with the Pearson correlation.
To determine the reliability of ICGs using the geNorm analysis the average expression stability values (M) of 6 ICG genes in all samples were analyzed with geNorm Visual basic application (V 3.5, Biogazelle NV, Zwijnaarde, Belgium) according to the manufactures manual and the procedures described by Vandesompele et al.  with respect to days of gestation, gender and treatment. The raw RT-qPCR data were then uploaded into geNorm, Bestkeeper, or Normfinder and ICGs ranked based on algorithm-specific values. For geNorm the ICGs were ranked based on the M value, which refers to the constancy of the expression ratio between two genes among all samples tested . At least two ICGs are required for stability analysis and subsequent M values are calculated prior to the stepwise exclusion of the least stable ICGs. The lowest M value indicates the pair of ICGs with the most stable expression.
In order to determine the minimum number of ICGs which has to be used to obtain a reliable normalization factor (NF), geNorm calculates the pairwise variation Vn/n + 1 between NFn and NFn + 1 from the NF calculated using the two most stable genes (i.e., with lower M) and the NF calculated using the additional more stable gene (i.e., n + 1) and so on until the stepwise inclusion of subsequent less stable ICGs has no significant contribution to the newly calculated NF (Figure 2). For example, V2/3 shows the variation of the NF of two genes in relation to three genes. geNorm suggests a cut-off for the pairwise variation V of ≤ 0.15. A relatively large decrease in pairwise variation V indicates that the added gene has a significant effect on the final NF and should be included in the calculation of NF. However, the proposed cut off at ≤0.15 of the V is considered to be indicative and the lower the V the higher the reliability of NF . Ideally the NF should be calculated by the geometrical mean of the combination of ICGs with the lower V-value. According to geNorm protocol the use of the three best ICGs is in most cases a valid normalization strategy, and results are much more accurate and reliable compared to the use of only one ICG .
It is known, that the geNorm approach has the problem that it tends to select genes with similar expression profiles . Therefore, the pairwise comparison approach of geNorm to determine ICG stability is highly biased by potential co-regulation between selected ICGs [1,65]. The potential ICGs should not be regulated through common upstream effectors or should not directly regulate each other . We therefore analyzed in a post-run check potential co-regulation between the 6 potential ICGs with Ingenuity Pathway Analysis web-based software, (Ingenuity System accessed 12.12.2012, default settings were used, one relationship step was analyzed; www.ingenuity.com, Version 14400082, Build 192063, Redwood City, CA). The software allows the discovery, exploration and visualization of gene interactions and co-regulations. Networks are generated relying on known published relationships among human, mouse and rat genes. Briefly, the network analysis was performed using only relationship types including “expression”, “transcription” and “any direction”. The “Build-Grow” option was used and all the up-stream transcription factors of our ICGs were added. The following options were selected: Interactions = “direct” and “indirect”; Grow out”All the molecules” that are “Upstream of the selected molecules” and molecules were limited to “Use Ingenuity Knowledge Base”; Relationship Types = “expression” and “transcription”; Molecule Types = “ligand-dependent nuclear receptor” and “transcription regulator”. All the other options were left as default. Once all the upstream molecules had been added by IPA, we manually deleted all the transcription factors that had only one relationship with our ICGs.
Quantities of the exemplary target gene (NC3R1) was calculated by: Q = efficiency^(ct min - ct sample), rescaled normalized expression levels Qnormalized/rescaled = (Qsample/NFsample)/Min(Qsample/NFsample) (geNorm v3.5 manual ), where NF was calculated using the number of most reliable ICGs as uncovered by geNorm.
18S ribosomal rRNA
28S ribosomal rRNA
FBJ murine osteosarcoma viral
Hypoxanthine phosphoribosyltransferase 1
Internal control gene
Average expression stability value
v-myc cytelomatosis viral oncogene homolog
v-myc myelocytomatosis viral
Peptidylprolyl isomerase A
Ribosomal protein, large, P0
Reverse transcription quantitative polymerase chain reaction
Tumor protein p53
Defines the pairwise variation between sequential NFs
Bionaz M, Loor JJ. Identification of reference genes for quantitative real-time PCR in the bovine mammary gland during the lactation cycle. Physiol Genomics. 2007;29:312–9.
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 biology 2002, 3:RESEARCH0034.
Drewlo S, Levytska K, Sobel M, Baczyk D, Lye SJ, Kingdom JC. Heparin promotes soluble VEGF receptor expression in human placental villi to impair endothelial VEGF signaling. J Thrombosis Haemostasis: JTH. 2011;9:2486–97.
Huggett J, Dheda K, Bustin S, Zumla A. Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 2005;6:279–84.
Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55:611–22.
Zhong H, Simons JW. Direct comparison of GAPDH, beta-actin, cyclophilin, and 28S rRNA as internal standards for quantifying RNA levels under hypoxia. Biochem Biophys Res Commun. 1999;259(3):523–6.
Thellin O, Zorzi W, Lakaye B, De Borman B, Coumans B, Hennen G, et al. Housekeeping genes as internal standards: use and limits. J Biotechnol. 1999;75(2):291–5.
Schmittgen TD, Zakrajsek BA. Effect of experimental treatment on housekeeping gene expression: validation by real-time, quantitative RT-PCR. J Biochem Biophys Methods. 2000;46(1–2):69–81.
Wu YY, Rees JL. Variation in epidermal housekeeping gene expression in different pathological states. Acta Derm Venereol. 2000;80(1):2–3.
Heinemeier KM, Olesen JL, Haddad F, Schjerling P, Baldwin KM, Kjaer M. Effect of unloading followed by reloading on expression of collagen and related growth factors in rat tendon and muscle. J Appl Physiol. 2009;106(1):178–86.
Ishii T, Wallace AM, Zhang X, Gosselink J, Abboud RT, English JC, et al. Stability of housekeeping genes in alveolar macrophages from COPD patients. Eur Respir J. 2006;27(2):300–6.
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:509–15.
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:5245–50.
Garcia-Crespo D, Juste RA, Hurtado A. Selection of ovine housekeeping genes for normalisation by real-time RT-PCR; analysis of PrP gene expression and genetic susceptibility to scrapie. BMC Vet Res. 2005;1(1):3.
Ohl F, Jung M, Xu C, Stephan C, Rabien A, Burkhardt M, et al. Gene expression studies in prostate cancer tissue: which reference gene should be selected for normalization? J Mol Med (Berl). 2005;83(12):1014–24.
Erkens T, Van Poucke M, Vandesompele J, Goossens K, Van Zeveren A, Peelman LJ. Development of a new set of reference genes for normalization of real-time RT-PCR data of porcine backfat and longissimus dorsi muscle, and evaluation with PPARGC1A. BMC Biotechnol. 2006;6:41.
Ayers D, Clements DN, Salway F, Day PJ. Expression stability of commonly used reference genes in canine articular connective tissues. BMC Vet Res. 2007;3:7.
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.
Monaco E, Bionaz M, Sobreira de Lima A, Hurley WL, Loor JJ, Wheeler MB. Selection and reliability of internal reference genes for quantitative PCR verification of transcriptomics during the differentiation process of porcine adult mesenchymal stem cells. Stem Cell Res Ther. 2010;1:7.
Piantoni P, Bionaz M, Graugnard DE, Daniels KM, Akers RM, Loor JJ. Gene expression ratio stability evaluation in prepubertal bovine mammary tissue from calves fed different milk replacers reveals novel internal controls for quantitative polymerase chain reaction. J Nutr. 2008;138:1158–64.
Tramontana S, Bionaz M, Sharma A, Graugnard DE, Cutler EA, Ajmone-Marsan P, et al. Internal controls for quantitative polymerase chain reaction of swine mammary glands during pregnancy and lactation. J Dairy Sci. 2008;91:3057–66.
Kadegowda AK, Bionaz M, Thering B, Piperova LS, Erdman RA, Loor JJ. Identification of internal control genes for quantitative polymerase chain reaction in mammary tissue of lactating cows receiving lipid supplements. J Dairy Sci. 2009;92:2007–19.
Szabo A, Perou CM, Karaca M, Perreard L, Palais R, Quackenbush JF, et al. Statistical modeling for selecting housekeeper genes. Genome Biol. 2004;5(8):R59.
Warrington JA, Nair A, Mahadevappa M, Tsyganskaya M. Comparison of human adult and fetal expression and identification of 535 housekeeping/maintenance genes. Physiol Genomics. 2000;2:143–7.
Roberts D, Dalziel S. Antenatal corticosteroids for accelerating fetal lung maturation for women at risk of preterm birth. Cochrane Database Syst Rev. 2006;3, CD004454.
Li S, Sloboda DM, Moss TJ, Nitsos I, Polglase G, Doherty CB, et al. Effects of glucocorticoid treatment given in early or late gestation on growth and development in sheep. JDOHaD. 2013;4(2):146–56.
Moss TJ, Sloboda DM, Gurrin LC, Harding R, Challis JR, Newnham JP. Programming effects in sheep of prenatal growth restriction and glucocorticoid exposure. Am J Physiol Regul Integr Comp Physiol. 2001;281(3):R960–970.
Seckl JR. Prenatal glucocorticoids and long-term programming. Eur J Endocrinol/Eur Federation Endocrine Soc. 2004;151 Suppl 3:U49–62.
Sloboda DM, Newnham JP, Challis JR. Repeated maternal glucocorticoid administration and the developing liver in fetal sheep. J Endocrinol. 2002;175:535–43.
Hyatt MA, Budge H, Symonds ME. Early developmental influences on hepatic organogenesis. Organogenesis. 2008;4(3):170–5.
Nyirenda MJ, Lindsay RS, Kenyon CJ, Burchell A, Seckl JR. Glucocorticoid exposure in late gestation permanently programs rat hepatic phosphoenolpyruvate carboxykinase and glucocorticoid receptor expression and causes glucose intolerance in adult offspring. J Clin Invest. 1998;101(10):2174–81.
Sloboda DM, Newnham JP, Challis JR. Repeated maternal glucocorticoid administration and the developing liver in fetal sheep. J Endocrinol. 2002;175(2):535–43.
Lajic S, Nordenstrom A, Hirvikoski T. Long-term outcome of prenatal treatment of congenital adrenal hyperplasia. Endocr Dev. 2008;13:82–98.
Li S, Nitsos I, Polglase GR, Braun T, Moss TJ, Newnham JP, et al. The effects of dexamethasone treatment in early gestation on hypothalamic-pituitary-adrenal responses and gene expression at 7 months of postnatal age in sheep. Reprod Sci. 2012;19:260–70.
Braun T, Li S, Sloboda DM, Li W, Audette MC, Moss TJM, et al. Effects of Maternal Dexamethasone Treatment in Early Pregnancy on Pituitary-Adrenal Axis in Fetal Sheep. Endocrinology. 2009;150(12):5466–77.
Ashley RL, Arreguin-Arevalo JA, Nett TM. Binding characteristics of the ovine membrane progesterone receptor alpha and expression of the receptor during the estrous cycle. Reprod Biol Endocrinol. 2009;7(1):42.
Fletcher CJ, Roberts CT, Hartwich KM, Walker SK, McMillen IC. Somatic cell nuclear transfer in the sheep induces placental defects that likely precede fetal demise. Reproduction. 2007;133(1):243–55.
Ma Y, Zhu MJ, Uthlaut AB, Nijland MJ, Nathanielsz PW, Hess BW, et al. Upregulation of growth signaling and nutrient transporters in cotyledons of early to mid-gestational nutrient restricted ewes☆. Placenta. 2011;32(3):255–63.
Boiko A, Porteous S, Razorenova O, Krivokrysenko V, Williams B, Gudkov A. A systematic search for downstream mediators of tumor suppressor function of p53 reveals a major role of BTG2 in suppression of Ras-induced transformation. Genes Dev. 2006;20(2):236–52.
Ginsberg D, Mechta F, Yaniv M, Oren M. Wild-type p53 can down-modulate the activity of various promoters. Proc Natl Acad Sci U S A. 1991;88:9979–83.
Santhanam U, Ray A, Sehgal P. Repression of the interleukin 6 gene promoter by p53 and the retinoblastoma susceptibility gene product. Proc Natl Acad Sci U S A. 1991;88(17):7605–9.
Cohen S, Gabel H, Hemberg M, Hutchinson A, Sadacca L, Ebert D, et al. Genome-Wide Activity-Dependent MeCP2 Phosphorylation Regulates Nervous System Development and Function. Neuron. 2011;72:72–85.
Di Como C, Gaiddon C, Prives C. p73 function is inhibited by tumor-derived p53 mutants in mammalian cells. Mol Cell Biol. 1999;19:1438–49.
Zhai W, Comai L. Repression of RNA polymerase I transcription by the tumor suppressor p53. Mol Cell Biol. 2000;20:5930–8.
Ng H, Zhang Y, Hendrich B, Johnson C, Turner B, Erdjument-Bromage H, et al. MBD2 is a transcriptional repressor belonging to the MeCP1 histone deacetylase complex. Nat Genet. 1999;23:58–61.
Ghoshal K, Majumder S, Datta J, Motiwala T, Bai S, Sharma S, et al. Role of Human Ribosomal RNA (rRNA) Promoter Methylation and of Methyl-CpG-binding Protein MBD2 in the Suppression of rRNA Gene Expression. J Biol Chem. 2004;279:6783–93.
Valentijn L, Koppen A, van Asperen R, Root H, Haneveld F, Versteeg R. Inhibition of a New Differentiation Pathway in Neuroblastoma by Copy Number Defects of N-myc, Cdc42, and nm23 Genes. Cancer Res. 2005;65:3136–45.
Boon K, Caron H, van Asperen R, Valentijn L, Hermus M, van Sluis P, et al. N-myc enhances the expression of a large set of genes functioning in ribosome biogenesis and protein synthesis. EMBO J. 2001;20:1383–93.
Johnston I, Spence H, Winnie J, McGarry L, Vass J, Meagher L, et al. Regulation of a multigenic invasion programme by the transcription factor, AP-1: re-expression of a down-regulated gene, TSC-36, inhibits invasion. Oncogene. 2000;19:5348–58.
Martin B, Brenneman R, Becker K, Gucek M, Cole R, Maudsley S. iTRAQ Analysis of Complex Proteome Alterations in 3xTgAD Alzheimer's Mice: Understanding the Interface between Physiology and Disease. Plos One. 2008;3:e2750.
Lyahyai J, Serrano C, Ranera B, Badiola JJ, Zaragoza P, Martin-Burriel I. Effect of scrapie on the stability of housekeeping genes. Anim Biotechnol. 2010;21:1–13.
Yakubu DP, Mostyn A, Hyatt MA, Kurlak LO, Budge H, Stephenson T, et al. Ontogeny and nutritional programming of mitochondrial proteins in the ovine kidney, liver and lung. Reproduction. 2007;134(6):823–30.
Gupta S, Alfaidy N, Holloway AC, Whittle WL, Lye SJ, Gibb W, et al. Effects of cortisol and oestradiol on hepatic 11beta-hydroxysteroid dehydrogenase type 1 and glucocorticoid receptor proteins in late-gestation sheep fetus. J Endocrinol. 2003;176(2):175–84.
Bonefeld K, Moller S. Insulin-like growth factor-I and the liver. Liver Int: Off J Int Assoc Study Liver. 2011;31:911–9.
Marcolongo P, Fulceri R, Gamberucci A, Czegle I, Banhegyi G, Benedetti A. Multiple roles of glucose-6-phosphatases in pathophysiology: state of the art and future trends. Biochim Biophys Acta. 1830;2013:2608–18.
Parent R, Beretta L. Translational control plays a prominent role in the hepatocytic differentiation of HepaRG liver progenitor cells. Genome Biol. 2008;9:R19.
Snykers S, Henkens T, De Rop E, Vinken M, Fraczek J, De Kock J, et al. Role of epigenetics in liver-specific gene transcription, hepatocyte differentiation and stem cell reprogrammation. J Hepatol. 2009;51:187–211.
Schrem H, Klempnauer J, Borlak J. Liver-enriched transcription factors in liver function and development. Part I: the hepatocyte nuclear factor network and liver-specific gene expression. Pharmacol Rev. 2002;54:129–58.
Gentili S, Morrison JL, McMillen IC. Intrauterine Growth Restriction and Differential Patterns of Hepatic Growth and Expression of IGF1, PCK2, and HSDL1 mRNA in the Sheep Fetus in Late Gestation. Biol Reprod. 2009;80(6):1121–7.
Passmore M, Nataatmadja M, Fraser JF. Selection of reference genes for normalisation of real-time RT-PCR in brain-stem death injury in Ovis aries. BMC Mol Biol. 2009;10:72.
Limesand SW, Rozance PJ, Smith D, Hay Jr WW. Increased insulin sensitivity and maintenance of glucose utilization rates in fetal sheep with placental insufficiency and intrauterine growth restriction. Am J Physiol Endocrinol Metab. 2007;293:E1716–1725.
Williams Jr RC, Shah C, Sackett D. Separation of tubulin isoforms by isoelectric focusing in immobilized pH gradient gels. Anal Biochem. 1999;275:265–7.
Yiallourides M, Sebert SP, Wilson V, Sharkey D, Rhind SM, Symonds ME, Budge H: The differential effects of the timing of maternal nutrient restriction in the ovine placenta on glucocorticoid sensitivity, uncoupling protein 2, peroxisome proliferator-activated receptor-gamma and cell proliferation. Reproduction 2009, 138:601-608. table 3: [26, 59, 64–66]
Huang H, Cheng F, Wang R, Zhang D, Yang L. Evaluation of Four Endogenous Reference Genes and Their Real-Time PCR Assays for Common Wheat Quantification in GMOs Detection. PLoS One. 2013;8:e75850.
He JQ, Sandford AJ, Wang IM, Stepaniants S, Knight DA, Kicic A, et al. Selection of housekeeping genes for real-time PCR in atopic human bronchial epithelial cells. Eur Respiratory J: Off J Eur Soc Clin Respiratory Physiol. 2008;32:755–62.
geNorm manual. [http://medgen.ugent.be/~jvdesomp/genorm/geNorm_manual.pdf].
We would like to thank Mr. Steven Sievers, Mrs. Karen Schellong and Mr. Thomas Ziska for their valuable assistance with the real time PCR experiments and Dr. Rebecca Rancourt for her critical reading of the manuscript. This research was supported by DFG (BR2925/3-1, 3–2 and PL241/8-2), the Canadian Institutes of Health Research, The Raine Medical Research Foundation of Western Australia, the Australian National Health and Medical Research Council (303261), Women and Infants Research Foundation of Western Australia and the Child Health Research Foundation of Western Australia Inc.
The authors declare that they have no competing interests.
HX and MB contributed equally to this work. carried out the molecular genetic studies and participated in the sequence alignment. HX, MB, and LE drafted the manuscript. MB and DS carried out the bioinformatics analyses. DS helped to draft the manuscript. SL participated in the study design, carried out the animal experiments and helped to draft the manuscript. JC, JN and TB conceived of the study, and participated in its design and coordination and drafted the manuscript. JWD, AP and WH helped to draft the manuscript and helped with data interpretation. All authors read and approved the final manuscript.
Huaisheng Xu and Massimo Bionaz are first authors.
Comparison of relative quantities of 6 ICGs across gestation. Kruskal-Wallis One Way Analysis of Variance on Ranks: different numbers indicate significant differences in RQ of ICGs (clour coded) across gestation with p<0.05. Relative Quantity of gene expression = efficiency^(ct min - ct sample)/ RQmin.
About this article
Cite this article
Xu, H., Bionaz, M., Sloboda, D.M. et al. The dilution effect and the importance of selecting the right internal control genes for RT-qPCR: a paradigmatic approach in fetal sheep. BMC Res Notes 8, 58 (2015). https://doi.org/10.1186/s13104-015-0973-7