Identification of optimal endogenous reference RNAs for RT-qPCR normalization in hindgut of rat models with anorectal malformations

Background Quantitative real-time polymerase chain reaction (RT-qPCR) is a sensitive method for quantifying mRNA abundance. With relative expression analysis, however, reliable data output is dependent on stably expressed reference genes across the samples being studied. In anorectal malformations (ARMs), there is limited data on the selection of appropriate reference genes. Purpose This study was aimed to investigate the optimal reference genes for PCR in ARM rat models. Methods We selected 15 commonly used reference genes (Rps18, Actb, B2m, Gapdh, Ppia, Hprt1, Pgk1, Ywhaz, Tbp, Ubc, Rps16, Rpl13a, Rplp1, Sdha, and Hmbs) as candidate reference genes and detected their mRNA expression in ARM samples by RT-qPCR. The expression stability and variability of these transcripts were subsequently evaluated using four methods (geNorm, NormFinder, comparative ΔCt, and BestKeeper). Results The abundance of the candidate reference genes was qualified by RT-qPCR and the cycle threshold (Ct) values ranged between 14.07 (Rplp1) and 21.89 (Sdha). In the overall candidate genes, different variations existed across the different algorithms. A comprehensive analysis revealed that Rpl13a ranked first among the relatively stable genes, followed by Ywhaz, Rps18, Sdha, and Hmbs. Conclusions The most stable reference genes for RT-qPCR were Rpl13a, Ywhaz, and Rps18 in ETU-induced ARMs in rat fetus. This study provided a foundation for reference gene selection for future gene expression analyses.


INTRODUCTION
Anorectal malformations (ARMs) are common congenital gastrointestinal malformations manifested by anal stenosis, ectopic anus, and urethrorectal fistula (Bai et al., 2004;De Blaauw et al., 2013;Wijers et al., 2014). The specific mechanism of ARMs still remains unclear although numerous studies have revealed that various genes are involved in the development of this disease (Draaken et al., 2012;Wijers et al., 2014). As clinical samples are difficult to obtain, animal models are commonly used to study ARMs. The most commonly used model is ETU-induced ARMs in rat fetus (Macedo, Martins & Meyer, 2007;Qi, Beasley & Frizelle, 2002).
In ARM-related studies, it is inevitable to analyze gene expression. Quantitative real-time polymerase chain reaction (RT-qPCR) is commonly used because of its high efficiency and sensitivity. While analyzing target gene expression by RT-qPCR, there are many strategies by which relative quantification is often adopted (Yuan et al., 2006). In this course, normalization to extensively and stably expressed endogenous reference genes is commonly used to reduce the variations caused by input RNA amount or reverse transcription efficiency (Bustin et al., 2009;Giulietti et al., 2001).
It is noteworthy that this approach is based on the presupposition that the transcript abundance of selected reference genes is stable in all samples under various conditions (Giulietti et al., 2001). However, researchers have found that what was previously thought as "stably expressed" reference genes are not actually strictly invariant. In fact, the reference genes can be affected by various factors (Harris, Reeves & Phillips, 2009;Mughal et al., 2018). For instance, diseases or experimental treatments can alter the expression stability of the reference genes (Aggarwal et al., 2018;Giulietti et al., 2001;Neerukonda et al., 2016). Physical or chemical factors can also influence the variability of the reference genes (Giulietti et al., 2001;Mahanty et al., 2017;Svingen et al., 2015;Yuanyuan et al., 2015). Moreover, the development stage or even seasonal factor has an effect on the optimal reference gene selection (Liu et al., 2014;Macabelli et al., 2014;Mughal et al., 2018). The reference gene expression instability may be caused by organ composition and cell types (Guenin et al., 2009).
Based on the above conditions, some researchers have suggested that normalization against a reference gene should not be acceptable unless it has been clearly confirmed that the reference gene is invariantly expressed in the experimental and control conditions (Aggarwal et al., 2018;Bustin et al., 2009). In order to obtain reliable results in RT-qPCR, systematic validation of reference genes is essential (Guenin et al., 2009).
To the best of our knowledge, there is no study identifying the optimal reference genes in ETU-induced ARM rats. Therefore, systematic evaluation is essential for obtaining accurate RT-qPCR results. In this study, we analyzed the expression of 15 candidate reference genes in the hindgut of normal and ETU-induced ARMs in rat fetuses during critical time points for anorectal development (Bai et al., 2004;Faria, De Jesus Simoes & Martins, 2015;Macedo, Martins & Meyer, 2007). In addition, we systematically evaluated and ranked the expression stability of these genes aiming to identify an optimal reference gene for RT-qPCR in ARM rat models.

