Validation of Reference Genes via qRT-PCR in Multiple Conditions in Brandt’s Voles, Lasiopodomys brandtii

Simple Summary This study validated the stability of the expression profiles of nine common candidate reference genes (Gapdh, Hprt1, β-actin, PPIA, Rpl13a, Tbp, Sdha, Hmbs, and B2M) using qRT-PCR in different tissues, developmental stages, and photoperiods. None of these genes were suitable as optimal reference genes at 4 weeks postnatal in different tissues. Under different developmental stages in the hypothalamus, B2M for males and Rpl13a for females were suitable as reference genes. Under different photoperiods in the hypothalamus, none of the selected genes were suitable as reference genes at 6 weeks postnatal, β-actin and PPIA were the optimal reference genes at 12 weeks postnatal, while Hprt1, β-actin, PPIA, Hmbs, and B2M were excellent reference genes at 24 weeks postnatal. Abstract The choice of optimal reference gene is challenging owing to the varied expression of reference genes in different organs, development stages, and experimental treatments. Brandt’s vole (Lasiopodomys brandtii) is an ideal animal to explore the regulatory mechanism of seasonal breeding, and many studies on this vole involve gene expression analysis using quantitative real-time polymerase chain reaction (qRT-PCR). In this study, we used the method of the coefficient of variation and the NormFinder algorithm to evaluate the performance of nine commonly used reference genes Gapdh, Hprt1, β-actin, PPIA, Rpl13a, Tbp, Sdha, Hmbs, and B2M using qRT-PCR in eight different tissues, five developmental stages, and three different photoperiods. We found that all nine genes were not uniformly expressed among different tissues. B2M and Rpl13a were the optimal reference genes for different postnatal development stages in the hypothalamus for males and females, respectively. Under different photoperiods in the hypothalamus, none of the selected genes were suitable as reference genes at 6 weeks postnatal; β-actin and PPIA were the optimal reference genes at 12 weeks postnatal; Hprt1, β-actin, PPIA, Hmbs, and B2M were excellent reference genes at 24 weeks postnatal. The present study provides a useful basis for selecting the appropriate reference gene in Lasiopodomys brandtii.


Introduction
The variation in gene expression level is a direct biomarker that can be used to capture individual responses to a changing environment. Quantitative real-time polymerase chain reaction (qRT-PCR) is a widely utilized method [1,2] with the advantage of easyaccessibility, fast processing, high accuracy, and high sensitivity in detecting the gene expression level [3]. In qRT-PCR, at least one reference gene is used as an internal control gene for the normalization of gene expression using 2 −∆∆Ct method [4], which can minimize the variations in RNA concentration and quantity, the amplification reaction, and a variety of treatments. An ideal reference gene is considered to be expressed at a constant level under all different conditions; such genes are often referred to as housekeeping genes. However, even housekeeping genes are differentially expressed across various tissues, developmental stages, and treatments [5]. Barber (2005) reported that a 15-fold difference in Gaphd mRNA copy number was observed between the highest-and lowest-expressing human tissues [6]. β-actin and Gapdh displayed a significantly variable expression in bronchoalveolar lavage fluid cells and endobronchial biopsy tissue between controls and patients [7,8]. The expression of β-actin by qRT-PCR showed a dose-dependent inhibition in matrigel treatment [9]. Therefore, a common tactic for selecting an optimal reference gene should be aimed at particular experimental conditions [10].
Currently, β-actin, Gapdh, PPIA, Rpl13a, Tbp, Sdha, B2M, and Hprt1 are commonly used as reference genes for qRT-PCR in humans [11], animals [12,13], and plants [3]. In seasonal breeding rodents, the photoperiod is considered as the most predictable indicator to mediate seasonal reproduction. The hypothalamus plays a vital role in photoperiodic response in seasonal breeders. Type 2 and type 3 iodothyronine deiodinase (Dio2 and Dio3) in the hypothalamus balance the local thyroid hormone levels to regulate the seasonal shifts of gonadal activity. Brandt's vole (Lasiopodomys brandtii) is an ideal animal to explore the regulatory mechanism of seasonal breeding, and many studies on this vole involve gene expression analyses using qRT-PCR. β-actin, Gapdh, and Hprt1 have been used for gene expression studies of seasonal breeding [14] and energy homeostasis [15] via qRT-PCR in Brandt's vole. However, suitable reference genes in different tissues, in the hypothalamus under various development stages, and in the hypothalamus under different photoperiods have not been identified for the normalization of the target gene stage's expression levels by qRT-PCR.
In the present study, we evaluated the stability of nine candidate reference genes in different tissues, developmental stages, and other photoperiod conditions. To further validate our results, type 2 iodothyronine deiodinase expression profiles were analyzed under different photoperiod conditions.

