Validation of endogenous reference genes in rat cerebral cortex for RT-qPCR analyses in developmental toxicity studies

Introduction

Exposure to thyroid hormone (TH) disruptors during pregnancy can have detrimental consequences for the child’s brain development. In the developing brain, THs exert their action by regulating gene transcription required for normal development and function. If TH levels or action is disturbed, it can lead to dysregulated gene expression and ultimately cause structural and functional changes to the brain (Bernal, 2017; Gilbert & Zoeller, 2011). Worryingly, numerous compounds that humans are exposed to can also disturb TH levels or action. Understanding how these compounds can cause adverse health effects is thus of great importance and will help us conceive new testing strategies designed to capture harmful compounds before they potentially are used in ways that may harm the developing child.

Various effect biomarkers are currently being explored in an effort to devise more sensitive and reliable testing strategies to detect adverse effects on brain developmental caused by developmental TH disruption. For instance, measurements of TH levels have recently been included in several OECD test guidelines used for toxicity testing of environmental chemicals, such as TG 421 (OECD, 2016a), TG 422 (OECD, 2016b) and TG 414 (OECD, 2018). Adverse outcome pathways (AOPs) are also being developed to improve our understanding of how different molecular initiating events can cause TH insufficiency in the brain and lead to adverse neurodevelopmental effects. One knowledge gap in these AOPs is the specific relationship between TH-dependent gene expression and subsequent effects on the developing nervous system (Crofton et al., 2017). Thus, future studies should aim to characterize this relationship at the molecular level.

Despite great advances in different transcriptomics approaches (Darde, Chalmel & Svingen, 2018), quantitative RT-PCR (RT-qPCR) remains a popular method for determining relative gene expression levels, especially when only a few target genes are examined. RT-qPCR is considered a reliable and cheap way of quickly determining relative transcript abundance without having to perform more complex bioinformatics post-run analyses. Most often, RT-qPCR results are simply converted into relative fold expression data by the comparative Ct method (Schmittgen & Livak, 2008). Although theoretically a sound method, it does rely on certain assumptions that must be adhered to if the data are to be true representations of the transcriptional landscape of the cells or tissues. These assumptions include stable RNA integrity across samples, correct primer design, amplification efficiencies close to 100% unless absolute quantification methodologies are implemented, and correct normalization of input RNA levels (Nolan, Hands & Bustin, 2006). These are all critical parameters, but normalization continues to be an overlooked issue, especially when using endogenous reference genes (RGs) for normalization (Piller, Decosterd & Suter, 2013).

Before exploring the utility of using changes to gene expression levels in the rodent brain as an adverse effect readout for TH disruption, it is imperative that we first identify RGs that are relatively unchanged by exposure to putative TH disrupting chemicals, not least if relative quantification methods such as RT-qPCR are to be used. Although there are studies that have examined stability of common RGs in rat brains at various developmental stages or following insult (Al-Bader & Al-Sarraf, 2005; Cook et al., 2009; Julian et al., 2014; Swijsen et al., 2012; Zhou et al., 2016), a systematic evaluation of suitable RGs in rat brains following developmental TH disruption is lacking. In this study, we have tested the stability of 10 common RGs in the cerebral cortex of rat offspring after exposure to the TH synthesis inhibitor drug propylthiouracil (PTU), as well as two putative TH disruptors: the UV filter octyl methoxycinnamate (OMC) and the flame retardant pentabromodiphenyl ether mixture, DE-71.

Materials and Methods

Results

In RT-qPCR experiments, it is advisable to first normalize input RNA levels by using the same amount of RNA for cDNA synthesis across all samples to be included in the same analysis. We used 500 ng of total RNA for each sample and synthesized cDNA from the same master mix. Following RT-qPCR amplification, we acquired the raw mean Ct values from technical duplicate reactions and did a first-round assessment of relative mRNA abundance and stability across samples for 10 putative endogenous RGs. Looking across all the 28 samples representing four different exposure groups, the SD did not exceed 0.42 cycles, and with a maximum SD of 0.5 in the PTU group (Table 2). This indicates a high level of expression stability for all included RGs in the cortex of young rats.