Ethical statement
This study was approved by the China Medical University Ethics Committee (no. 2015 PS213K). All procedures were performed according to the guidelines for the care and use of laboratory animals.

Animal models and sample preparation
A total of 18 mature female Wistar rats (age: 7-9 weeks old, body weight: 250-300 g) were provided by the Experimental Animal Center of the Shengjing Hospital of China Medical University (Shenyang, Liaoning, China). The animals were raised in a specific pathogen-free environment of The Key Laboratory of Health Ministry for Congenital Malformations (Shenyang, Liaoning, China). Surgeries were performed after euthanasia by intraperitoneally injecting sodium pentobarbital, and all efforts were made to minimize the pain and suffering of the animals.
Anorectal malformations models were developed according to an earlier report (Bai et al., 2004;Long et al., 2018). Briefly, nine pregnant rats were administered 125 mg/kg of 1% ETU (Sigma-Aldrich, Merck Millipore, Darmstadt, Germany) by single gavage on GD10 (gestational day (GD); GD0, presence of sperm in vaginal smear after overnight mating), while the other nine rats received normal saline as control. The presence of ARMs was determined by microscopy. The hindguts of the embryos were removed from the surrounding tissues (Fig. S1). The tissues were then immediately frozen in liquid nitrogen for RT-qPCR.
A total of 123 ETU-induced ARM fetuses (A) and 115 saline-treated normal fetuses (N) were obtained by cesarean section on GD16. The incidence of ARMs in the ETU-treated fetuses was 85.8% (97/113). All hindgut samples from each age group were immediately frozen in liquid nitrogen for subsequent RNA extraction.

RNA extraction and cDNA synthesis
The total RNA was isolated using the miRNeasy Mini kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer's instructions. RNA was qualified and quantified using Tecan Infinite 200 multimode plate reader (Tecan AG, Hombrechtikon, Switzerland) and formaldehyde agarose gel electrophoresis. Only samples with RNA absorption ratios of A260/280 ranging from 1.9 to 2.1 and RNA integrity number 7 from the Agilent 2200 RNA assay (Agilent Technologies, Inc., Santa Clara, CA, USA) were used for cDNA synthesis. The cDNA was synthesized from 500 ng of total RNA in a final volume of 10 ml using PrimeScript RT reagent kit (Takara Biotechnology Co. Ltd., Dalian, China) according to the manufacturer's instructions.

Amplification by RT-qPCR
The PCR efficiency of each primer pair was tested using calibration curves, and the efficiency rate was calculated according to the following equation: PCR amplification efficiency ¼ 10 -1/slope -1 (Bustin et al., 2009). The gene expression levels were assessed using SYBR Premix Ex Taq II kit (Takara Biotechnology Co. Ltd., Dalian, China). Each reaction system contained 2 ml cDNA template, 0.8 ml of each forward and reverse primers (10 nM), 6 ml RNase-free water, 0.4 ml ROX II, and 10 ml SYBR Green premix Ex Taq II (Tli RNaseH Plus) in a final volume of 20 ml. The amplification cycles were set as follows: denaturation at 95 C for 30 s followed by 40 cycles of 95 C for 5 s and 60 C for 34 s on an ABI 7500 detection system (Applied Biosystems, Carlsbad, CA, USA). Nine biological replications were used for each group and three technical replications for each biological replicate (Table S1).

