Evaluation of qPCR reference genes in GH-overexpressing transgenic zebrafish (Danio rerio)

Reference genes (RGs) must have a stable expression in tissues in all experimental conditions to normalize real-time quantitative reverse transcription PCR (qRT-PCR) data. F0104 is a highly studied lineage of zebrafish developed to overexpress the growth hormone (GH). It is assumed that the transgenic process may influence the expression levels of commonly used RGs. The objective of the present study was to make a comprehensive analysis of stability of canditade RGs actb1, actb2, b2m, eif2s2, eef1a1, gapdh, rplp2, rpl7, rpl13α, tuba1, and rps18, in gh-transgenic and non-transgenic zebrafish. Liver, brain, intestine and muscle samples from both groups had qRT-PCR results analyzed by dCt, geNorm, NormFinder, BestKeeper, and RefFinder softwares. Consensus analyses among software concluded that rpl13α, rpl7, and eef1a1 are the most stable genes for zebrafish, considering the studied groups and tissues. Gapdh, rps18, and tuba1 suffered variations in stability among different tissues of both groups, and so, they were listed as the genes with lowest stability. Results from an average pairwise variations test indicated that the use of two RGs would generate reliable results for gene expression analysis in the studied tissues. We conclude that genes that are commonly used in mammals for qRT-PCR assays have low stability in both non-transgenic and gh-transgenic zebrafish reinforcing the importance of using species-specific RGs.

tissue collection. The animals were euthanized using overdose of tricaine methanesulfonate (MS-222, Sigma-Aldrich, USA) at 400 mg L −1 . The animals were weighed and measured. Samples of liver, brain, intestine, and caudal muscle of adult zebrafish were collected from both groups. The samples collected from 3 animals were mixed to compose a tissue pool, totalizing seven replicate tissue pools (from liver, brain, intestine, and muscle) for each experimental group. After the collection, samples were stored in liquid nitrogen (N 2 ) until RNA extraction.
RNA extraction and cDNA synthesis. RNA extraction and cDNA synthesis were performed according to Silveira and collaborators 4 . Briefly, total RNA was extracted from samples using TRIZOL Reagent (Thermo Fisher Scientific, USA), following the manufacturer's instructions. RNA was treated with the DNAse-free kit (Ambion, USA) to remove genomic DNA contamination. Subsequently, RNA concentration and purity were measured using a NanoVue Plus spectrophotometer (GE Healthcare Life Science, USA) and after, the samples were stored at − 80 °C. Complementary DNA (cDNA) was obtained using High Capacity cDNA Reverse Transcription (Applied Biosystems, USA) according to the manufacturer's recommendation at a concentration of 500 ng. Finally, cDNA was stored at − 20 °C until use.
Gene expression analysis by qRT-PCR. The primers used in this study are presented in Table 1. The qRT-PCR was performed using the SYBR Green PCR Master Mix (Applied Biosystems, USA), and carried out in the Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, EUA). The amplification conditions were 95 °C for 10 min, 40 cycles at 95 °C for 15 s, and 60 °C for 1 min, followed by the conditions necessary for the calculation of the melting curve. The study used primers that have been validated for zebrafish (Table 1) for 11 candidate RGs (actb1, actb2, b2m, eef1a1, eif2s2, gapdh, rpl13α, rpl7, rplp2, rps18, and tuba1) belonging to different physiological pathways and representing different functions to evaluate expression stability. All reactions were performed in duplicate. Data analysis. The normality distribution of biometric data was evaluated by the Shapiro-Wilk test. The unpaired Student's t test (one-tailed) was used to compare the mean values between transgenic and non-transgenic groups. The obtained values were expressed as g ± standard error of the mean (SEM) and mm ± SEM for length and mass, respectively. Analyses of gene expression stability were performed using the comparative method delta-Ct (dCt) 23 and geNorm 24 , NormFinder 25 , BestKeeper 26 , and RefFinder statistical approaches as previously reported 27,28 . The results were expressed as mean ± standard deviation (SD) for the dCt method; Scientific RepoRtS | (2020) 10:12692 | https://doi.org/10.1038/s41598-020-69423-y www.nature.com/scientificreports/ average expression stability value for GeNorm and NormFinder; Pearson correlation coefficient (r) for the Best-Keeper; and geometric mean of rating values for RefFinder. Average pairwise variations (V) were calculated by geNorm between the normalization factors NFn and NFn + 1 to indicate whether inclusion of an extra RG would add to the stability of the normalization factor. The seven best RG candidates out of the eleven evaluated were used to calculate the V values. The use of seven genes was, independently of the results, due the inclusion of unstable genes (represented by the four worst genes that were excluded) which would be worse than the use of a few suitable RG.