We used three different algorithms to assess the relative expression stability of the 10 RGs: geNorm (Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Ørntoft, 2004) and BestKeeper (Pfaffl et al., 2004). With some exceptions, the stability ranking was similar between the three algorithms, all scoring Rps18, Rpl13a and Sdha in the top five and Gapdh, B2m and Actb in the bottom five (Table 3). Notably, all 10 RGs were deemed suitable (“stable enough”) for use as single endogenous normalizing genes when performing RT-qPCR analysis with the comparative Ct method on rat cortex tissue under the experimental conditions as described in this study. For relatively homogenous tissues, the cut-off M-value for geNorm is <0.5, whereas the worst performer tested herein, Actb, had an M-value of 0.023. Likewise, NormFinder recommends a cut-off S-value of <0.5 and the highest S-values obtained in our experiment was 0.081 for B2m. Finally, BestKeeper recommends that RGs should have a SD value <1.0 to be deemed suitable for normalization and herein the highest SD value was 0.32 for Actb.

It is common to use only one endogenous RG to normalize RT-qPCR data, although this has a higher risk of resulting in erroneous normalized expression data. Using two or more RGs for normalization minimizes this problem and is recommended, provided they are all stably expressed across the examined samples (Vandesompele et al., 2002). The geNorm and NormFinder algorithms also provide suggestions for best combinations of RGs, which not necessarily comprise those that are ranked at the very top of the list. Across the samples tested in this study, geNorm suggests using a combination of Rps18 and Rpl13a, whereas NormFinder suggests combining Rps18 and Ubc (Fig. 1). Notably, Rps18 was ranked second best by all three algorithms as single-use RG and was included in both RG pairs.

To test the performance of the 10 RGs, we performed relative fold expression analysis on two genes, Pvalb and Dio2, as they could potentially be dysregulated after exposure to PTU and other chemicals. Both genes were normalized by the comparative Ct method using all of the 10 RGs individually (Fig. 2). Pvalb was significantly downregulated in the cortex of PTU-exposed offspring (around threefold) regardless of which of the 10 RGs that was used for normalization of input RNA. Dio2 was unchanged in the exposed tissue relative to controls when performing the analysis with most of the RGs; however, it was significantly up-regulated (up to about 1.4-fold) when Rpl13a, Sdha or Tbp was used for normalization. When normalizing data using the best pair combinations suggested by geNorm or NormFinder (Fig. 3), Pvalb was again significantly downregulated (about threefold), whereas Dio2 was unchanged between control and PTU-exposed cortex tissue, albeit the mean values showed a trend toward being upregulated in the exposed group. Comparable expression patterns were observed between the two pairs of RGs for OMC and DE-71 exposed groups.

Discussion

Quantitative RT-PCR remains a popular method for analyzing gene expression levels in cells and tissues. It has developed into a fast, cheap and reliable method to quantify the abundance of RNA species. However, data analysis relies heavily on specific assumptions, not least the use of stable internal normalizers for the popular comparative Ct method (Nolan, Hands & Bustin, 2006). Since using endogenous RGs is the most common way of normalizing inter-sample variations, there has been an increasing focus on testing RG stability across different tissues under various experimental conditions; and this has proven to be extremely important, as many “historical housekeeping genes” have in many instances been shown to be highly unstable (Piller, Decosterd & Suter, 2013).

In this study, we tested 10 common endogenous RGs for their stability in cerebral cortex of young male rats after in utero exposure to three TH-disrupting chemicals: PTU, OMC and DE-71. Expression stability was tested across the three exposure groups in addition to control animals. Following PCR cycling, the relative expression stability was remarkably good for all of the 10 RGs. According to three different algorithms designed to calculate relative expression stability—geNorm, NormFinder and BestKeeper—all of the 10 RGs obtained stability scores well inside recommended cut-off values, in essence predicting all of the RGs to be suitable for RT-qPCR normalization. This is in contrast to previous studies testing across different rat tissues (Svingen et al., 2015) or brain tissue from rats (Cook et al., 2009; Yang et al., 2008; Zhou et al., 2016), mice (Crans et al., 2019) or humans (Rydbirk et al., 2016). Notably, brain tissue from a rat sleep apnea model displayed high stability of all tested RGs except 18S rRNA (Julian et al., 2014), a stability also seen in rat brains across fetal developmental stages (Al-Bader & Al-Sarraf, 2005), which agrees with our findings of relative high stability among RGs in this tissue (Julian et al., 2014). There are, however, studies showing unstable expression of common RGs depending on experimental parameters and what brain regions are examined (Cook et al., 2009; Swijsen et al., 2012; Yang et al., 2008; Zhou et al., 2016).