Assessment of gene expression
To reduce bias resulting from a single method, four software including geNorm, NormFinder, comparative DCt, and BestKeeper were used to analyze the expression stability of the candidate genes. For input data of geNorm and NormFinder, the raw cycle threshold (Ct) values were transformed to quantify using 2 -DCt (DCt, Ct value of each gene minus the lowest Ct value of the corresponding gene in different samples). The geNorm applet calculates the gene expression stability value (M value) of the reference genes during stepwise exclusion of the least stable reference gene (Vandesompele et al., 2002). The genes were ranked by increasing expression stability, ending with the two most stably expressed genes. Next, the optimal number of genes was identified according to a pairwise variation value (V value) between the candidate genes. The optimal number of reference genes was considered acceptable when the V value was <0.15. NormFinder calculates the stability value of each gene as well as a combination of two genes to analyze the expression stability of the candidate genes (Andersen, Jensen & Orntoft, 2004). The corresponding estimated intragroup and intergroup variations are also provided. Using the DCt method, we calculated the DCt values for the pairwise genes and assessed the expression stability using the standard deviations (SD) of DCt values for each reference gene (Pfaffl et al., 2004). BestKeeper calculates the expression stability based on coefficient variance (CV) and SD based on the untransformed Ct values. BestKeeper also ranks the candidate genes according to the stability (Silver et al., 2006).

Statistical analyses
All statistical analyses were performed using SPSS 21.0 software (IBM Corporation, Armonk, NY, USA). The data are presented as the mean ± SD unless indicated otherwise. All experiments were performed in at least three replicates.

Selection of the candidate reference genes
In this study, 15 commonly used candidate RNA reference genes (He et al., 2013;Kim et al., 2011;Klenke et al., 2016;Ruedrich et al., 2013;Taki, Abdel-Rahman & Zhang, 2014) were investigated in the hindgut of both normal and ETU-induced ARMS in rat fetuses. Detailed information including full names, NCBI accession numbers,  and position and PCR product lengths of these candidate genes are listed in Table 1. The sequences of forward and reverse primers used for RT-qPCR are listed in Table 2.

Expression profiles of the candidate reference genes
To preliminarily evaluate the expression of the candidate reference genes, the transcript abundances of these genes were estimated in normal, ARMs, and total samples. Only single peaks were found in RT-qPCR melting curves indicating that the specificity of the primers was good (Fig. 1). Based on the Ct values obtained from RT-qPCR, the expression profiles of the candidate genes in different samples were presented as the mean ± SD of the Cts in Fig. 2

Stability analysis based on geNorm
The expression stability of the candidate reference genes was further evaluated using four different methods: geNorm, NormFinder, comparative DCt, and BestKeeper. geNorm was used to calculate the stability values (M values) of each gene based on logarithmically transformed expression ratios and stepwise exclusion of the most unstable genes. As the M value decreased, the gene expression stability increased with improved ranking. Thus, the gene pairs with the lowest M values in the rank were most stably expressed. As shown in Fig. 3

Stability analysis based on NormFinder
NormFinder was used to calculate the stability of genes as well as intra-and inter-group variations. Similar to geNorm, the input data of NormFinder was also logarithmically transformed. The gene stability values calculated by NormFinder are presented in Fig. 4A. As the stability value decreased, the gene expression stability increased with high ranking order. Similar to geNorm, B2m (stability value ¼ 0.091), Pgk1 (stability value ¼ 0.058), and Rps16 (stability value ¼ 0.056) were the most unstable genes in NormFinder, further confirming their respective unstable expression. The four most stable genes were Rpl13a (stability value ¼ 0.031), Rps18 (stability value ¼ 0.032), Rplp1 (stability value ¼ 0.032), and Ywhaz (stability value ¼ 0.034). Slight differences occurred between geNorm and NormFinder. For example, Ywhaz ranked first in geNorm but fourth in NormFinder; and Sdha was included in four most stable genes in geNorm but not in NormFinder. Because of this discrepancy, various methods were necessary for the conjoint analysis while determining the optimal reference genes. The intra-and inter-group variations of each gene were also provided by NormFinder. As shown in Fig. 4B, Rps18 and Sdha also showed good stability with low inter-and intra-group variations of which variability values were close to zero.

Stability analysis based on comparative DCt
The comparative DCt method gives ranking based on the relative pair, and the average SD was used to evaluate the expression stability of the genes. The larger the average SD, the lower was the stability ranking. As shown in Table 3