Sample Collection
Brandt's voles were collected from a laboratory colony maintained at the Institute of Plant Protection, Chinese Academy of Agricultural Sciences (CAAS). Food and filtered tap water were provided ad libitum; the cotton-nesting material was available in the cage. Ambient temperature and relative humidity were continuously held at 23 ± 2 • C and 50 ± 10%, respectively. Sampled tissues were immediately frozen and stored at −80 • C until RNA extraction. The sample information is listed in Table 1.

Selection of Candidate Reference Genes
According to previous studies in rodents, nine candidate reference genes were selected [16][17][18][19]. Gapdh and β-actin are classic reference genes in various species, and Hprt1, Rpl13a, Tbp, Sdha, Hmbs, PPIA, and B2M are high-frequency reference genes in rodents. The cDNAs of these genes were cloned according to the respective cDNAs of Microtus ochrogaster and used for further primer design. GeneBank Accession Numbers were

qRT-PCR
Primers for qRT-PCR were designed by the Oligo 7.0 (OLIGO team, Colorado Springs, CO, USA) software and synthesized by Sangon Biotech Company (Beijing, China). After the quality test of primers according to the standard curve in qRT-PCR using SYBR Green PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA), we used the BioMark™ HD System (Fluidigm Sciences Inc., South San Francisco, CA, USA) following Wang et al.'s method to assess mRNA levels of these candidate reference genes [14]. We used 48.48 Dynamic Array IFC (Fluidigm Sciences Inc., South San Francisco, CA, USA), which can load 48 samples and 48 assays. To exclude possible system errors, we performed a test of the same sample set using the same chip.