geNorm suggests that the stability value M should be lower than 0.5 for homogenous tissues in order to be used as endogenous normalizing gene. In our study, Actb was the worst performer yet obtained an M-value of 0.023, which is suggestive of high stability between groups. Similar scores were obtained with the other two algorithms. Based on these scores, it can be concluded that all the 10 RGs tested herein are suitable for normalization of rat cortex under the described experimental conditions. However, we would still recommend using at least two RGs when performing RT-qPCR analyses, also in accordance with previous recommendations (Vandesompele et al., 2002). This point is particularly apparent in view of the data obtained for the proof-of-principle gene expression profiles that we carried out for Pvalb and Dio2.

To scrutinize the relative expression data obtained by using the comparative Ct method depending on which endogenous RG that was used, we performed RT-qPCR experiments on the same brain tissues for Pvalb and Dio2. These two genes were chosen based on previous reports showing dysregulation in rat brain following exposure to TH disrupting chemicals. For instance, to compensate for low brain T3 Dio2 expression levels may increase for increased T4 deiodination to maintain T3 (Morse et al., 1993; Peeters et al., 2001; Sharlin et al., 2010), whereas Pvalb may be a suitable biomarker for TH disruption in the brain; Pvalb expression has been correlated to serum T4 levels (O’Shaughnessy et al., 2018). Correspondingly, we found that Pvalb was significantly downregulated to about threefold in the PTU-exposed rats compared to controls. This result is unequivocal irrespective of which of the 10 RGs was used for normalization, albeit with small variations in relative means and p-values. Dio2, on the other hand, was significantly upregulated when normalizing data with either Rpl13a, Sdha or Tbp, whereas the tending upregulation was not statistically significant if the data were normalized with any of the remaining seven RGs; nor with the recommended pairs of RGs: Rpl13a/Rps18 or Rps18/Ubc. This can seem counter-intuitive since all of the RGs should be stably expressed across the samples in question, but the answer may lie more with the uncertainty of determining significant changes in gene expression if they are of a smaller magnitude.

When performing RT-qPCR experiments, they are typically run in duplicate or triplicate reactions, so-called technical replicates. It is common to allow for variation of 0.5–1.0 Ct within technical replicates, which in itself will impact on the relative sensitivity of the assay. When transforming data, as is done with the comparative Ct method, a difference of 0.5 Ct (comparable to delta-Ct between samples) will equal about 1.4-fold difference in relative expression. Thus, it can become problematic to reliably determine changes in relative expression less than 50% up or down regardless of statistical p-value. This sensitivity parameter should be considered when reporting data with small changes to relative fold expression, as for instance is the case with the Dio2 data presented in Fig. 2. Notably, Dio2 was not statistically significantly altered relative to control when normalizing the data with the recommended pairs of RGs.

For future studies looking at the effects on gene expression in cerebral cortex of rat offspring, it is not certain that the RGs evaluated to be stably expressed herein remains stable under new experimental conditions. But by closely monitoring RT-qPCR performance, and in particular any changes to mean Ct values for RGs, it is feasible to detect when the RGs are affected by the experimental conditions. In such cases, alternative RGs should be verified to be unaffected by exposure and used instead. By including two to three RGs for the experiments, the chances of correcting for smaller changes to RG regulation is greatly improved contra only using one RG for normalization.

Conclusions

In conclusion, we recommend using either of the combinations Rps18/Rpl13a or Rps18/Ubc for RT-qPCR experiments on rat cortex following in utero exposure to TH disruptors.

Supplemental Information