BACKGROUND
the main quantitative method for the study of gene expression is reverse transcription real-time quantitative pcr , which is considered a highly sensitive technique.
in qpcr, the amount of amplified product is monitored during the course of the reaction by measuring the fluorescence during the annealing phase of each amplification cycle.
fluorescence is produced by dyes or probes which bind to dna, and so it is proportional to the amount of synthesized product.
the dna intercalating dye sybr green i is one of the most widely applied systems, since the fluorescence readings can be obtained from any pcr amplicon, irrespective of its sequence
<cit> .
two types of qpcr can be performed: the expression levels of the genes can represent either an absolute quantification that relates the pcr signal to the initial copy number using a calibration curve or, as in our work, a relative quantification which measures the relative change in rna expression level.
a number of technical parameters such as rna and cdna quality, primer specificity, pcr efficiency and the genes used for normalization heavily condition the quality of qpcr results.
despite the widespread popularity of qpcr, there is a worrying lack of consensus on how it should be performed and how its results should be analyzed.
the publication of the miqe guidelines
<cit>  represented a landmark towards qpcr standardization, but not only are there many publications which still ignore the miqe guidelines, but also new controversies have arisen which require further discussion.
due to the quantitative nature of qpcr, an appropriate normalization method is critical to achieve reliable results.
the purpose of normalization is to remove sampling noise  in order to estimate gene expression accurately
<cit> .
ideally, reference genes used for this purpose should show the same level of expression in all cells and tissues, and remain stable under different experimental conditions.
as pointed out in several publications, there is no universal reference gene, and housekeeping gene expression can vary considerably
<cit> , the best reference gene probably varying in the same species according to the tissue and the experimental conditions
<cit> .
so, as mentioned in the miqe guidelines, normalization against a single reference gene is not recommended unless a clear evidence of its invariant expression is described for the specific experimental conditions of the study.
the optimal number and choice of reference genes should be experimentally determined
<cit> , yet many publications employ a single normalization gene without appropriate validation.
several methods and software have been described to determine the optimum reference genes, however which method is the most suitable has still not been addressed.
four reference gene determination methods are commonly used in qpcr studies: the comparative delta-ct method
<cit> , bestkeeper
<cit> , genorm
<cit>  and normfinder
<cit> .
gene expression stability is evaluated differently in each of the four methods.
briefly, the comparative delta-ct method calculates the stability of each gene by obtaining the standard deviation of cq differences  within each sample for each pairwise comparison with the other genes and averaging them.
normfinder takes into account both intra-group and inter-group gene variation to evaluate its stability.
bestkeeper ranks the genes according to the standard deviation  of their cqs, but the output includes more information, for example the coefficient of variation , which was proposed as a validation method for the results offered by normfinder and genorm
<cit> .
genorm determines the pairwise standard deviation of cq values of all genes, and then excludes the one with the lowest stability, repeating the process until only two genes remain, which are then considered the most stable ones.
another topic, which has recently focused the attention of specialist on this ground, is the kinetics of qpcr and the efficiency determination associated to it.
traditionally, standard curves have been the gold standard to calculate qpcr efficiency.
however, pipetting errors or poorly calibrated pipettes can greatly affect the accuracy of the standard curves due to the cumulative nature of error
<cit> .
also, cdna may include pcr inhibitors which diminish the efficiency of the qpcr reaction.
these inhibitors often remain in the samples from steps prior to qpcr amplification.
the dilution steps involved in standard curve construction, which also dilute inhibitors, might lead to efficiency overestimation
<cit> .
this can be easily confirmed by the existence of efficiencies above 100% and the usual practice of accepting a pair of primers as valid if its efficiency is between 90-110%.
theoretically, it is impossible to obtain qpcr efficiencies above 100%.
more recently, several mathematical models have been published describing the kinetics of the qpcr reaction and trying to estimate qpcr efficiency from a single reaction.
many different models have been proposed, ranging from exponential
<cit>  to logistic ones employing up to five parameters
<cit> ; even more complex models, which take into account the efficiency of each of the steps of the qpcr reaction, have been tackled
<cit> .
here we analyzed four methods which allow an easy determination of efficiency for each reaction and amplicon: i) linregpcr
<cit> , ii) lreanalyzer
<cit> , iii) dart
<cit>  and iv) pcr-miner
<cit> , all publicly available and implemented in user-friendly software or online applications.
marine flatfish represent a valuable group of teleosts because of their highly appreciated white flesh
<cit> .
turbot is a marine flatfish species with a notable aquaculture projection in europe.
it is predicted that by  <dig> its production will duplicate that of  <dig>  .
also, since turbot was introduced in china in  <dig>  the farming industry of this species has developed into one of the main mariculture industries with a production of  <dig> tons per year
<cit> .
the main trait targets for genetic breeding programs in this species are growth rate, sex ratio and disease resistance
<cit> .
turbot shows one of the largest sex-dependant size dimorphism in marine aquaculture
<cit> : females outgrow males by 50% when they reach commercial size.
some studies have demonstrated a zz/zw system in turbot
<cit>  and identified the main sex determining region in linkage group  5
<cit> , but these authors also suggest the existence of other minor genetic and environmental factors, for example temperature, which might affect sex determination.
however, expression analyses have only been carried out in immune tissues so far
.
reference genes for qpcr have been characterized in different tissues of turbot
<cit>  and in other flatfish
<cit> , but gonads have not been included in these studies.
in this study, we evaluated the main factors which might compromise qpcr results, reference gene choice and qpcr efficiency determination, using gonads of turbot reared at different temperatures and along the development process.
our results suggest that for research purposes, normfinder and linregpcr implement the best approaches for reference gene selection and efficiency determination, respectively, ant that ubq and rps <dig> would be the best reference genes for the normalization of gonad development in turbot from  <dig> up to 135 days post fertilization.
to our knowledge, this is the first qpcr evaluation in turbot gonads and no similar studies have been carried out in fish to date.
our approach, although applied in a particular tissue in turbot could be used as a guideline for qpcr development in other tissues or species.
RESULTS
amplification
amplification of each reference gene in  <dig> samples  produced a  <dig> cq values dataset.
samples with missing cq values or inconsistencies between replicates  in any of the reference genes were removed from the analysis.
after averaging duplicates a total of  <dig> samples were kept  and we obtained descriptive statistics and kolmogorov-smirnov tests to check for normality for each of the assayed genes .
a single amplification product for each primer pair was confirmed by a single peak in the melting curve analysis and also by pcr product sequencing.table  <dig>
descriptive statistics of the reference genes cq values
actb
b2m
gapdh
rpl17
rps4
ubq
number of samples , mean, standard deviation , minimum cq value , maximum cq value  and p value of the kolmogorov-smirnov test  are shown for each candidate reference gene.
actb showed the highest expression , amplification being more than two cycles earlier than any other gene.
on the other side, rps <dig> showed the lowest expression .
ubq standard deviation  was the lowest  while b2m presented the largest variation between cq values .
also, reference gene cq distributions were normal in every case but that of b2m .
according to the experimental design, samples were divided in groups according to fish age in days post fertilization  and rearing temperature.
this produced a total of  <dig> groups , with a minimum of six samples per group and a maximum of ten.
a boxplot of all the groups for the six reference genes can be observed in figure
<dig>  we also considered grouping our samples by degree-days, however, since groups remained basically the same , which did not alter neither the results nor the discussion, we decided to name the groups by their age and rearing temperature since we considered it clearer .figure  <dig>
reference gene cq value distributions.
boxplots of the cq values in each experimental group  for each of the six reference genes.
each group is named with a number, which indicates age in days post fertilization, and either “high”, “normal” or “low” which indicates rearing temperature.
a similar dataset of cq values for six sex differentiation related genes was obtained and their descriptive statistics are presented in table
<dig>  these genes are involved in gonad differentiation and were used to check normalization and efficiency correction effects.table  <dig>
descriptive statistics of the genes involved in gonad differentiation cq values
cyp19a
amh
sox19
sox9
vasa
sox17
number of samples , mean, standard deviation , minimum cq value , maximum cq value  and p value of the kolmogorov-smirnov test  are shown for each target gene.
analysis of the reference genes
we analyzed the  <dig> cq values obtained for each of the reference genes with comparative delta ct method, bestkeeper, normfinder and genorm.
for each method and gene a ranking of stability values is shown with the most stable gene at the top and the least stable at the bottom .
due to the importance of gene-to-gene correlations for comparative delta-ct method and genorm, correlations were graphically represented .
finally, the average intergroup and intragroup variation for each gene is shown in table
<dig> as reported by normfinder.
groups were formed as specified above according to fish age and rearing temperature.table  <dig>
stability rankings obtained with the different reference gene determination methods
1
2
3
4
5
6
stability values obtained by each method are shown in parenthesis for each candidate reference gene.
both standard deviation  and coefficient of variation  rankings are shown for bestkeeper.
the genes are ranked from most stable  to least stable .
correlation between reference genes.
legend: correlation between reference genes cq values.
the highest correlations are colored in red, medium correlations in green and the lowest in yellow.
correlation coefficient  values are shown, p value <  <dig> .
intra-group and inter-group variation estimates by normfinder
actb
b2m
gapdh
rpl17
rps4
ubq
intra-group
inter-group
average intra-group and inter-group cq variation estimates obtained by normfinder.
groups were constituted by fish of the same age and rearing temperature.
in our analysis, ubq appears ranked first by bestkeeper, genorm and comparative delta-ct method, and ranked second by normfinder.
b2m, a frequently used reference gene, is consistently ranked the last by all methods and does not show a strong correlation with any of the other genes  .
this gene cq distribution also deviated from normality.
the inconsistency of b2m as reference gene has also been previously reported in human tissues
<cit> .
the ranking between these two extremes varies depending on the method.
bestkeeper ranked ubq as the most stable gene with  <dig> , a value above the recommended cutoff of 1
<cit> .
however, considering that our dataset includes samples coming from different tissues , experimental conditions and development stages, a low standard deviation was not expected.
ubq also shows the lowest cv and a high correlation with rps <dig> .
actb is ranked as the second best reference gene by bestkeeper , however, it does not show a high correlation with any of the other genes and it is ranked 4th by other methods .
a possible explanation for this discrepancy might be the high expression shown by actb  .
this renders a lower error when measuring the fluorescence values and a lower copy number difference between samples, which does not imply higher stability between the different experimental conditions.
this should be taken into account when choosing a reference gene since they are usually highly expressed genes.
the most important characteristic of a reference gene is that its sample-to-sample variation must be representative of the technical error produced by the sampling, extraction and retrotranscription steps in order to reduce target gene cq values error.
bestkeeper ranks rpl <dig> third , rps <dig> fourth  and gapdh fifth  .
an interesting fact is that only rps <dig> shows high correlations with other genes  .
like actb but at the other extreme, rps <dig> shows the lowest amplification cycle , which can result in higher technical error.
so, according to bestkeeper, ubq is the most stable gene, followed by actb and rpl <dig>  as mentioned, probably actb is not a good reference gene and it seems that rps <dig> might deserve a better ranking given its correlation with other genes.
normfinder, which assesses inter-group variation  in order to discard regulated genes, ranked rps <dig> as the most stable gene , ubq appears second  and rpl <dig> third  .
rps <dig> showed the lowest inter-group variation , which explains its ranking .
gapdh is clearly pointed as less stable than rps <dig> and rpl <dig> by normfinder, with an overall score of  <dig>  and an inter-group variation of  <dig> .
surprisingly, normfinder suggests actb and rps <dig> as the most stable two gene combination.
actb is ranked 4th by normfinder  .
the information provided by bestkeeper and normfinder should be enough to decide among the most stable genes, to say, those with an acceptable low level of overall variation, low inter-group variation and high correlation between them.
the last two methods, comparative delta-ct method and genorm, both follow pairwise approaches but with different procedures and outcomes; the first one ranks the genes following an average pairwise sd while the second follows a progressive exclusion of the least stable gene after pairwise comparison.
both methods agree with the results obtained by normfinder.
genorm recommends the couple of reference genes ubq/rps <dig> with a value of  <dig>  , and points that adding another reference gene  would not improve normalization .
pairwise comparison methods tend to select those genes with the highest degree of similarity across the sample set, implying that the candidates with minimal expression variation do not necessarily become top ranked
<cit> .
while both approaches are based in pairwise comparisons, the progressive exclusion of genes by genorm increases the tendency to select the most correlated genes.
since we have obtained inter-group variation estimates by normfinder  which points toward rps <dig>  ubq, rpl <dig> and possibly actb not being differentially expressed between groups, and also due to the unexpected presence of actb in the best normalization indexes calculated by normfinder, we checked how the use of different normalization factors, geometric mean of different reference genes, affected six sex-differentiation related genes  .
we normalized the samples by , , , , ,  and .
we checked the intragroup and intergroup standard deviation for each of the six target genes and, since the samples were genetically sexed using the smausc-e <dig> marker according to
<cit> , we also checked the standard deviation of male and female groups .
this way, we can have an independent measure out of the fish age/rearing temperature groups we have used to check the stability of the reference genes.
interestingly, ubq and rps <dig> seem to behave differently.
while rps <dig> renders lower sd values for fish age/rearing temperature groups, ubq normalizes male and female groups better.
however, the lowest sd estimates were obtained when both  were used for normalization, except for average intergroup sd where rps <dig> alone performed better.
the addition of rpl <dig> or actb in the index for normalization did not yield lower sd estimates.
the use of just  for normalization is in agreement with the results of genorm and also with the rankings produced by normfinder and comparative delta-ct method.
the use of  as suggested by normfinder does not perform better.table  <dig>
standard deviation for target genes when normalized by different gene combinations
ubq
rps4
ubq + rsp4
actb + rps4
ubq + rps4 + rpl17
ubq + rps4 + actb
ubq + rps4 + rpl17 + actb
intragroup and intergroup normalized cq standard deviations  averaging the results for the six target genes are shown for fish age + rearing temperature groups when normalized by different candidate reference gene combinations.
standard deviations  for males and females when normalized by the same combinations are also shown.
to assess the robustness of each method, we repeated the stability calculations in fifty subsets of the samples ,  <dig> subsets include  <dig> samples of each experimental group  and another  <dig> include  <dig> samples of each experimental group .
we evaluated the  <dig> subsets together since the results show similar trends both with three and two samples per group.
since in many studies three genes are used for normalization, we compared not only the whole ranking but also the top <dig> genes.
the most robust method is clearly bestkeeper sd, which renders an identical ranking as that obtained with the whole data set for a 44% of the subsets and, in 40% of the remaining subsets, it ranks the top three genes correctly .
normfinder selected the same top <dig> genes also in 88% of the subsets, however the rank order was altered most of the times.
on the contrary, the pairwise approaches showed a higher degree of variation, the top <dig> genes were different from those in the full dataset in 60% of the subsets for comparative delta ct method and in 66% for genorm.table  <dig>
robustness of the gene stability determination method
similarity of  <dig> subsets stability rankings by each method and the ranking obtained with the whole dataset.
<dig> subsets are formed by  <dig> samples per group  and another  <dig> subsets have  <dig> samples per group.
efficiency determination analysis
we obtained mean gene efficiencies by lingregpcr, lreanalyzer, dart and pcr-miner for each primer pair  and correlations between mean efficiencies by each method for each gene .
there is around a 10% difference between the efficiencies calculated by linear fit methods  and non linear fit models , meaning that exponential methods might be underestimating efficiency or non linear methods overestimating it .
two lreanalyzer efficiency estimates are over 100% , which is theoretically impossible for a pcr reaction, so lreanalyzer is likely overestimating qpcr efficiency.
however, despite this 10% efficiency difference, mean efficiencies calculated by the four methods are correlated for each gene, indicating that although they are using different algorithms they are rendering similar relative results.
best correlation coefficient and p value are observed between linregpcr and pcr-miner, which might be highlighting the importance of baseline correction since both methods use iterative approaches instead of relying on a fluorescence correction based on the average fluorescence of the first qpcr cycles.
also, linregpcr and pcr-miner include several functions to remove outliers, so filtering the reactions before efficiency calculation might also be important to obtain more precise efficiency estimations.table  <dig>
efficiency values for each gene with each efficiency determination method
mean efficiency values for each reference gene with lingregpcr, lreanalyzer, dart and pcr-miner.
correlation between efficiency determination methods
pearson correlation coefficients and p values  for mean gene efficiencies with each of the four efficiency determination methods are shown.
normalization and efficiency correction on target genes
six target genes  involved in gonad differentiation were efficiency corrected and normalized by four different combinations of efficiency determination methods and reference gene combinations .
for each combination, first, efficiency correction was performed on every cq value of both reference and target genes.
afterwards, each target efficiency-corrected cq value was normalized by the reference gene/s efficiency-corrected cq values, obtaining efficiency-corrected delta cq values.
we computed mean efficiency-corrected delta cq values and standard deviations for the three temperature groups   and also for males and females .
two different patterns are shown in the tables, one caused by normalization and the other by efficiency correction.
a higher standard deviation is obtained in most of the b2m normalized dataset compared to the ubq + rps <dig> normalized ones, which is expected when a gene is not stable.
however, this is not true for the amh normal temperature group neither for the gene vasa, suggesting some type of co-regulation.
the other trend is observed when comparing the linregpcr efficiency corrected datasets with the pcr-miner corrected ones.
pcr-miner produces higher mean cqs  increasing the difference between groups.table  <dig>
efficiency-corrected delta cqs by temperature group with each efficiency + reference gene combination
mean efficiency-corrected delta cqs and sd values for the three rearing temperatures : high, normal and low; in the four datasets produced after efficiency correction with linregpcr or pcr-miner and later normalization with ubq + rps <dig> or b2m.
efficiency-corrected delta cqs by sex group with each efficiency + reference gene combination
mean efficiency-corrected delta cqs and standard deviation  values for males and females in the four datasets produced after efficiency correction with linregpcr or pcr-miner and later normalization with ubq + rps <dig> or b2m.
furthermore, the use of a gene which presents systematic differences between groups for normalization can lead to changes in the mean cq values of some genes.
for example, amh gene expression in each temperature group is severely affected by normalization with b2m, varying from  <dig>  to - <dig>  at high temperature , from  <dig>  to  <dig>  at normal temperature and from - <dig>  to - <dig>  at low temperature .
discussion
reference gene analysis
the four methods commonly used to check the stability of reference genes, comparative delta-ct method, normfinder, bestkeeper and genorm, represent viable strategies, although none of them is currently considered the best one and some problems can arise in certain experimental scenarios.
the bestkeeper method is apparently the “common sense” solution to measure stability since standard deviation is a direct measure of variation.
however, a gene might show a low standard deviation but still not be a good reference gene if its variation does not reflect the errors produced by sampling, rna extraction and retrotranscription steps.
this problem could be circumvented by analyzing the correlations between genes, assuming that the reference genes are not co-regulated.
this means that sampling point differences  affecting one of the genes should not affect the others, and so the correlations between them would reflect the inter-sample variation produced by the sample processing steps and not by co-regulation due to the experimental conditions.
nevertheless, it is risky to assume that genes are not co-regulated because this cannot be easily demonstrated.
the genorm and the comparative delta-ct method approaches present the same problem but in addition these methods rank genes mainly by their correlations, to say, genorm establishes the most stable genes by assuming “that the control reference genes are not co-regulated”
<cit> , and the same happens to the comparative delta-ct method which follows a very similar approach.
as a consequence, two co-regulated genes could fully spoil the analysis leading to wrong reference genes.
finally, normfinder is not affected by the co-regulated gene drawback since it takes into account intergroup variation , which should be as lower as possible for a good reference gene; however, similarly to bestkeeper, a low overall intergroup and intragroup variation does not necessarily mean that it is a good reference gene.
the advantages and disadvantages of each strategy should be taken into account when analyzing putative reference genes according to the experimental scenario.
normfinder and genorm are the most extended methodologies to find the optimum reference genes.
in many cases, normfinder and genorm algorithms render very similar results, however, discrepancies between the output of normfinder and genorm have been previously described
<cit> .
in these works, the cv has been used to decide which genes should be used for normalization, confirming normfinder results in every case.
while normfinder results are non-biased, genorm stepwise exclusion can lead to awkward results by selecting reference genes which in fact are not the most stable.
normfinder, bestkeeper and comparative delta-ct stability method results have also been reported to be more consistent among them than with those of genorm
<cit> , although in other study bestkeeper was reported as the least consistent method
<cit> .
our results with the whole dataset support the high consistency between normfinder and comparative delta-ct method, while bestkeeper results seem to be the least consistent and only correlation values between reference genes seem to suggest a similar ranking.
however, when working with different subsets which include a lower number of samples, the pairwise approaches results vary significantly between subsets.
this lack of robustness has been described previously: it was shown that the exclusion of a single sample could change the status of one gene from unstable to 2nd most stable gene by genorm
<cit> .
genorm lack of robustness can most likely be explained by the removal of the least correlated gene by pairwise comparison with all the others until only two genes are left, which can lead to stable genes being removed of the analysis early on.
robustness is a critical parameter.
since experiments are budget limited, it is important to be able to determine correct reference genes with a low number of qpcr reactions.
bestkeeper and normfinder appear to be more robust than comparative delta-ct method and genorm in our study.
there is not a method to check how much normalization has improved our gene expression data.
in principle, a reduction in the cq variability of the gene of interest should be expected, however the highest reduction of this variability would also occur if the gene of interest and the reference gene are co-regulated, so this is a risky strategy.
an example of this is observed in vasa  standard deviation after normalization, obtaining a lower sd when normalized by a clearly not stable gene  than when normalized by ubq + rps <dig>  suggesting co-regulation between vasa and b2m.
the same applies to detecting significant/non-significant results depending on the reference gene used for normalization.
this is only useful to stress the importance of choosing a good reference gene, not to choose between one or another since a co-regulated reference gene would lead to non-significant results even if there are expression differences between groups.
the fact that there is not a post-control which enables us to check if we have chosen the correct reference gene makes the choice even more critical.
every experiment and dataset is different, so the analysis has to be done carefully.
given the huge importance of normalization and its great impact on the conclusions, it would be recommended to analyze each case separately, paying attention to details; using any method as a black box can lead up to low confident results.
several studies have solved the disagreement between the four methods by ranking them according to the geometric mean of the four ranking numbers for each gene, the lower the mean a gene gets the most stable it is
<cit> .
however, attributing the same weight to every method is arguable, especially because some of these methods include redundant information.
this is a practical option without any biological meaning.
if the four methods disagree, we recommend instead relying on the ranking provided by normfinder, while ignoring its suggested combination, supported by descriptive statistics like mean, standard deviation and correlations, information offered by bestkeeper or any common statistical package.
this approach would enable to assess the two most important and complementary issues: absence of inter-group variation and correlation between reference genes.
this approach does not make any previous assumption and has proven to be robust when only a few samples are assayed.
to our knowledge this is the first experiment to analyze the stability of reference genes during the gonad development in fish.
even in mature organs, there is only one study carried out in zebrafish were testis and ovaries were analyzed separately
<cit> .
however, studies have been carried out in other organs.
the stability of several genes was studied in the liver, spleen, kidney, heart, brain, gill and muscle of turbot subjected to edwardsiella tarda infection
<cit> .
gene stability was checked before infection in all the organs together by normfinder and genorm.
in that study, out of eight genes, normfinder ranked rpsd as the most stable one, followed by actb, rpl <dig>  b2m and gapdh among those genes shared with our study, although primer pairs for gapdh were different.
we tried to develop primers for rpsd but they were discarded due to late amplification cycle in gonad .
ubq and rps <dig>  were not assayed in that work.
gapdh, which has been classically used as a reference gene but recently classified as unstable in several studies
<cit> , performed badly in both our study and
<cit> .
however, gapdh is ranked as the most stable gene in heart and liver in
<cit> , which emphasizes the importance of checking reference gene stability in each study separately, reference genes cannot be “exported”.
there are two gapdh isoforms in diploid teleost fish as a result of the fish-specific genome duplication event, however the same variant has been analyzed in both studies .
there are two similar qpcr studies carried out in flatfish.
the first one studied six reference genes during hippoglossus hippoglossus development in sixteen different tissues using bestkeeper, normfinder and genorm
<cit> .
gonads were not included.
they assayed actb which was found as one of the least stable genes.
the most stable genes found were ef1a <dig> and rpl <dig>  the second study was carried out during larval development in solea senegalensis and hippoglossus hippoglossus.
the stability of twelve genes, including ubq, rps <dig>  actb and gapdh- <dig>  was checked by genorm and normfinder
<cit> .
the combined stability index of the two species ranked ubq, rps <dig> and eef1a <dig> as the three best normalization genes by both methods, while actb appears in 5th and 6th place, respectively.
interestingly, in solea senegalensis, gapdh- <dig> appears ranked 3rd and 1st by genorm and normfinder, respectively.
efficiency determination analysis
efficiency determination is an essential step in qpcr.
constant amplification efficiency in all compared samples is a very important criterion for reliable comparison between samples.
it is also crucial for an accurate quantification of gene expression.
ideally, the efficiency of an assay should be 100%, which means that during the logarithmic phase of the reaction the pcr product is doubling each cycle.
each of the four tested efficiency determination methods differ in their baseline fluorescence determination, type of fit to the log-linear phase of the qpcr reaction, and preprocessing steps to remove outliers.
dart
<cit>  is based on a linear regression of the exponential phase of the qpcr reaction.
baseline subtraction is determined by fitting a saturation function to the first 2- <dig> cycles of the qpcr reaction.
then, a linear regression is performed in a 10-fold range around the middle point of the exponential phase, which is calculated using the maximum fluorescence and standard deviation of the fluorescence in the first  <dig> cycles.
linregpcr
<cit>  is also based on a linear regression fit to the log-linear phase of the amplification curve.
linregpcr determines baseline fluorescence through an iterative algorithm to get the best fit of the linear regression to 4- <dig> points in the log-linear phase of the reaction.
then, after baseline subtraction, these points are used for efficiency determination.
pcr-miner
<cit>  uses a non-linear regression fit.
as linregpcr, baseline fluorescence is determined by an iterative fit to a four-parameter logistic model which also determines the exponential phase of the reaction.
then, a three-parameter exponential model is fitted to the exponential phase to determine efficiency.
finally, lreanalyzer
<cit>  uses a sigmoidal fit approach.
baseline subtraction is determined by averaging 6- <dig> cycles fluorescence values.
then, efficiency estimates are calculated for each cycle of the qpcr reaction.
an lre window is selected with those cycle efficiencies which fit to a linear regression.
finally, a derivative of the boltzmann sigmoidal function is used for the calculation of the maximum efficiency of the reaction.
logistic models are pure empirical models not designed to be kinetically realistic
<cit>  and rely purely in their good fit to the real-time pcr curve.
as previously reported by other study, qpcr curves are not symmetric since they do not have the same curvature at both sides of the inflection point, implying the existence of two or more different mechanisms affecting the efficiency of the reaction
<cit>  and so, making this good-fit models hardly reliable.
furthermore, in some reactions sybr green depletion might be the main mechanism leading to the plateau phase of the curve
<cit> .
while sybr green has an impact in the visualization of the real-time reaction, it does not have a connection with the kinetics of the pcr reaction.
a recent qpcr study has tried a new approach to assess qpcr efficiency, defining the global efficiency as the sum of denaturing efficiency, annealing efficiency, polymerase binding efficiency and elongation efficiency
<cit> .
the polymerase binding efficiency and the elongation efficiency can be constant provided that there is an excess of polymerase and a long elongation time.
however, the denaturing efficiency is constantly decreasing each cycle at the same rate due to thermal damage in both the dna and the polymerase.
the annealing efficiency is also decreasing and depends on the proportion of ssdna bound to the primers during this step and total ssdna present.
some ssdna chains might bind to its complementary strand instead of to the primers
<cit> .
this efficiency varies from cycle to cycle.
this theoretical study was validated in
<cit> .
so, at first, a constant efficiency should not be assumed.
still, the qpcr curve shows a large exponential component, since, as confirmed by a previous study, in most cases the best fit to the log-linear region of the qpcr reaction is exponential
<cit> , suggesting that before and at the log-linear region the qpcr efficiency reduction is low.
a recent study analyzed all publicly available efficiency determination methods
<cit>  in a large dataset, included four-point 10-fold dilution series, which allows calculation of the bias of each method.
similarly to our results, they report linregpcr and dart to produce an underestimation of efficiency and pcr-miner and lre analyzer an overestimation.
they also analyze different parameters and find linregpcr and pcr-miner amongst the best methods for most of the evaluated characteristics, for example precision and resolution, performing better than lreanalyzer and dart.
the reader is encouraged to consult
<cit>  to learn more about the different efficiency determination methods, their characteristics and performance differences.
this study
<cit>  is the most complete on qpcr efficiency determination methods done so far.
both linregpcr and pcr-miner performed similarly and produced highly correlated efficiency estimates in our study.
the main difference was that while linregpcr underestimates efficiency, pcr-miner overestimated it.
knowing this, linregpcr is probably the best choice for the average qpcr researcher since it will not produce erroneous significant differences between groups  or an overestimation of the fold change.
however, linregpcr might not be the best option for clinical purposes, where the method of choice should be considered depending on the consequences of a false positive/overestimation or a false negative/underestimation.
there are a good number of alternative efficiency estimation algorithms, however they are implemented as extensions to the open source statistical programming environment r
<cit>  and probably not available for most researchers, so they have not been analyzed here.
still, linregpcr and pcr-miner perform as well or better than all other methods, as shown in
<cit> .
as a final remark, although
<cit>  clearly improved our understanding of the different efficiency determination methods available, currently there is not a clear best method for estimating qpcr efficiencies.
however, the publication of a theoretical study based on the pcr kinetics which defines the overall pcr efficiency as the product of the efficiency of each of the separate steps
<cit>  and its experimental validation
<cit>  is a good step towards finding a biologically meaningful solution.
similar approaches will be likely applied in more studies in the near future provided they are implemented in appropriate user friendly softwares for the whole research community.
normalization and efficiency correction on target genes
the effects of normalization with wrong genes are important; a high standard deviation will produce higher p values and, so, possibly lead to missing biologically relevant differences.
even worse, the use of a regulated gene  for normalization, will lead to changes in gene values which can end in misguided results.
the effect of efficiency correction, though not so dramatic, is also important since it can lead to overestimation  of differences between groups.
CONCLUSIONS
we found the ranking produced by normfinder method as the most reliable one to choose reference genes for qpcr analysis when results differ between gene stability determination methods.
normfinder information should be complemented by the descriptive statistics offered by bestkeeper, especially the correlation coefficient.
accordingly, we found that ubq and rps <dig> should be used as reference genes to study turbot gonad development from  <dig> up to 135 days post fertilization.
we found pair-wise methods to be less robust than normfinder and bestkeeper and also the suggested normfinder two genes combination not reliable.
we also recommend the use of linregpcr for efficiency determination for research purposes, however, efficiency determination is still a matter of discussion and probably new improved models will be published in the upcoming years.

 the main quantitative method for the study of gene expression is reverse transcription real-time quantitative pcr, which is considered a highly sensitive technique. in qpcr, the amount of amplified product is monitored during the course of the reaction by measuring the fluorescence during the annealing phase of each amplification cycle. fluorescence is produced by dyes or probes which bind to dna, and so it is proportional to the amount of synthesized product. the dna intercalating dye sybr green i is one of the most widely applied systems, since the fluorescence readings can be obtained from any pcr amplicon, irrespective of its sequence. two types of qpcr can be performed: the expression levels of the genes can represent either an absolute quantification that relates the pcr signal to the initial copy number using a calibration curve or, as in our work, a relative quantification which measures the relative change in rna expression level. a number of technical parameters such as rna and cdna quality, primer specificity, pcr efficiency and the genes used for normalization heavily condition the quality of qpcr results. despite the widespread popularity of qpcr, there is a worrying lack of consensus on how it should be performed and how its results should be analyzed. the publication of the miqe guidelines represented a landmark towards qpcr standardization, but not only are there many publications which still ignore the miqe guidelines, but also new controversies have arisen which require further discussion. due to the quantitative nature of qpcr, an appropriate normalization method is critical to achieve reliable results. the purpose of normalization is to remove sampling noise in order to estimate gene expression accurately. ideally, reference genes used for this purpose should show the same level of expression in all cells and tissues, and remain stable under different experimental conditions. as pointed out in several publications, there is no universal reference gene, and housekeeping gene expression can vary considerably, the best reference gene probably varying in the same species according to the tissue and the experimental conditions. so, as mentioned in the miqe guidelines, normalization against a single reference gene is not recommended unless a clear evidence of its invariant expression is described for the specific experimental conditions of the study. the optimal number and choice of reference genes should be experimentally determined, yet many publications employ a single normalization gene without appropriate validation. several methods and software have been described to determine the optimum reference genes, however which method is the most suitable has still not been addressed. four reference gene determination methods are commonly used in qpcr studies: the comparative delta-ct method, bestkeeper, genorm and normfinder. gene expression stability is evaluated differently in each of the four methods. briefly, the comparative delta-ct method calculates the stability of each gene by obtaining the standard deviation of cq differences within each sample for each pairwise comparison with the other genes and averaging them. normfinder takes into account both intra-group and inter-group gene variation to evaluate its stability. bestkeeper ranks the genes according to the standard deviation of their cqs, but the output includes more information, for example the coefficient of variation, which was proposed as a validation method for the results offered by normfinder and genorm. genorm determines the pairwise standard deviation of cq values of all genes, and then excludes the one with the lowest stability, repeating the process until only two genes remain, which are then considered the most stable ones. another topic, which has recently focused the attention of specialist on this ground, is the kinetics of qpcr and the efficiency determination associated to it. traditionally, standard curves have been the gold standard to calculate qpcr efficiency. however, pipetting errors or poorly calibrated pipettes can greatly affect the accuracy of the standard curves due to the cumulative nature of error. also, cdna may include pcr inhibitors which diminish the efficiency of the qpcr reaction. these inhibitors often remain in the samples from steps prior to qpcr amplification. the dilution steps involved in standard curve construction, which also dilute inhibitors, might lead to efficiency overestimation. this can be easily confirmed by the existence of efficiencies above 100% and the usual practice of accepting a pair of primers as valid if its efficiency is between 90-110%. theoretically, it is impossible to obtain qpcr efficiencies above 100%. more recently, several mathematical models have been published describing the kinetics of the qpcr reaction and trying to estimate qpcr efficiency from a single reaction. many different models have been proposed, ranging from exponential to logistic ones employing up to five parameters; even more complex models, which take into account the efficiency of each of the steps of the qpcr reaction, have been tackled. here we analyzed four methods which allow an easy determination of efficiency for each reaction and amplicon: i) linregpcr, ii) lreanalyzer, iii) dart and iv) pcr-miner, all publicly available and implemented in user-friendly software or online applications. marine flatfish represent a valuable group of teleosts because of their highly appreciated white flesh. turbot is a marine flatfish species with a notable aquaculture projection in europe. it is predicted that by its production will duplicate that of. also, since turbot was introduced in china in the farming industry of this species has developed into one of the main mariculture industries with a production of tons per year. the main trait targets for genetic breeding programs in this species are growth rate, sex ratio and disease resistance. turbot shows one of the largest sex-dependant size dimorphism in marine aquaculture: females outgrow males by 50% when they reach commercial size. some studies have demonstrated a zz/zw system in turbot and identified the main sex determining region in linkage group 5, but these authors also suggest the existence of other minor genetic and environmental factors, for example temperature, which might affect sex determination. however, expression analyses have only been carried out in immune tissues so far. reference genes for qpcr have been characterized in different tissues of turbot and in other flatfish, but gonads have not been included in these studies. in this study, we evaluated the main factors which might compromise qpcr results, reference gene choice and qpcr efficiency determination, using gonads of turbot reared at different temperatures and along the development process. our results suggest that for research purposes, normfinder and linregpcr implement the best approaches for reference gene selection and efficiency determination, respectively, ant that ubq and rps would be the best reference genes for the normalization of gonad development in turbot from up to 135 days post fertilization. to our knowledge, this is the first qpcr evaluation in turbot gonads and no similar studies have been carried out in fish to date. our approach, although applied in a particular tissue in turbot could be used as a guideline for qpcr development in other tissues or species.  amplification amplification of each reference gene in samples produced a cq values dataset. samples with missing cq values or inconsistencies between replicates in any of the reference genes were removed from the analysis. after averaging duplicates a total of samples were kept and we obtained descriptive statistics and kolmogorov-smirnov tests to check for normality for each of the assayed genes. a single amplification product for each primer pair was confirmed by a single peak in the melting curve analysis and also by pcr product sequencing.table descriptive statistics of the reference genes cq values actb b2m gapdh rpl17 rps4 ubq number of samples, mean, standard deviation, minimum cq value, maximum cq value and p value of the kolmogorov-smirnov test are shown for each candidate reference gene. actb showed the highest expression, amplification being more than two cycles earlier than any other gene. on the other side, rps showed the lowest expression. ubq standard deviation was the lowest while b2m presented the largest variation between cq values. also, reference gene cq distributions were normal in every case but that of b2m. according to the experimental design, samples were divided in groups according to fish age in days post fertilization and rearing temperature. this produced a total of groups, with a minimum of six samples per group and a maximum of ten. a boxplot of all the groups for the six reference genes can be observed in figure  we also considered grouping our samples by degree-days, however, since groups remained basically the same, which did not alter neither the results nor the discussion, we decided to name the groups by their age and rearing temperature since we considered it clearer.figure reference gene cq value distributions. boxplots of the cq values in each experimental group for each of the six reference genes. each group is named with a number, which indicates age in days post fertilization, and either “high”, “normal” or “low” which indicates rearing temperature. a similar dataset of cq values for six sex differentiation related genes was obtained and their descriptive statistics are presented in table  these genes are involved in gonad differentiation and were used to check normalization and efficiency correction effects.table descriptive statistics of the genes involved in gonad differentiation cq values cyp19a amh sox19 sox9 vasa sox17 number of samples, mean, standard deviation, minimum cq value, maximum cq value and p value of the kolmogorov-smirnov test are shown for each target gene. analysis of the reference genes we analyzed the cq values obtained for each of the reference genes with comparative delta ct method, bestkeeper, normfinder and genorm. for each method and gene a ranking of stability values is shown with the most stable gene at the top and the least stable at the bottom. due to the importance of gene-to-gene correlations for comparative delta-ct method and genorm, correlations were graphically represented. finally, the average intergroup and intragroup variation for each gene is shown in table  as reported by normfinder. groups were formed as specified above according to fish age and rearing temperature.table stability rankings obtained with the different reference gene determination methods 1 2 3 4 5 6 stability values obtained by each method are shown in parenthesis for each candidate reference gene. both standard deviation and coefficient of variation rankings are shown for bestkeeper. the genes are ranked from most stable to least stable. correlation between reference genes. legend: correlation between reference genes cq values. the highest correlations are colored in red, medium correlations in green and the lowest in yellow. correlation coefficient values are shown, p value < . intra-group and inter-group variation estimates by normfinder actb b2m gapdh rpl17 rps4 ubq intra-group inter-group average intra-group and inter-group cq variation estimates obtained by normfinder. groups were constituted by fish of the same age and rearing temperature. in our analysis, ubq appears ranked first by bestkeeper, genorm and comparative delta-ct method, and ranked second by normfinder. b2m, a frequently used reference gene, is consistently ranked the last by all methods and does not show a strong correlation with any of the other genes. this gene cq distribution also deviated from normality. the inconsistency of b2m as reference gene has also been previously reported in human tissues. the ranking between these two extremes varies depending on the method. bestkeeper ranked ubq as the most stable gene with, a value above the recommended cutoff of 1. however, considering that our dataset includes samples coming from different tissues, experimental conditions and development stages, a low standard deviation was not expected. ubq also shows the lowest cv and a high correlation with rps. actb is ranked as the second best reference gene by bestkeeper, however, it does not show a high correlation with any of the other genes and it is ranked 4th by other methods. a possible explanation for this discrepancy might be the high expression shown by actb. this renders a lower error when measuring the fluorescence values and a lower copy number difference between samples, which does not imply higher stability between the different experimental conditions. this should be taken into account when choosing a reference gene since they are usually highly expressed genes. the most important characteristic of a reference gene is that its sample-to-sample variation must be representative of the technical error produced by the sampling, extraction and retrotranscription steps in order to reduce target gene cq values error. bestkeeper ranks rpl third, rps fourth and gapdh fifth. an interesting fact is that only rps shows high correlations with other genes. like actb but at the other extreme, rps shows the lowest amplification cycle, which can result in higher technical error. so, according to bestkeeper, ubq is the most stable gene, followed by actb and rpl as mentioned, probably actb is not a good reference gene and it seems that rps might deserve a better ranking given its correlation with other genes. normfinder, which assesses inter-group variation in order to discard regulated genes, ranked rps as the most stable gene, ubq appears second and rpl third. rps showed the lowest inter-group variation, which explains its ranking. gapdh is clearly pointed as less stable than rps and rpl by normfinder, with an overall score of and an inter-group variation of. surprisingly, normfinder suggests actb and rps as the most stable two gene combination. actb is ranked 4th by normfinder. the information provided by bestkeeper and normfinder should be enough to decide among the most stable genes, to say, those with an acceptable low level of overall variation, low inter-group variation and high correlation between them. the last two methods, comparative delta-ct method and genorm, both follow pairwise approaches but with different procedures and outcomes; the first one ranks the genes following an average pairwise sd while the second follows a progressive exclusion of the least stable gene after pairwise comparison. both methods agree with the results obtained by normfinder. genorm recommends the couple of reference genes ubq/rps with a value of, and points that adding another reference gene would not improve normalization. pairwise comparison methods tend to select those genes with the highest degree of similarity across the sample set, implying that the candidates with minimal expression variation do not necessarily become top ranked. while both approaches are based in pairwise comparisons, the progressive exclusion of genes by genorm increases the tendency to select the most correlated genes. since we have obtained inter-group variation estimates by normfinder which points toward rps ubq, rpl and possibly actb not being differentially expressed between groups, and also due to the unexpected presence of actb in the best normalization indexes calculated by normfinder, we checked how the use of different normalization factors, geometric mean of different reference genes, affected six sex-differentiation related genes. we normalized the samples by,,,,, and. we checked the intragroup and intergroup standard deviation for each of the six target genes and, since the samples were genetically sexed using the smausc-e marker according to, we also checked the standard deviation of male and female groups. this way, we can have an independent measure out of the fish age/rearing temperature groups we have used to check the stability of the reference genes. interestingly, ubq and rps seem to behave differently. while rps renders lower sd values for fish age/rearing temperature groups, ubq normalizes male and female groups better. however, the lowest sd estimates were obtained when both were used for normalization, except for average intergroup sd where rps alone performed better. the addition of rpl or actb in the index for normalization did not yield lower sd estimates. the use of just for normalization is in agreement with the results of genorm and also with the rankings produced by normfinder and comparative delta-ct method. the use of as suggested by normfinder does not perform better.table standard deviation for target genes when normalized by different gene combinations ubq rps4 ubq + rsp4 actb + rps4 ubq + rps4 + rpl17 ubq + rps4 + actb ubq + rps4 + rpl17 + actb intragroup and intergroup normalized cq standard deviations averaging the results for the six target genes are shown for fish age + rearing temperature groups when normalized by different candidate reference gene combinations. standard deviations for males and females when normalized by the same combinations are also shown. to assess the robustness of each method, we repeated the stability calculations in fifty subsets of the samples, subsets include samples of each experimental group and another include samples of each experimental group. we evaluated the subsets together since the results show similar trends both with three and two samples per group. since in many studies three genes are used for normalization, we compared not only the whole ranking but also the top genes. the most robust method is clearly bestkeeper sd, which renders an identical ranking as that obtained with the whole data set for a 44% of the subsets and, in 40% of the remaining subsets, it ranks the top three genes correctly. normfinder selected the same top genes also in 88% of the subsets, however the rank order was altered most of the times. on the contrary, the pairwise approaches showed a higher degree of variation, the top genes were different from those in the full dataset in 60% of the subsets for comparative delta ct method and in 66% for genorm.table robustness of the gene stability determination method similarity of subsets stability rankings by each method and the ranking obtained with the whole dataset. subsets are formed by samples per group and another subsets have samples per group. efficiency determination analysis we obtained mean gene efficiencies by lingregpcr, lreanalyzer, dart and pcr-miner for each primer pair and correlations between mean efficiencies by each method for each gene. there is around a 10% difference between the efficiencies calculated by linear fit methods and non linear fit models, meaning that exponential methods might be underestimating efficiency or non linear methods overestimating it. two lreanalyzer efficiency estimates are over 100%, which is theoretically impossible for a pcr reaction, so lreanalyzer is likely overestimating qpcr efficiency. however, despite this 10% efficiency difference, mean efficiencies calculated by the four methods are correlated for each gene, indicating that although they are using different algorithms they are rendering similar relative results. best correlation coefficient and p value are observed between linregpcr and pcr-miner, which might be highlighting the importance of baseline correction since both methods use iterative approaches instead of relying on a fluorescence correction based on the average fluorescence of the first qpcr cycles. also, linregpcr and pcr-miner include several functions to remove outliers, so filtering the reactions before efficiency calculation might also be important to obtain more precise efficiency estimations.table efficiency values for each gene with each efficiency determination method mean efficiency values for each reference gene with lingregpcr, lreanalyzer, dart and pcr-miner. correlation between efficiency determination methods pearson correlation coefficients and p values for mean gene efficiencies with each of the four efficiency determination methods are shown. normalization and efficiency correction on target genes six target genes involved in gonad differentiation were efficiency corrected and normalized by four different combinations of efficiency determination methods and reference gene combinations. for each combination, first, efficiency correction was performed on every cq value of both reference and target genes. afterwards, each target efficiency-corrected cq value was normalized by the reference gene/s efficiency-corrected cq values, obtaining efficiency-corrected delta cq values. we computed mean efficiency-corrected delta cq values and standard deviations for the three temperature groups and also for males and females. two different patterns are shown in the tables, one caused by normalization and the other by efficiency correction. a higher standard deviation is obtained in most of the b2m normalized dataset compared to the ubq + rps normalized ones, which is expected when a gene is not stable. however, this is not true for the amh normal temperature group neither for the gene vasa, suggesting some type of co-regulation. the other trend is observed when comparing the linregpcr efficiency corrected datasets with the pcr-miner corrected ones. pcr-miner produces higher mean cqs increasing the difference between groups.table efficiency-corrected delta cqs by temperature group with each efficiency + reference gene combination mean efficiency-corrected delta cqs and sd values for the three rearing temperatures: high, normal and low; in the four datasets produced after efficiency correction with linregpcr or pcr-miner and later normalization with ubq + rps or b2m. efficiency-corrected delta cqs by sex group with each efficiency + reference gene combination mean efficiency-corrected delta cqs and standard deviation values for males and females in the four datasets produced after efficiency correction with linregpcr or pcr-miner and later normalization with ubq + rps or b2m. furthermore, the use of a gene which presents systematic differences between groups for normalization can lead to changes in the mean cq values of some genes. for example, amh gene expression in each temperature group is severely affected by normalization with b2m, varying from to - at high temperature, from to at normal temperature and from - to - at low temperature. discussion reference gene analysis the four methods commonly used to check the stability of reference genes, comparative delta-ct method, normfinder, bestkeeper and genorm, represent viable strategies, although none of them is currently considered the best one and some problems can arise in certain experimental scenarios. the bestkeeper method is apparently the “common sense” solution to measure stability since standard deviation is a direct measure of variation. however, a gene might show a low standard deviation but still not be a good reference gene if its variation does not reflect the errors produced by sampling, rna extraction and retrotranscription steps. this problem could be circumvented by analyzing the correlations between genes, assuming that the reference genes are not co-regulated. this means that sampling point differences affecting one of the genes should not affect the others, and so the correlations between them would reflect the inter-sample variation produced by the sample processing steps and not by co-regulation due to the experimental conditions. nevertheless, it is risky to assume that genes are not co-regulated because this cannot be easily demonstrated. the genorm and the comparative delta-ct method approaches present the same problem but in addition these methods rank genes mainly by their correlations, to say, genorm establishes the most stable genes by assuming “that the control reference genes are not co-regulated”, and the same happens to the comparative delta-ct method which follows a very similar approach. as a consequence, two co-regulated genes could fully spoil the analysis leading to wrong reference genes. finally, normfinder is not affected by the co-regulated gene drawback since it takes into account intergroup variation, which should be as lower as possible for a good reference gene; however, similarly to bestkeeper, a low overall intergroup and intragroup variation does not necessarily mean that it is a good reference gene. the advantages and disadvantages of each strategy should be taken into account when analyzing putative reference genes according to the experimental scenario. normfinder and genorm are the most extended methodologies to find the optimum reference genes. in many cases, normfinder and genorm algorithms render very similar results, however, discrepancies between the output of normfinder and genorm have been previously described. in these works, the cv has been used to decide which genes should be used for normalization, confirming normfinder results in every case. while normfinder results are non-biased, genorm stepwise exclusion can lead to awkward results by selecting reference genes which in fact are not the most stable. normfinder, bestkeeper and comparative delta-ct stability method results have also been reported to be more consistent among them than with those of genorm, although in other study bestkeeper was reported as the least consistent method. our results with the whole dataset support the high consistency between normfinder and comparative delta-ct method, while bestkeeper results seem to be the least consistent and only correlation values between reference genes seem to suggest a similar ranking. however, when working with different subsets which include a lower number of samples, the pairwise approaches results vary significantly between subsets. this lack of robustness has been described previously: it was shown that the exclusion of a single sample could change the status of one gene from unstable to 2nd most stable gene by genorm. genorm lack of robustness can most likely be explained by the removal of the least correlated gene by pairwise comparison with all the others until only two genes are left, which can lead to stable genes being removed of the analysis early on. robustness is a critical parameter. since experiments are budget limited, it is important to be able to determine correct reference genes with a low number of qpcr reactions. bestkeeper and normfinder appear to be more robust than comparative delta-ct method and genorm in our study. there is not a method to check how much normalization has improved our gene expression data. in principle, a reduction in the cq variability of the gene of interest should be expected, however the highest reduction of this variability would also occur if the gene of interest and the reference gene are co-regulated, so this is a risky strategy. an example of this is observed in vasa standard deviation after normalization, obtaining a lower sd when normalized by a clearly not stable gene than when normalized by ubq + rps suggesting co-regulation between vasa and b2m. the same applies to detecting significant/non-significant results depending on the reference gene used for normalization. this is only useful to stress the importance of choosing a good reference gene, not to choose between one or another since a co-regulated reference gene would lead to non-significant results even if there are expression differences between groups. the fact that there is not a post-control which enables us to check if we have chosen the correct reference gene makes the choice even more critical. every experiment and dataset is different, so the analysis has to be done carefully. given the huge importance of normalization and its great impact on the conclusions, it would be recommended to analyze each case separately, paying attention to details; using any method as a black box can lead up to low confident results. several studies have solved the disagreement between the four methods by ranking them according to the geometric mean of the four ranking numbers for each gene, the lower the mean a gene gets the most stable it is. however, attributing the same weight to every method is arguable, especially because some of these methods include redundant information. this is a practical option without any biological meaning. if the four methods disagree, we recommend instead relying on the ranking provided by normfinder, while ignoring its suggested combination, supported by descriptive statistics like mean, standard deviation and correlations, information offered by bestkeeper or any common statistical package. this approach would enable to assess the two most important and complementary issues: absence of inter-group variation and correlation between reference genes. this approach does not make any previous assumption and has proven to be robust when only a few samples are assayed. to our knowledge this is the first experiment to analyze the stability of reference genes during the gonad development in fish. even in mature organs, there is only one study carried out in zebrafish were testis and ovaries were analyzed separately. however, studies have been carried out in other organs. the stability of several genes was studied in the liver, spleen, kidney, heart, brain, gill and muscle of turbot subjected to edwardsiella tarda infection. gene stability was checked before infection in all the organs together by normfinder and genorm. in that study, out of eight genes, normfinder ranked rpsd as the most stable one, followed by actb, rpl b2m and gapdh among those genes shared with our study, although primer pairs for gapdh were different. we tried to develop primers for rpsd but they were discarded due to late amplification cycle in gonad. ubq and rps were not assayed in that work. gapdh, which has been classically used as a reference gene but recently classified as unstable in several studies, performed badly in both our study and. however, gapdh is ranked as the most stable gene in heart and liver in, which emphasizes the importance of checking reference gene stability in each study separately, reference genes cannot be “exported”. there are two gapdh isoforms in diploid teleost fish as a result of the fish-specific genome duplication event, however the same variant has been analyzed in both studies. there are two similar qpcr studies carried out in flatfish. the first one studied six reference genes during hippoglossus hippoglossus development in sixteen different tissues using bestkeeper, normfinder and genorm. gonads were not included. they assayed actb which was found as one of the least stable genes. the most stable genes found were ef1a and rpl the second study was carried out during larval development in solea senegalensis and hippoglossus hippoglossus. the stability of twelve genes, including ubq, rps actb and gapdh- was checked by genorm and normfinder. the combined stability index of the two species ranked ubq, rps and eef1a as the three best normalization genes by both methods, while actb appears in 5th and 6th place, respectively. interestingly, in solea senegalensis, gapdh- appears ranked 3rd and 1st by genorm and normfinder, respectively. efficiency determination analysis efficiency determination is an essential step in qpcr. constant amplification efficiency in all compared samples is a very important criterion for reliable comparison between samples. it is also crucial for an accurate quantification of gene expression. ideally, the efficiency of an assay should be 100%, which means that during the logarithmic phase of the reaction the pcr product is doubling each cycle. each of the four tested efficiency determination methods differ in their baseline fluorescence determination, type of fit to the log-linear phase of the qpcr reaction, and preprocessing steps to remove outliers. dart is based on a linear regression of the exponential phase of the qpcr reaction. baseline subtraction is determined by fitting a saturation function to the first 2- cycles of the qpcr reaction. then, a linear regression is performed in a 10-fold range around the middle point of the exponential phase, which is calculated using the maximum fluorescence and standard deviation of the fluorescence in the first cycles. linregpcr is also based on a linear regression fit to the log-linear phase of the amplification curve. linregpcr determines baseline fluorescence through an iterative algorithm to get the best fit of the linear regression to 4- points in the log-linear phase of the reaction. then, after baseline subtraction, these points are used for efficiency determination. pcr-miner uses a non-linear regression fit. as linregpcr, baseline fluorescence is determined by an iterative fit to a four-parameter logistic model which also determines the exponential phase of the reaction. then, a three-parameter exponential model is fitted to the exponential phase to determine efficiency. finally, lreanalyzer uses a sigmoidal fit approach. baseline subtraction is determined by averaging 6- cycles fluorescence values. then, efficiency estimates are calculated for each cycle of the qpcr reaction. an lre window is selected with those cycle efficiencies which fit to a linear regression. finally, a derivative of the boltzmann sigmoidal function is used for the calculation of the maximum efficiency of the reaction. logistic models are pure empirical models not designed to be kinetically realistic and rely purely in their good fit to the real-time pcr curve. as previously reported by other study, qpcr curves are not symmetric since they do not have the same curvature at both sides of the inflection point, implying the existence of two or more different mechanisms affecting the efficiency of the reaction and so, making this good-fit models hardly reliable. furthermore, in some reactions sybr green depletion might be the main mechanism leading to the plateau phase of the curve. while sybr green has an impact in the visualization of the real-time reaction, it does not have a connection with the kinetics of the pcr reaction. a recent qpcr study has tried a new approach to assess qpcr efficiency, defining the global efficiency as the sum of denaturing efficiency, annealing efficiency, polymerase binding efficiency and elongation efficiency. the polymerase binding efficiency and the elongation efficiency can be constant provided that there is an excess of polymerase and a long elongation time. however, the denaturing efficiency is constantly decreasing each cycle at the same rate due to thermal damage in both the dna and the polymerase. the annealing efficiency is also decreasing and depends on the proportion of ssdna bound to the primers during this step and total ssdna present. some ssdna chains might bind to its complementary strand instead of to the primers. this efficiency varies from cycle to cycle. this theoretical study was validated in. so, at first, a constant efficiency should not be assumed. still, the qpcr curve shows a large exponential component, since, as confirmed by a previous study, in most cases the best fit to the log-linear region of the qpcr reaction is exponential, suggesting that before and at the log-linear region the qpcr efficiency reduction is low. a recent study analyzed all publicly available efficiency determination methods in a large dataset, included four-point 10-fold dilution series, which allows calculation of the bias of each method. similarly to our results, they report linregpcr and dart to produce an underestimation of efficiency and pcr-miner and lre analyzer an overestimation. they also analyze different parameters and find linregpcr and pcr-miner amongst the best methods for most of the evaluated characteristics, for example precision and resolution, performing better than lreanalyzer and dart. the reader is encouraged to consult to learn more about the different efficiency determination methods, their characteristics and performance differences. this study is the most complete on qpcr efficiency determination methods done so far. both linregpcr and pcr-miner performed similarly and produced highly correlated efficiency estimates in our study. the main difference was that while linregpcr underestimates efficiency, pcr-miner overestimated it. knowing this, linregpcr is probably the best choice for the average qpcr researcher since it will not produce erroneous significant differences between groups or an overestimation of the fold change. however, linregpcr might not be the best option for clinical purposes, where the method of choice should be considered depending on the consequences of a false positive/overestimation or a false negative/underestimation. there are a good number of alternative efficiency estimation algorithms, however they are implemented as extensions to the open source statistical programming environment r and probably not available for most researchers, so they have not been analyzed here. still, linregpcr and pcr-miner perform as well or better than all other methods, as shown in. as a final remark, although clearly improved our understanding of the different efficiency determination methods available, currently there is not a clear best method for estimating qpcr efficiencies. however, the publication of a theoretical study based on the pcr kinetics which defines the overall pcr efficiency as the product of the efficiency of each of the separate steps and its experimental validation is a good step towards finding a biologically meaningful solution. similar approaches will be likely applied in more studies in the near future provided they are implemented in appropriate user friendly softwares for the whole research community. normalization and efficiency correction on target genes the effects of normalization with wrong genes are important; a high standard deviation will produce higher p values and, so, possibly lead to missing biologically relevant differences. even worse, the use of a regulated gene for normalization, will lead to changes in gene values which can end in misguided results. the effect of efficiency correction, though not so dramatic, is also important since it can lead to overestimation of differences between groups.  we found the ranking produced by normfinder method as the most reliable one to choose reference genes for qpcr analysis when results differ between gene stability determination methods. normfinder information should be complemented by the descriptive statistics offered by bestkeeper, especially the correlation coefficient. accordingly, we found that ubq and rps should be used as reference genes to study turbot gonad development from up to 135 days post fertilization. we found pair-wise methods to be less robust than normfinder and bestkeeper and also the suggested normfinder two genes combination not reliable. we also recommend the use of linregpcr for efficiency determination for research purposes, however, efficiency determination is still a matter of discussion and probably new improved models will be published in the upcoming years. 