Verification.
A protocol was established to demonstrate the differences found in gene expression stability between RGs by algorithms. The number of RGs used as 'the best' was determined by the pairwise variation analysis (the result was that two genes were necessary to normalize reactions efficiently, see in results section). The worst and the third best RG candidates of each tissue were used as target genes in gene expression analysis. Thus, the comparison of the implications of the use of unsuitable and suitable RGs in experimental analysis was possible. This procedure was performed in two parallel ways. In one way the qRT-PCR reactions were normalized with the two bests RGs of each tissue. In other way the reactions were normalized by the two bests RGs, considering all tissues grouped. This made it possible to compare if there are differences between the use of a single set of RGs for all tissues and the use of various sets of RGs, one set for each tissue. Because the third best liver RG was already committed with the normalization of 'the first way' , this gene was replaced by the fourth best RG. Furthermore, this fact precluded 'the second way' of analysis because a redundancy would originate, since the same gene would be used as target gene and RG. The qRT-PCR data from the target genes were analyzed using the 2 −ΔΔCt method, according to Livak and Schmittgen (2001) 5 . The mean between the Cts of the RGs that were used was applied to the Livak and Schmittgen method. The normality distribution of the data was evaluated by the Shapiro-Wilk test. The unpaired Mann-Whitney test (two-tailed) was used to compare the mean values. The results were expressed as mean ± SD.
Ethics approval. All the protocols used in this study were approved by the Ethics Committee of the Federal University of Rio Grande (FURG), Brazil, under the code 23116.008403/2018-32.

Results
Transgenic phenotype expression. The transgenic group expressed the phenotype attributed by the transgenesis, being statistically longer and heavier (p < 0.05) (Fig. 1a,b, respectively) in comparison to the length and body mass of the non-transgenic group. The non-transgenic group presented 31.35 ± 0.68 mm and the  F: TGC ATG GCC GTT CTT AGT TG  R: AGT CTC GTT CGT TAT CGG  AATGA   FJ915075   4   Glyceraldehyde-3-phosphate dehydrogenase  gapdh  Glycolytic enzyme   F: GAT GGT CAT GCA ATC ACA  GTCTA  R: ATC ATA CTT GGC AGG TTT  CTCAA   BC095386   15   Beta-2-microglobulin  b2m  Beta chain of major histocompatibility  F: GCC TCC ACC CCA GAG AAA GG  R: GCG GTT TTT ATT TAC ATG TTG

RefFinder analysis.
The RefFinder is a tool used to construct a consensus comprehensive ranking of stability of RGs among comparative dCt method and the geNorm, NormFinder, and BestKeeper algorithms using the calculation of geometric mean for the ranks calculated by each of these other methods. Candidate genes with the lowest geometric mean are most stable. A consensus graph resulted from RefFinder analysis considering all the analyzed methods, tissues, and combining both non-transgenic and gh-transgenic groups was generated and is presented together with the results of the other methods that supported the consensus analysis ( Fig. 2a- (Fig. 2a). Furthermore, consensus results of RefFinder, individually discriminated by non-transgenic and gh-transgenic groups but combining all tissues, are also shown (Fig. 3a,b). In the non-transgenic group, the comprehensive consensus classification constructed by RefFinder (Fig. 3a) (10.74). In the gh-transgenic group the classification order ( Fig. 3b) was: rpl13a (1.57) > rpl7 (2.59) > eef1a1 (3.13) > rplp2 (3.13) > actb2 (3.34) > actb1 (4.56) > eif2s2 (7.24) > b2m (7.74) > tuba1 (9.00) > rps18 (10.00) > gapdh (11.00).
Finally, in addition to graphic rankings, RefFinder consensus results are presented for each group and by stability in each tissue (liver, brain, intestine, and muscle) (Tables 2, 3, 4). According to the consensus analysis of the tissues of both groups together, the most stable gene was considered the rpl13α for brain and intestine, occupying the second and third place for muscle and liver, respectively (Table 2). Also, in consensus analysis, the gapdh gene was the most unstable and unsuitable gene, except for the intestine, when gapdh was the second most unstable (Table 2).

Non-transgenics
Gh-transgenics  (Table 3), the best specific RG for brain and liver was actb1. The most stable genes in intestine and muscle were rpl13a and actb2, respectively. The worst genes found were rps18 in brain and muscle and gapdh in intestine and liver. In the gh-transgenic group (Table 4), the most stable genes were the actb1 in intestine and muscle; eef1a1 in the brain; and rplp2 in the liver. In this group, rps18 was the most inappropriate gene for intestine and liver, while for brain and muscle were eif2s2 and b2m, respectively ( Table 4).
The comparative dCt method analysis. The comparative dCt method was used to select the most stable RG. A low average of SD value represented a low expression variance, or high stability. The dCt general analysis of both groups revealed the follow ranking order:   Fig. 1b).