Statistical Analysis
The Shapiro-Wilk test was used for the normality test. The threshold cycle (Cq) value distribution of a candidate reference gene should be theoretically normal. Variations in the Cq value of each gene among different tissues, developmental stages, or photoperiods in individual assays were tested by one-way ANOVA or the Kruskal-Wallis test using SPSS version 19.0 (IBM Corp., Armonk, NY, USA) according to the normality test. p < 0.05 was considered statistically significant for all of the tests. The coefficient of variation (CV) was used to show the extent of the variability of candidates in all treatments in one sample set, with genes with lower CVs considered to be more stable. The NormFinder algorithm (Aarhus University Hospital, Aarhus, Denmark) provides a stability value that estimates the most stable genes based on the estimation of intra-and intergroup variations. Intragroup variation means variation among individuals of the same treatment, and intergroup variation means variation between treatments. NormFinder compares the expression stability of genes across tissues, developmental stages, and photoperiods. This software focuses on finding reference genes with lower intra-and inter-group variations. We used the intergroup variation of 0.5 recommended by NormFinder as the stability threshold [20]. The gene with the lowest stability value is considered to show a stably expressed pattern and vice versa.
The stability of gene expression was analyzed by both CV and the NormFinder algorithm [20]. The result of the Normfinder algorithm was obtained by R software version 3.6.1 (R core Team, Vienna, Austria) with codes (https://www.moma.dk/normfindersoftware. accessed on 14 June 2020).

Specificity and Amplification Efficiency of Primers
The single lane with the expected size on an agarose gel and the single peak in the melting curve of qRT-PCR, as shown in Figures S1 and S2, indicated the high specificity of each primer pair. The PCR efficiency ranged from 95.443 to 107.581%, the melting temperatures were all 62 • C, and the R-squared value ranged from 0.983 to 1 ( Table 2), indicating that these primers satisfied the standard requirements of qRT-PCR.

Expression Features in Different Tissues
We analyzed the intergroup variations and Coefficients of Variation (CV) of candidate genes and combined them with the NormFinder algorithm to screen out optimal reference genes under different treatments. After analyses of the intergroup Cq value variations, candidates with no significant difference in intergroup Cq values were further analyzed using the NormFinder algorithm. We also considered the results of the normality test, CV, and stability value provided by the NormFinder algorithm to determine the final suitable reference gene(s). In this assay, we analyzed the expression feature of candidate reference genes among different organs or tissues. Of nine candidates, seven genes, including Gapdh, β-actin, PPIA, Rpl13a, Sdha, Hmbs, and B2M, were not detected, with a significant intergroup difference (p > 0.05; Table 3), while the genes PPIA, Hmbs, and B2M displayed a significantly skewed distribution (p < 0.05; Table 3). The gene Sdha showed the lowest CV and the highest rank in the NormFinder algorithm (Table 4). These results indicated that Sdha was the best of the candidates when tissues were collected at postnatal day (PND) 4 w under the natural photoperiod condition. However, compared to the expected intergroup variation value (0.5) for an optimal reference gene using the NormFinder algorithm, the value of 1.95 indicated that Sdha could not be considered a good reference gene.

Expression Features at Different Developmental Stages
We tested the expression levels of nine candidates in male hypothalamus tissues collected at different developmental stages. Of the nine candidates, only B2M did not exhibit a significant intergroup difference (p > 0.05; Table 5). This indicates that most candidates were variably expressed at different developmental stages in the male hypothalamus and were not suitable for use as reference genes. B2M displayed a normal distribution and a low value of CV and could be selected as a reference gene under these circumstances. Here, we did not test its performance using NormFinder, as it requires at least three genes to perform the test.
In female samples, we detected three of nine candidates, including PPIA, Rpl13a, and B2M, that displayed no significant intergroup difference (p > 0.05; Table 6). Of the three genes, B2M displayed the lowest CV and the worst intergroup variation by the NormFinder algorithm, and Rpl13a displayed a similar CV and the highest rank by the NormFinder algorithm (Tables 6 and 7). This indicates that Rpl13a is an optimal reference gene in females in different developmental stages.

Expression Features under Different Photoperiod Conditions
We analyzed the expression features of candidate genes in hypothalamus tissues at the developmental stage of PND 6w, which experienced different photoperiod conditions. Of nine candidates, only Tbp displayed no significant intergroup difference (p > 0.05; Table 8). These results indicate that most candidate genes are sensitive to photoperiod conditions at the developmental stage of PND 6w. The Shapiro-Wilk test indicated a significantly skewed distribution of Tbp. This result indicates that even Tbp is not a good choice as a reference gene under different photoperiod conditions, although it displays a relatively small intergroup variation. At the developmental stage of PND 12w, all candidate genes displayed no significant intergroup difference under different photoperiod conditions and no significantly skewed distribution (Table 9). This indicates that these genes are not sensitive to photoperiod conditions when Brandt's voles develop to PND 12w. PPIA and β-actin displayed a lower intergroup variation than the threshold of 0.5 using the Normfinder algorithm (Table 10). This indicates that these two genes can be selected as the optimal reference genes under different photoperiod conditions. At the developmental stage of PND 24w, all candidate genes displayed no significant intergroup difference under different photoperiod conditions and no significantly skewed distribution (Table 11). By the NormFinder algorithm, five of nine candidates, including Hprt1, β-actin, PPIA, Hmbs, and B2M, displayed a lower intergroup variation than the  (Table 12). These genes could be selected as the reference genes in future studies. Compared to the developmental stage of PND 12w, more candidate genes were expressed stably under different photoperiod conditions when Brandt's voles developed to PND 24w.

Validation of the Selected Reference Genes under Different Photoperiod Conditions
To validate the stability of the candidate reference genes, the mRNA relative expression levels of a target gene Dio2 were normalized by the top two most stable and least stable candidates in their individual developmental stages. Under different photoperiods, no optimal reference genes were selected at PND 6. At PND 12 ( Figure 1A), there was no significant difference between Dio2 expression normalization by the top two most stable reference genes, PPIA and β-actin, in long-day photoperiod (LD), natural-day photoperiod (ND), or short-day photoperiod (SD) condition. One-way ANOVA indicated a significant difference between PPIA, β-actin, Sdha, and B2M under the ND condition (p = 0.015), and the result normalized using Sdha was significantly lower than B2M under ND condition (independent t-test, p = 0.018). Although no significant difference was found under the SD condition, the fold ratio of the mean of the expression level of the result between using PPIA and B2M was up to 2.35. Similarly, at PND 24w ( Figure 1B), there was no significant difference between the top two most stable reference genes, Hmbs and Hprt1, in LD, ND, or SD condition. No significant difference was found analyzed by one-way ANOVA, whereas the fold ratio of the mean of the expression level of the result between using Hmbs and Rpl13a was up to 1.76. These results indicated the qPCR results normalized by the top two most stable reference genes were beyond those normalized by the top two least stable reference genes under different photoperiod conditions. whereas the fold ratio of the mean of the expression level of the result between using Hmbs and Rpl13a was up to 1.76. These results indicated the qPCR results normalized by the top two most stable reference genes were beyond those normalized by the top two least stable reference genes under different photoperiod conditions.

Discussion
In this present study, we evaluated the expression levels of nine genes that are usually selected as the reference genes in qRT-PCR. We found that all the tested genes displayed a higher intergroup difference in various tissues at PND 4w. B2M and Rpl13a were the optimal reference genes for different postnatal development stages in the hypothalamus of males and females, respectively. Under different photoperiods in the hypothalamus, none of the selected genes were suitable as reference genes at PND 6w; β-actin and PPIA were the optimal reference genes in the hypothalamus after PND 12w, more optimal reference genes (Hprt1, β-actin, PPIA, Hmbs, and B2M) were found at PND 24w.
We analyzed the data using the coefficient of variation (CV) and NormFinder, which assessed the candidate reference genes' stability using different algorithms. By using CV and intergroup statistical analysis, the candidates with a normal distribution, no significant intergroup variation, and a lower CV were found to be relatively stable, which can be further analyzed using NormFinder. The NormFinder algorithm combines intragroup and intergroup variations to calculate the candidate reference genes' stability value for validation. The reference gene with the lowest intergroup variation and stability value is the most stable. Only the candidates who pass the evaluation using these two algorithms can be selected as an optimal reference gene.
None of these genes were suitable as optimal reference genes in different tissues at PND 4w. This may be due to the fact that most dynamic changes in gene expression occur before puberty, during which different tissues have different gene expression profiles and types [21][22][23]. For example, Tbp was highly expressed in the testis compared with other organs at PND 4w in our study because it is involved in the transcription of the androgen receptor [24]. Androgen receptors are highly expressed in the testis during puberty [25], and PND 4w is around puberty in Brandt's voles, which reach sex maturation after 4-6 weeks [26]. These results remind us that other candidate reference genes besides the selected candidates in this study need to be selected for validation.
We next investigated the stability of the nine selected genes in the hypothalamus under different developmental stages and photoperiod conditions. In this study, B2M in males and Rpl13a in females were not developmentally dynamic genes and had housekeeping functions in the postnatal development stages. Although B2M and Rpl13a are not classic reference genes, they are also listed as suitable reference genes in other studies. For instance, B2M is the most stable reference gene in the developmental stage of the halfsmooth tongue sole (Cynoglossus semilaevis) [27], and Rpl13a is the most stable reference

Discussion
In this present study, we evaluated the expression levels of nine genes that are usually selected as the reference genes in qRT-PCR. We found that all the tested genes displayed a higher intergroup difference in various tissues at PND 4w. B2M and Rpl13a were the optimal reference genes for different postnatal development stages in the hypothalamus of males and females, respectively. Under different photoperiods in the hypothalamus, none of the selected genes were suitable as reference genes at PND 6w; β-actin and PPIA were the optimal reference genes in the hypothalamus after PND 12w, more optimal reference genes (Hprt1, β-actin, PPIA, Hmbs, and B2M) were found at PND 24w.
We analyzed the data using the coefficient of variation (CV) and NormFinder, which assessed the candidate reference genes' stability using different algorithms. By using CV and intergroup statistical analysis, the candidates with a normal distribution, no significant intergroup variation, and a lower CV were found to be relatively stable, which can be further analyzed using NormFinder. The NormFinder algorithm combines intragroup and intergroup variations to calculate the candidate reference genes' stability value for validation. The reference gene with the lowest intergroup variation and stability value is the most stable. Only the candidates who pass the evaluation using these two algorithms can be selected as an optimal reference gene.
None of these genes were suitable as optimal reference genes in different tissues at PND 4w. This may be due to the fact that most dynamic changes in gene expression occur before puberty, during which different tissues have different gene expression profiles and types [21][22][23]. For example, Tbp was highly expressed in the testis compared with other organs at PND 4w in our study because it is involved in the transcription of the androgen receptor [24]. Androgen receptors are highly expressed in the testis during puberty [25], and PND 4w is around puberty in Brandt's voles, which reach sex maturation after 4-6 weeks [26]. These results remind us that other candidate reference genes besides the selected candidates in this study need to be selected for validation.
We next investigated the stability of the nine selected genes in the hypothalamus under different developmental stages and photoperiod conditions. In this study, B2M in males and Rpl13a in females were not developmentally dynamic genes and had housekeeping functions in the postnatal development stages. Although B2M and Rpl13a are not classic reference genes, they are also listed as suitable reference genes in other studies. For instance, B2M is the most stable reference gene in the developmental stage of the halfsmooth tongue sole (Cynoglossus semilaevis) [27], and Rpl13a is the most stable reference gene in the developmental stage of Diaphania caesalis (Lepidoptera Pyralidae) [28]. Under different photoperiods, the expression of all nine genes was affected by photoperiods for males going through puberty, while these candidates' stability increased along with the development. It was supported by the increased number of suitable reference genes from PND 12w to PND 24w. We verified these results using a photoperiod-sensitive gene, Dio2.