we also analyzed the efficiency calculation softwares linregpcr, lreanalyzer, dart and pcr-miner and we recommend linregpcr for research purposes since it does not systematically overestimate efficiency. we analyzed gene stability of six reference genes: rps rpl gapdh, actb, ubq and b2m using the comparative delta-ct method, bestkeeper, normfinder and genorm approaches in gonad samples of turbot. supported by descriptive statistics, we found normfinder to be the best method, while on the other side, genorm results proved to be unreliable. our dataset is composed by gonad samples of turbot at different development times reared at different temperatures. there is no agreement on which method is the best to detect the most stable genes neither on how to perform efficiency determination. according to our analysis, ubq and rps were the most stable genes, while b2m was the least stable gene. since females largely outgrow males, identification of genes related to sex determination, gonad development and reproductive behavior, and analysis of their expression profiles are of primary importance for turbot industry. gene expression analysis by reverse transcription quantitative pcr is the most widely used method for analyzing the expression of a moderate number of genes and also for the validation of microarray results. our results indicate that normfinder and linregpcr are the best approaches for reference gene selection and efficiency determination, respectively. several issues are crucial for a successful qpcr study, particularly the selection of internal reference genes for normalization and efficiency determination. qpcrreference genesamplification efficiencyturbotscophthalmus maximusgonadissue-copyright-statement© the author 2014 we also recommend the use of ubq and rps for normalization of gonad development samples in turbot. in this study we offer a comprehensive evaluation of the characteristics of reference gene selection methods and how to decide which one is more reliable when they show discordant outcomes.
 gene expression analysis by reverse transcription quantitative pcr is the most widely used method for analyzing the expression of a moderate number of genes and also for the validation of microarray results. several issues are crucial for a successful qpcr study, particularly the selection of internal reference genes for normalization and efficiency determination. there is no agreement on which method is the best to detect the most stable genes neither on how to perform efficiency determination. in this study we offer a comprehensive evaluation of the characteristics of reference gene selection methods and how to decide which one is more reliable when they show discordant outcomes. also, we analyze the current efficiency calculation controversy. our dataset is composed by gonad samples of turbot at different development times reared at different temperatures. turbot is a relevant marine aquaculture european species with increasing production in the incoming years. since females largely outgrow males, identification of genes related to sex determination, gonad development and reproductive behavior, and analysis of their expression profiles are of primary importance for turbot industry.  we analyzed gene stability of six reference genes: rps rpl gapdh, actb, ubq and b2m using the comparative delta-ct method, bestkeeper, normfinder and genorm approaches in gonad samples of turbot. supported by descriptive statistics, we found normfinder to be the best method, while on the other side, genorm results proved to be unreliable. according to our analysis, ubq and rps were the most stable genes, while b2m was the least stable gene. we also analyzed the efficiency calculation softwares linregpcr, lreanalyzer, dart and pcr-miner and we recommend linregpcr for research purposes since it does not systematically overestimate efficiency.  our results indicate that normfinder and linregpcr are the best approaches for reference gene selection and efficiency determination, respectively. we also recommend the use of ubq and rps for normalization of gonad development samples in turbot. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords qpcrreference genesamplification efficiencyturbotscophthalmus maximusgonadissue-copyright-statement© the author 2014 