Comparative dCt method r p l 1 3 a e e f 1 a 1 r p l 7 r p l p 2 a c t b 2 a c t b 1 b 2 m e i f 2 s 2 t u b a 1 g a p d h
GeNorm analysis. The geNorm analysis was carried out for the resulting data set following transformation of the Ct values into relative quantities through the 2 (minimum Ct value in a set sample − Ct value of a sample) formula and compared pairwise variation (SD values) for each gene pair. Then, the geometric mean of SD values was used to calculate the M-value. The generation of low average expression stability represents a low variance. The geNorm general ranking order for both groups was: (Fig. 2c). In the non-transgenic     Fig. 2a). The geNorm analysis revealed the following order of stability classification of genes in the gh-transgenic group, independent of tissue: Fig. 2b).

NormFinder analysis.
NormFinder was applied to analyze the most stable evaluated genes. A low mean expression stability represents a low variance. The general consensus result for both groups generated by  Fig. 4b).
Pairwise variation. The analysis of average pairwise variations showed that the lowest V value was between 2 and 3 RG (V2/V3 = 0.202) indicating that the inclusion of a third RG would not be necessary to normalize gene expression. All the other V values resulting of the analysis were higher than the V value of 0.202, indicating that the inclusion of a third, fourth, fifth, sixth, or seventh RG would be superfluous to obtaining reliable results (Fig. 4).  24 . The manual of the geNorm software suggests that the 0.15 value must not be taken as a strict cut-off and may vary according to each experimental design. Therefore, it is proposed that two of the best candidate RG should be used to obtain accurate and reliable results using the different tissues of the gh-transgenic zebrafish and its non-transgenic counterpart from the F0104 lineage.
Verification. The verification was made using two selected RGs, as this number has been the best established by the average pairwise variation analysis. Considering the bests RGs of all grouped tissues, rpl13a and rpl7 were selected for qRT-PCR data normalization. To analyze the tissues separately, the results found for each tissue were considered: for brain and intestine, rpl13a and eef1a1; for muscle, rpl13a and actb2.
For the brain analysis, gapdh was chosen to be tested as a target gene, as it was found to be the less stable RG in this tissue, while b2m was choosen as a suitable RG, as it was the third best in previous analyses (stressing that the two most stable genes were used to normalize the test) (Fig. 5). It is possible to observe that even considering only the dataset taken the brain as for all tissues grouped, gapdh presented a high variation of gene expression between non-transgenic and gh-transgenic groups, being significantly more expressive (p < 0.05) in the gh-transgenic group (Fig. 5a,b), while for b2m, this expression did not present statistical variation (p > 0.05) between the groups (Fig. 5c,d). . Verification of the suitability of the reaction when normalized with the candidate reference genes ranked as the best for brain tissue (a,b) and for all tissues grouped (c,d), using brain gene expression results. As target genes, it was chosen the worst (a) and the third best (b) for brain. In all the graphics, the Y-axis represents mRNA relative expression, and X-axis indicates the experimental group: non-transgenic and gh-transgenic zebrafish of F0104 lineage. Data are expressed as mean ± standard error of the mean. Asterisk indicate significant difference mean values between groups (n = 21; p < 0.05).
Regarding the intestine, the target genes chosen were rps18 (the worst) and actb1 (third best) (Fig. 7). In this analysis, both genes we found to not vary between groups (p > 0.05), considering the results of only this tissue (Fig. 7a,b) as for all grouped tissues (Fig. 7c,d).
As for the liver, the target genes selected were gapdh (the worst) and eef1a1 (second best) (Fig. 8). The gapdh gene was significantly overexpressed (p < 0.05) in the gh-transgenic group (Fig. 8a), demonstrating its instability, which had already been verified by the stability analysis. On the other hand, the eef1a1 gene presents expressions levels without significant difference (p > 0.05) between the non-transgenic and gh-transgenic groups (Fig. 8b).

