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

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 μl using PrimeScript RT reagent kit (Takara Biotechnology Co. Ltd., Dalian, China) according to the manufacturer’s instructions.

Selection of the candidate reference genes

A total of 15 commonly used candidate genes (Rps18, Actb, B2m, Gapdh, Ppia, Hprt1, Pgk1, Ywhaz, Tbp, Ubc, Rps16, Rpl13a, Rplp1, Sdha, and Hmbs) were assessed to identify the most stably expressed and optimal reference genes in studies of ARMs. All the primers were designed and synthesized by Sangon (Sangon Biotech Co. Ltd., Shanghai, China).

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 μl cDNA template, 0.8 μl of each forward and reverse primers (10 nM), 6 μl RNase-free water, 0.4 μl ROX II, and 10 μl SYBR Green premix Ex Taq II (Tli RNaseH Plus) in a final volume of 20 μl. 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 ΔCt, 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^−ΔCt (ΔCt, 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 ΔCt method, we calculated the ΔCt values for the pairwise genes and assessed the expression stability using the standard deviations (SD) of ΔCt 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. Generally, the mean Ct values ranged from 14 to 22, mainly 14–17. Among all, Rplp1 was most abundantly expressed with lowest mean Ct value (14.07 ± 0.36), followed by Rpl13a (14.40 ± 0.41), Rps16 (14.87 ± 0.65), Rps18 (14.89 ± 0.29), Ppia (15.10 ± 0.53), Actb (15.10 ± 0.56), Gapdh (16.06 ± 0.44), Ubc (16.57 ± 0.25), Ywhaz (17.29 ± 0.43), Pgk1 (19.01 ± 0.63), B2m (19.20 ± 0.67), Hprt1 (21.36 ± 0.39), Hmbs (21.38 ± 0.40), Tbp (21.64 ± 0.62), and Sdha (21.89 ± 0.48). However, the stability of the candidate reference genes needs to be further explored.

Stability analysis based on geNorm

The expression stability of the candidate reference genes was further evaluated using four different methods: geNorm, NormFinder, comparative ΔCt, 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, B2m (M value = 0.533), Pgk1 (M value = 0.494), and Rps16 (M value = 0.459) were the most unstable genes, while Rpl13a/Ywhaz (M value = 0.266), Sdha (M value = 0.286), and Rps18 (M value = 0.311) were the most stable ones. Among these, Ywhaz/pl13a was the most stable gene pair.

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 ΔCt

The comparative ΔCt 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, Rpl13a (0.423) was top-ranked, followed by Ywhaz (0.437), Rps18 (0.439), Hmbs (0.454), Sdha (0.458), Ubc (0.461), Rplp1 (0.471), Gapdh (0.490), Actb (0.517), and Hprt1 (0.527). In contrast, B2m (0.785), Pgk1 (0.726), Rps16 (0.703), Tbp (0.556), and Ppia (0.542) genes were ranked the last.

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 ΔCt). Seven genes showed strong correlation, including Rpl13a (r = 0.921), Actb (r = 0.892), Ywhaz (r = 0.891), Sdha (r = 0.876), Rps18 (r = 0.786), Hmbs (r = 0.726), Rplp1 (r = 0.707), and Gapdh (r = 0.696) (Table 4). The other genes that showed moderate correlation were Ubc (r = 0.55) and Hprt1 (r = 0.51).

Comprehensive stability analysis of the candidate reference genes

Due to different algorithms, different methods used in this study might have led to different ranking order. To avoid bias resulting from a single method and to obtain a more reliable result of gene stability, the comprehensive ranking orders based on four 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).