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.
Materials and Methods
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).
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.
|Genes||Full name||Accession No.||Position||Length (bp)||Amplification efficiency|
|Rps18||18S Ribosomal RNA||NM_213557||exon2/exon3||140||1.02|
|Ppia||Peptidylprolyl isomerase A||NM_017101||exon1/exon3||133||1.04|
|Hprt1||Hypoxanthine phosphoribosyl transferase 1||NM_012583||exon8/exon9||105||1.02|
|Pgk1||Phosphoglycerate kinase 1||NM_053291||exon8/exon9||104||1.03|
|Tbp||TATA binding protein||NM_001004198||exon4/exon5||123||1.01|
|Rps16||16S Ribosomal RNA||NM_001169146.1||exon1/exon3||147||1.05|
|Rpl13a||Ribosomal protein L13A||NM_173340||exon5/exon7||131||1.01|
|Rplp1||60S acidic ribosomal protein large PI||NM_001007604||exon2/exon3||87||0.98|
|Sdha||Succinate dehydrogenase complex, subunit A, flavoprotein (Fp)||NM_130428||exon9/exon10||105||0.99|
Accession No., NCBI accession number. PCR amplification efficiency = 10−1/slope − 1.
|Genes||Forward primer sequence (5′-3′)||Reverse primer sequence (5′-3′)|
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.
Avg SD, average standard deviation (SD) based on ΔCt of gene pairs for each reference gene; Ct, cycle threshold value; Rank, the stability of genes are ranked by Avg SD, and the higher the value, the lower the ranking.
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).
|Genes||N||GM (Ct)||AM (Ct)||Min (Ct)||Max (Ct)||SD (±Ct)||CV (% Ct)||r with BKI||p-value||Ranking order|
N, sample numbers; Ct, cycle threshold value; GM, geometric; AM, arithmetic; SD, standard deviation; CV, coefficient of variation; r, Pearson correlation coefficiency; BKI, BestKeeper Index; p-value, p-values associated with Pearson correlation.
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).
|Ranking||ΔCt||geNorm||NormFinder||BestKeeper||Comprehensive ranking (mean rank value)|
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 ΔCt, 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 ΔCt, 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.
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.
Anatomic part of the samples used for RT-qPCR.
The removed specimens are shown in the red dotted box. N, normal; A, anorectal malformations; U, urethra; R, rectum.