Discussion
In this study, we present the most extensive evaluation of candidates for normalizing genes for qRT-PCR in gh-transgenic and non-transgenic zebrafish groups. In both groups, we evaluated 11 genes in 4 different tissues, and the results were analyzed using 5 different bioinformatic approaches. In addition, consensus analyses were performed between the organs and between the groups. Figure 6. Verification of the suitability of the reaction when normalized with the candidate reference genes ranked as the best for muscle tissue (a,b) and for all tissues grouped (c,d), using muscle gene expression results. As target genes, it was chosen the worst (a) and the third best (b) for muscle. In all the graphics, the Y-axis represents mRNA relative expression, and X-axis indicates the experimental group: non-transgenic and ghtransgenic zebrafish of F0104 lineage. Data are expressed as mean ± standard error of the mean. Asterisk indicate significant difference mean values between groups (n = 7; p < 0.05).
Scientific RepoRtS | (2020) 10:12692 | https://doi.org/10.1038/s41598-020-69423-y www.nature.com/scientificreports/ F0104 is a zebrafish transgenic lineage, developed in 2004, that overexpresses the GH from silverside (Odontesthes argentinensis) in a constitutive promoter (actb), and also carries the gene of green fluorescent protein (GFP) as a marker 17 . This lineage has been very important as an animal model, not only for elucidating the effect of GH overexpression in tissues and its effect in different physiologic systems 17-19, 21, 29-33 , but also to evaluate the interaction of this hormone with other proteins in the organism [34][35][36][37] . To achieve all those elucidations, quantitative PCR analysis has been applied in diverse studies using this lineage. Different RGs have been used to normalize and validate the results, but the most widely used was eef1a1 18,19,21,[30][31][32][35][36][37][38] . As eef1a1 was the third most stable RG found here (Fig. 2a) for the same lineage, those previous results corroborate our data. In general, our analyses found the eef1a1 gene to be highly stable, occupying between the first and at least the third position in the rankings-except in the BestKeeper and NormFinder analysis in the gh-transgenic group-demonstrating to be one of the most specific and stable gene for analysis in zebrafish in both groups. Other studies have also examined the stability of RGs in non-transgenic zebrafish, showing that eef1a1 is a stable gene for qRT-PCR normalization 11-13, 15, 16 .
However, eef1a1 was not ranked as the most stable gene among all analyzed in this study. Here, we could verify that in a RefFinder consensus analysis between groups and tissues, the most stable gene was rpl13α, followed by rpl7 (Fig. 2). However, when the consensus analysis was verified considering the groups and the tissues individually, the positions in stability rankings can be altered (Tables 2, 3, 4; Supplementary Figs. 1-4). In this case, rpl13α was not the most stable in all the tissues from the different groups, but only in the intestine in the non-transgenic group ( Table 3). The rpl7 gene, which occupied the second position in the consensus of general Figure 7. Verification of the suitability of the reaction when normalized with the candidate reference genes ranked as the best for intestine tissue (a,b) and for all tissues grouped (c,d), using intestine gene expression results. As target genes, it was chosen the worst (a) and the third best (b) for intestine. In all the graphics, the Y-axis represents mRNA relative expression, and X-axis indicates the experimental group: non-transgenic and gh-transgenic zebrafish of F0104 lineage. Data are expressed as mean ± standard error of the mean. Asterisk indicate significant difference mean values between groups (n = 7; p < 0.05).
Scientific RepoRtS | (2020) 10:12692 | https://doi.org/10.1038/s41598-020-69423-y www.nature.com/scientificreports/ stability, was shown to be the third most stable in the liver of the non-transgenic group (Table 3), and in the muscle of the gh-transgenic group ( Table 4). The high stability of rpl13α in zebrafish liver has also been demonstrated by Xu et al. (2013) 39 when fish were stressed with arsenate. The rpl13α was also shown to be a more stable and reliable gene for chemical treatments in a recent study using zebrafish 1 , consistent with previous studies 11,40,41 . Considering the analysis by tissues, some studies corroborate our findings. Jaramillo et al. (2018) 16 found that the best RGs in the brain were rpl8 and actb1; the former was not investigated in our study, but the latter was also found by us as the best for that tissue in the non-transgenic group (Table 3), and the fourth best in the gh-transgenic group ( Table 4). The results obtained by Lang et al. (2016) 15 were compatible with our findings, where rpl13α, rpl7, actb2 and eef1a1 were the most suitable RGs in liver when exposed in different concentration groups to cadmium and subsequently infected with Aeromonas hydrophila. The same study found that in the intestine the best RGs were ef1a1, rnf7, rplp2 and rpl13α. In our study, we did not investigate rnf7, and the other genes appeared in different ranking positions: rpl13α was the best RG for intestine in non-transgenic zebrafish; ef1a1 was at 5th and rplp2 was at 8th place (Table 3).
Gapdh and rps18 showed to be the most inappropriate genes for both groups in terms of stability-with the exception of the analysis carried out using the algorithm BestKeeper, where surprisingly gapdh occupied the first place in the gh-transgenic group when considering all tissues together ( Supplementary Fig. 4). In general, both gapdh and rps18, when analyzed separately in the tissues, occupied the last positions in almost all analyses and statistical programs. This result demonstrates that they are unsuitable RGs, and when they are used, they can seriously compromise the analysis and the accuracy of qRT-PCR results. Also considering the rankings generated by our analysis, tuba1 was the third less stable gene in both groups separately (Fig. 3) and together in most cases (Fig. 2a-d). For wild zebrafish, gapdh 11,12 and tuba1 12 have already been reported as unstable genes and unsuitable RG.
Another result found in our study was that the use of two RGs to interpret qRT-PCR results was enough (Fig. 4). This result is positive, as needing less genes to evaluate the analysis makes it cheaper and simpler. Furthermore, this result reinforces the methods and results obtained in previous studies with gh-transgenic zebrafish, which also used two RGs for qRT-PCR analyses 30,31,33,[36][37][38] . Some authors of these previous studies also analyzed candidate RGs with computational algorithms (i.e. geNorm) 31,32,36,38 . However, this is the first study to use an integrative approach with various methods to evaluate the variability of candidate RGs.
The verification analysis, which was made for each tissue, considering the best and worst candidate RGs found as target genes by each tissue and for all grouped tissues, has also shown the difference of gene expression between suitable and unsuitable RGs, highlighting the risk of choosing few and unstable genes to normalize the dataset. Besides it was possible to test if each tissue needs to be evaluated individually or if the grouped results are enough to access the bests RGs to use in a specific experiment with the F0104 lineage (Figs. 5, 6, 7). The second situation (grouped results) can make an experiment quicker and less costly because fewer RGs are needed to obtain reliable results. The results from brain, muscle, and liver clearly show the difference between a suitable and unsuitable (gapdh) RG, showing a stable expression of the second/third best RG and variable expression of the gapdh RG (Figs. 5, 6, 8). In the intestine, even though the gene expression did not presenting a statistic difference between the groups, a stronger trend of variability was observed in rps18 expression in both analysis, when compared to the actb1 results (Fig. 7).
Yet, the verification method revealed that the gapdh is overexpressed in muscle and liver (Figs. 6, 8), while it is downregulated in the brain of gh-transgenic zebrafishes (Fig. 5). This gene encodes to an enzyme with diverse Normalization with the bests candidates RGs of all tissues grouped Figure 8. Verification of the suitability of the reaction when normalized with the candidate reference genes ranked as the best for all tissues grouped (a,b), using liver gene expression results. As target genes, it was chosen the worst (a) and the second best (b) for liver. In the two graphics, the Y-axis represents mRNA relative expression, and X-axis indicates the experimental group: non-transgenic and gh-transgenic zebrafish of F0104 lineage. Data are expressed as mean ± standard error of the mean. Asterisk indicate significant difference mean values between groups (n = 7; p < 0.05).
Scientific RepoRtS | (2020) 10:12692 | https://doi.org/10.1038/s41598-020-69423-y www.nature.com/scientificreports/ function for maintaining cell homeostasis, but its main function is to break down glucose in the glycolysis. The GH hormone induces IGF-1 production and also promotes postnatal body growth 42 , which explains the impair of glycolysis in the liver and muscle of fish in the gh-transgenic group. Maybe, the downregulation of gapdh in the brain can be explained by the downregulation of glucose transport-related genes in the brain of fish form the F0104 lineage 37 . Furthermore, the maintenance of the response pattern of target gene expression independently of the set of genes (or the two best RGs of specific tissue or the two best consensual RGs after tissue grouping) selected to normalize the qRT-PCR (see in Figs. 5,6,7,8) indicates that a unique consensual set of genes can be used to normalize the data of various tissues. Hence, instead of selecting two RGs for the brain, two RGs for the intestine, two RGs for the liver, and two RGs for the muscle, the researcher could select only two RGs to normalize the results of qRT-PCR for these four different tissues.
Regarding the differences in stability of the candidate genes, probably this result happened due to the different functions and heterogeneity of the organs 43 . Besides that, the different results presented by each algorithm are expected due to the distinct statistical approach that was used to construct the rankings 44 . The discrepancies in the ranks of the most stable genes obtained by each software reflect the differences of the algorithms and procedures of each method and demonstrate the need to evaluate the RGs using several bioinformatic algorithms 16,23 .
Therefore, it is important to consider that RGs for qRT-PCR can be altered by different experimental conditions, organs, gender, stage of development and chemical concentration 11,16,24 . So, it is impossible to find suitable RGs that exhibit constant expression pattern for all species and under all experimental conditions 6,45 . According to the Minimum Information for Publication of Quantitative Real Time PCR Experiments (MIQE) 46 , a RG needs to be validated and established for each species and for different physiological conditions. Our results show that although the literature has already demonstrated alterations in the metabolism of the F0104 lineage 19, 21, 29-35, 38, 47, 48 , these animals, in comparison to non-transgenic animals, did not suffer notable variations in the expression of most of the analyzed genes.
In summary, for the first time, RGs for non-transgenic and gh-transgenic zebrafish of F0104 lineage were evaluated for gene expression stability and these results support the use of two stable RGs for proper interpretation of the qRT-PCR results among different tissues collected in specific experiments (brain, instestine, liver, and muscle in this case). The actb1 and actb2 genes were shown to be more stable and suitable for intestine and muscle, while in brain and liver, the eef1a1 gene was shown to be the most suitable. In addition to the rankings of the most suitable candidate RGs in each tissue, we presented a consensual analysis that showed that the genes rpl13α, rpl7 and eef1a1 were the most stable between tissues and groups. More than 30 studies have already been performed using zebrafish of the F0104 lineage and our conclusions can be used in countless future studies which analyze the results of qRT-PCR in this lineage, the gold standard technique to evaluate gene expression in experimental samples. In addition, the present study showed the capital importance of using various methods to evaluate stability. As seen in our case, if used individually, BestKeeper could have led to a complete misinterpretation, pointing gapdh as the most stable gene in the gh-transgenic group considering all tissues together. Moreover, the risks of a mistaken choice of RG and its consequences for the interpretation of results were demonstrated visually through the verification steps. Finally, this study allowed us to conclude that genes that are commonlly used in mammals for qRT-PCR assays have low stability in both non-transgenic and gh-transgenic zebrafish. This reinforces the importance of using species-specific genes, instead of 'consecrated' RGs for normalization and, wherever possible, using RGs that have been previously checked for stability based on the animal model used and the specific experimental conditions. Therefore, in a near future more screenings studies about RGs in other genetically engineered lineages will be needed to guarantee reliable interpretations about the expressions of the studied genes.