Stability analysis based on BestKeeper
BestKeeper determines the expression stability based on the pairwise correlation of all the candidate genes. Due to input limitation of BestKeeper, we ranked the candidate genes and selected the 10 most stable according to ranking orders of previous three methods (geNorm, NormFinder, and comparative DCt). Seven genes showed strong correlation,   analyses were calculated using geometric means of the corresponding rankings of the top 10 genes. To avoid bias, only the genes common for at least three methods were chosen for further analysis, and Tbp was excluded from the comprehensive ranking because it appeared only in the top 10 gene list of NormFinder (Fig. 5). Rpl13a ranked first in every analysis, and so its comprehensive ranking order is still the highest, followed by Ywhaz (2.50), Rps18 (3.50), Sdha (4.25), Hmbs (5.50), Rplp1 (6.00), Actb (7.00), Ubc (7.25), Hprt1 (7.50), and Gapdh (8.00) ( Table 5).

DISCUSSION
Real-time polymerase chain reaction is one of the most commonly used methods in gene expression analysis. It provides a simultaneous estimation of various gene expressions in different samples. However, many factors can affect the results of PCR, including the selection of the reference genes (Ferguson et al., 2010;Guo et al., 2013). An ideal reference gene should be stably expressed in all the samples without being affected by species, tissues, and development (Harris, Reeves & Phillips, 2009;Mughal et al., 2018). In fact, researches have revealed that no reference genes can be universally used under all conditions (Thellin et al., 1999). Therefore, analyses of the candidate reference genes are essential for accurate results in PCR. In this study, we analyzed 15 candidate RNA reference genes in the hindgut of normal and ARMs in rat fetuses aiming to determine the optimal reference gene for RT-qPCR. In previous ARM-related studies, Gapdh and Actb were most commonly used as reference genes in RT-qPCR (Geng et al., 2017;Tang et al., 2018;Xiao et al., 2018). However, this study revealed that these two genes could be replaced by better candidates. In fact, six candidate reference genes were ranked before Gapdh and Actb in the expression stability including Rpl13a, Ywhaz, Rps18, Sdha, Hmbs, and Rplp1. Therefore, reference genes should be systematically evaluated, and only those that are relatively stable should be selected for further experiments.
While analyzing the gene expression stability, four different methods were mostly used, namely, geNorm, NormFinder, comparative DCt, and BestKeeper (Andersen, Jensen & Orntoft, 2004;Pfaffl et al., 2004;Silver et al., 2006;Vandesompele et al., 2002). In this study, nine (Actb, Gapdh, Rplp1, Ubc, Sdha, Hmbs, Rps18, Ywhaz, and Rpl13a) out of the top 10 stable candidate reference genes obtained by four methods overlapped. Out of the two remaining genes, Hprt1 appeared in of geNorm, BestKeeper, and Comparable DCt, while Tbp replaced Hprt1 in the top 10 list of NormFinder. However, discrepancies in the ranking orders by the four methods might have occurred because of the different computational principles followed. Thus, a comprehensive analysis of the multiple methods was further performed to reduce errors in screening optimal internal reference genes by a single method. This comprehensive analysis was broadly used in various studies to identify the optimal reference genes (He et al., 2013;Klenke et al., 2016;Taki, Abdel-Rahman & Zhang, 2014). A comprehensive ranking order was generated by sorting the arithmetic mean of each gene's rank order obtained from the four methods. The smaller the mean, the higher was the comprehensive ranking. A comprehensive analysis revealed that Rpl13a was the most stable reference gene followed by Ywhaz and Rps18. The commonly used Actb and Gapdh gene ranked seven and 10, respectively. These results indicated that Rpl13a was probably a better candidate for RT-qPCR normalization.
To the best of our knowledge, this is the first study to identify the optimal reference genes present in the hindgut of normal rat fetuses and those with ARMs. The results suggest the importance of systematically evaluating the expression stability of the candidate reference genes. However, there are some limitations to the present study. For example, RT-qPCR performed in this study was based on SYBR Green, and analyses between different experimental methods were not analyzed. The above problems are expected to be resolved in future studies.

CONCLUSIONS
In conclusion, we identified 15 candidate RNA reference genes in the hindgut of normal and ETU-induced ARMs in rat fetuses and found that Rpl13a, Ywhaz, and Rps18 were most stably expressed genes. This result provides valuable information for future gene expression studies in ETU-induced ARMs.