Comparison of reference gene expression stability in mouse skeletal muscle via five algorithms

Real-time quantitative PCR (RT-qPCR) is a widely applied technique for relative quantification of gene expression. In this context, the selection of a suitable reference gene (RG) is an essential step for obtaining reliable and biologically relevant RT-qPCR results. The present study aimed to determine the expression stability of commonly used RGs in mouse skeletal muscle tissue. The expression pattern of eight RGs (ACTB, GAPDH, HPRT, YWHAZ, B2M, PPIA, TUBA and 18S) were evaluated by RT-qPCR in different sample groups classified based on genetic background, muscle tissue type, and growth stage, as well as in a C2C12 myoblast cell line model. Five computational programs were included in the study (comparative ΔCq value, NormFinder, BestKeeper, geNorm, RefFinder) to evaluate the expression stability of RGs. Furthermore, the normalization effects of RGs in soleus (SOL) and gastrocnemius (GAS) muscle tissue were evaluated. Collectively, ACTB, HPRT and YWHAZ were shown to be the most stable RGs, while GADPH and 18S were the least stable. Therefore, the combined use of ACTB, HPRT and YWHAZ is recommended for the normalization of gene expression results in experiments with murine skeletal muscle. The results discussed herein provide a foundation for gene expression analysis by RT-qPCR in mammalian skeletal muscle.


INTRODUCTION
In most mammals, nearly 40% of body weight is composed of skeletal muscle, which is a complex and heterogeneous tissue (Bentzinger, Wang & Rudnicki, 2012). Skeletal muscle plays an important role in locomotion and metabolic regulation of the body. Several physiological processes occur in skeletal muscle and involve changes in gene expression, such as myogenesis (Wang et al., 2019), atrophy (Yin et al., 2021) and regeneration (Tidball, excludes unstable RGs to enable individual calculation of the expression stability (M value) of RGs (Vandesompele et al., 2002), and inadequate RGs are then discarded until two optimal RGs are found. Moreover, Pfaffl et al. (2004) developed BestKeeper to identify the most stable RGs based on the standard deviation (SD) and coefficient of variance (CV) of Cq value of RGs. In addition, Silver et al. (2006) proposed a comparative Ct value algorithm that enabled the calculation of the SD of the difference between Cq values ( Cq) of two genes across all samples, and RGs are then ranked based on the average of SD (Andersen, Jensen & Øntoft, 2004). However, RefFinder is the most recent online tool for the analysis of the selection of RGs (Xie et al., 2012), which comprehensively considered multiple algorithms to provide a ranking of RGs. These are free tools that are readily accessible online, thus providing a great convenience for the analysis of RGs.
Previous studies have reported the selection of RGs in muscle tissue (Nakao et al., 2015;Hildyard, Finch & Wells, 2019). In the present study, the expression stability of eight RGs in mouse skeletal muscle tissue was compared using five different algorithms, i.e., comparative Cq value ( Cq), geNorm, NormFinder, BestKeeper and RefFinder. The selected RGs were ACTB, GAPDH, 18S and HPRT which are commonly used in RT-qPCR experiments in skeletal muscle; in addition, YWHAZ, PPIA, TUBA and B2M were included in the experiments based on previous studies (Masilamani, Loiselle & Sutherland, 2014;Nakao et al., 2015;Niu et al., 2016;Hildyard, Finch & Wells, 2019). Differences in the stability of RGs in mice of different genetic backgrounds, namely inbred C57BL/6 and outbred ICR was compared. Furthermore, differences in the stability of RGs in skeletal muscle in different body sections and growth stages were identified. In addition, murine C2C12 myoblast cell line was included as an in vitro model of skeletal muscle. Finally, the normalization effect of optimized combination of RGs were determined to RT-qPCR analysis of the following genes: troponin I1 (TNNI1), troponin I2 (TNNI2), troponin C1 (TNNC1), troponin C2 (TNNC2) in soleus (SOL) and gastrocnemius (GAS) muscle tissue. Taken together, the results discussed herein provided an optimized protocol for the selection of RGs to be applied in future RT-qPCR experiments in skeletal muscle samples.

Ethics statement
All procedures described herein involving animals were approved by the Animal Ethical and Welfare Committee of Sichuan Agricultural University, China (approval number: 20220207), according to the guidelines of the Regulations on the Management of Laboratory Animal License (Ministry of Science and Technology, China, 2004).

Sample collection
C57BL/6 and ICR male mice were purchased from Chengdu Dashuo Experimental Animal Co., Ltd. (Chengdu, Sichuan, China). All animals were housed individually in plastic cages and maintained in a dedicated animal room at 22 • C ± 3 • C and 40% humidity under a natural light cycle. During the experiments, the animals were housed on shavings and provided with free sufficient food and water. All mice were treated humanely and killed by ether asphyxiation to collect muscle samples. C2C12 cells line were purchased from NATIONAL INFRASTRUCTURE OF CELL-LINE RESOURCE (NICR) (http://www.cellresource.cn). C2C12 myoblasts in natural proliferation without any treatment were collected based on a previous study conducted by our group (Gan et al., 2022).
In order to investigate the effects of different factors on expression stability of RGs, RT-qPCR analysis were conducted on samples of four distinct groups: (i) group A included soleus muscle tissue (SOL) samples from 8-week-year-old mice of two mouse lineages (C57BL/6 mouse, n = 8; ICR mouse, n = 6); (ii) group B included samples from tibialis anterior (TA), longissimus dorsi (LD), gastrocnemius (GAS), SOL from 8-week-year-old ICR mice (n = 6); (iii) group C included samples of mice in different growth stages (based on days of age), namely immature ICR mice in Lactation (n = 12) and Adulthood ICR mice (all samples in group B, n = 24); (iv) group D was based on sample group, and differences between the expression stability of RGs in skeletal muscle tissue (SM) (above muscle tissue samples, n = 44) and C2C12 myoblasts line (n = 24) were explored.

RNA extraction and reverse transcription
Total RNA from skeletal muscle and C2C12 myoblasts were extracted using the TRIzol kit (TaKaRa, Dalian, China) according to the manufacturer's instructions. RNA concentration and quality were determined using NanoDrop 2000 (Thermo Scientific, Waltham, MA, USA) (Table S3). PrimeScript TM RT reagent kit with gDNA Eraser (TaKaRa) was used to remove genomic DNA from 1,000 ng total RNA samples subsequently used for reverse transcription according to the manufacturer's instructions. Reverse transcription was conducted in the following conditions: 42 • C for 2 min; then maintained at 4 • C; followed by reverse transcription at 37 • C for 15 min, then heated to 85 • C for 5 s, and maintained at 4 • C. All samples were stored at −80 • C until further analysis.

RT-qPCR
RT-qPCR was performed using the TB Green R Premix Ex Taq TM II (Tli RNaseH Plus) kit (TaKaRa, Code No.RR820Q) in a CFX96 real time PCR detection system (Bio-Rad, Richmond, CA, USA) according to the manufacturer's instructions. Amplification program was as follows: pre-denaturation at 95 • C for 5 min, followed by 40 cycles at 95 • C for 10 s, then at 60 • C for 30 s according to the manufacturer's instructions. For each sample, RT-qPCR was performed in duplicates to enable calculation of the average Cq value for each gene. Cq value in each sample and sample information are listed in Table S1. Primer sequences were designed using the Primer-BLAST software (https://www.ncbi.nlm.nih.gov/tools/primer-blast) (Table S2). Melting curve analysis was carried out to verify primer specificity (Fig. S1). Primer pairs were validated by making standard curves using a cDNA dilution series to assess amplification efficiency. Primer amplification efficiencies are shown in Table S2. All selected primer pairs displayed between 90-110% amplification efficiency. Relative gene expression of target genes was calculated using the 2 − Cq method according to Livak method (Livak & Schmittgen, 2001).

Statistical analysis methods
Spearman correlation analysis and data visualization were performed using GraphPad Prism 8 software (GraphPad Inc., San Diego, CA). Comparison of intergroup means was performed using the non-parametric Mann-Whitney test in GraphPad Prism 8 software.
Differences were reported as statistically significant when p values were <0.01. showed a lower extent of dispersion, while GADPH, TUBA, 18S showed a contrary pattern hence indicating that their expression pattern was unstable. Therefore, expression stability of RGs and the number of selected appropriate RGs were further explored. Then, we divided the samples into four groups to analyze the stability of RGs. Group A includes samples from different mouse strains. Group B includes samples from different skeletal muscular tissue types. Group C includes samples from different growth stages. The samples in group D were divided into skeletal muscle tissue or C2C12 myoblasts.

Comparative Cq value ( Cq)
According to the results of Cq analysis, stability values of RGs in each sample group are shown in Table 2. The top two best RGs based on Cq results for each sample group was as follows: (i)

NormFinder
NormFinder is an application based on Microsoft Excel that considers intragroup and intergroup variations to calculate the stability of RGs. Table 3   Notes.
Values into the parenthesis refer to ranking of stability value in each group. Where 1 is given highest priority and 8 is the lowest priority. The last row is the geometric mean of the reference gene ranking in all groups. The top four RGs were marked in bold.

Notes.
Values into the parenthesis refer to ranking of stability value in each group. Where 1 is given highest priority and 8 is the lowest priority. The last row is the geometric mean of the reference gene ranking in all groups. The top four RGs were marked in bold.  Table 4. We then reanalyzed the Pearson coefficients of correlation between RGs with the BestKeeper index after excluding tow RGs with top SD for each group. ACTB and GAPDH were ranked as the most stable RGs in sample group of different genetic background. In TA, LD, GAS and SOL samples, top ranked RGs were YWHA, PPIA, ACTB and GAPDH, respectively. YWHAZ and HPRT were the most stable gene in samples of mice in the two growth stage (lactation: SD = 1.02, r = 0.984; adulthood: SD = 0.314, r = 0.888). In skeletal muscle tissue and C2C12 myoblasts, HPRT (SD = 0.942, r = 0.957) and YWHAZ (SD = 0.856, r = 0.919) were the optimal RGs, respectively. The performance of RGs based on Best-Keeper analysis was compared with Cq and NormFinder results, which revealed that low-ranking RGs were similar.
geNorm geNorm evaluated the stability of M values based on stepwise exclusion of the least stable RGs, and results are summarized in Fig. 2. RGs with high stability had lower M value. Table 5 shows the ranking of stability of RGs based on geNorm analysis. RGs with the highest suitability were as follows: (i) group A: GAPDH and YWHAZ in ICR, and HPRT and PPIA in C57BL/6; (ii) group B: TUBA and B2M in TA, HPRT and ACTB in LD, HPRT and PPIA in GAS, and GAPDH and YWHAZ in SOL; (iii) group C: PPIA and ACTB in mice in lactation, and HPRT and B2M in adult mice; (iv) group D: TUBA and YWHAZ in skeletal muscle tissue, and PPIA and ACTB in C2C12 myoblasts. The ranking of RGs obtained with geNorm was similar compared to BestKeeper, NormFinder and Cq results. Pairwise variation of v values indicated that the use of two RGs is recommended for normalization of RT-qPCR results in each sample group.

Spearman's correlation analysis
To confirm the reliability of the five computational tools employed in the current study, Spearman's correlation analysis was performed on the results of the ranking of RGs for each sample group. A strong correlation was found between the results obtained with comparative Cq value, NormFinder and RefFinder tools. Correlation coefficient be-tween Cq and NormFinder was 0.843; between RefFinder and Cq was 0.942; between NormFinder and RefFinder was 0.797. A moderate correlation was found in the following comparisons: Cq vs. BestKeeper (0.756), RefFinder vs. geNorm (0.747), Cq vs. geNorm (0.732). The lowest correlation coefficient was found between NormFinder and geNorm (0.501) ( Table 6).

RGs validation
To evaluate the effect of normalization of RT-qPCR results based on the optimal RGs, a validation experiment was conducted in SOL and GAS. The expression of genes TNNI1,  TNNC1, TNNI2, TNNC2 was normalized using RGs with high (ACTB, HPRT, YWHAZ) and poor performance (GAPHD and 18S) as demonstrated in the stability ranking ob-tained previously (Fig. 3). The combined use of ACTB, YWHAZ, HPRT yielded a more stable normalization of RT-qPCR results of the four target genes. Considering the normal-ization of TNNI1 and TNNC1, fold changes were higher using GAPDH as the RG, whereas  an opposite trend was found when 18S was used as the RG (Fig. 4AB). No sig-nificant differences were found in the expression of TNNI2 and TNNC2 in SOL and GAS when using GAPDH as RG. A higher mRNA expression of TNNI2 and TNNC2 was observed in GAS when 18S was used alone as the RG (Fig. 4CD). Therefore, the combined use of ACTB, YWHAZ, HPRT as RGs was more suitable for the normalization of RT-qPCR results in murine skeletal muscle tissue than GADPH and 18S.

DISCUSSION
In the present study, five different algorithms were used to compared the expression stability of eight RGs in mouse skeletal muscle tissue. In addition, the effects of genetic background, skeletal muscle type and growth stages in the expression stability of selected RGs was compared. Furthermore, the C2C12 myoblast cell line, which has been widely used as an in vitro model to mimic skeletal muscle tissue, was included in the current study. Firstly, raw Cq values of RGs were evaluated; the higher the Cq values, the lower the abundance of nucleic acids. The distribution of Cq values of eight RGs indicated that 18S was the most abundant transcript in murine skeletal muscle and C2C12 myoblasts (Fig. 1). This is in agreement with another study employing mouse skeletal muscle tissue (Piazza et al., 2017). Cq values of all RGs evaluated herein were below 30, thus indicating that the expression levels of RGs were acceptable (Cedraz de Oliveira et al., 2017). Subsequently, RT-qPCR results were analyzed by Cq method and NormFinder. The results of Cq method and NormFinder were highly comparable, which showed that ACTB, YWHAZ, HPRT performed better (higher ranking) and 18S, PPIA, GAPDH performed poorly (lower ranking). Based on the results obtained with BestKeeper, ACTB, HPRT, YWHAZ and B2M were higher-ranked, whereas GAPDH and 18S were lower-ranked. These results were not in complete agreement with those obtained with Cq and NormFinder, although poorranked RGs were comparable. This could be linked to the differences in the algorithm problem solving method. geNorm provided the ranking of RGs as well as the number of recommended RGs. The present results indicate that the two RGs met the recommended threshold, i.e., V value <0.15 (Vandesompele et al., 2002). When V values are <0.15, it is not necessary to increase the number of RGs for normalization. However, V values >0.15 have also been reported in previous studies, which could be affected by the number of RGs and sample types (Kuijk et al., 2007). The accuracy of different tools influenced the selection of RGs. Thus, it can be hypothesized that there may have been a bias for the selection of RGs when using a single computational tool. Previous studies reported using multiple computational tools in combination of RGs selection (Fan et al., 2020;Jose et al., 2020). (De Spiegelaere et al., 2015) have reported that the result of RefFinder may be biased as it does not account for PCR amplification efficiences. Given this we calculate the geometric mean of the RGs ranking values obtained from above four methods to compare with RefFinder results (Fig. S2). This result is similar to the RefFinder result. The ranking of RGs based on the results of Cq, NormFinder, BestKeeper, geNorm, and comprehensively of RefFinder. Collectively, ACTB, HPRT and YWHAZ were considered the RGs with the best performance, whereas 18S, GAPDH and TUBA performed poorly.
Conversely, slight differences in ranking of RGs were observed between different sample groups. When considered the two different genetic backgrounds evaluated in the study, the most stable RGs were ACTB and GAPDH in ICR mouse, whereas ACTB and TUBA in C57BL/6 as indicated by Cq and NormFinder; in contrast, GAPDH and ACTB were found to be the most stable RGs both in ICR and C57BL/6 by BestKeeper. Using geNorm, ACTB was the most stable RG in the two genetic backgrounds. Kristen et al. (Thomas et al., 2014) reported that the stability ranking of RGs differed in three mouse strains (R129, C57BL/6J and C57BL/10). Thus, our results indicated that differences linked to mice genetic backgrounds were not significant when conducting RT-qPCR experiments on samples of the same tissue, and the stability of good RGs was similar among samples.
Additionally, the effect of muscle tissue type (TA, LD, GAS and SOL) on the stability of RGs was evaluated. In Cq analysis, the difference of stability between good RGs in TA and in GAS was low, HPRT = 0.45, B2M = 0.48, YWHAZ = 0.49, TUBA = 0.49, and HPRT = 0.69, B2M = 0.69, YWHAZ = 0.72, TUBA = 0.73, respectively. Thus, based on Cq and NormFinder results, HPRT and YWHAZ were considered, respectively, the optimal RGs in TA, LD and GAS, whereas ACTB and GAPDH were considered better RGs in SOL. SOL is a slow-twitch muscle fiber type, whereas TA, LD, GAS are fast-twitch muscle fiber types, which may have accounted for the observed discrepancies. Based on BestKeeper and geNorm results, high-ranked RGs were inconsistent in different skeletal muscles. Previous studies have focused on the selection of RGs in single skeletal muscle types (Thomas et al., 2014;Niu et al., 2016;Hildyard, Finch & Wells, 2019). The results discussed herein suggest that there is a difference on RGs stability in various murine skeletal muscle types.
Considering growth stage (mice in lactation and adult mice), HPRT, ACTB and YWHAZ were among the top four high-ranked among the five algorithms evaluated in the present study. Niu et al. (2016) reported candidate RGs in porcine skeletal muscle in 26 different developmental stages. Moreover, it has been suggested that the commonly used RGs GAPDH and ACTB may be not suitable for RT-qPCR experiments in skeletal muscle in different growth stages. Another study described that PPIA and HPRT were the most stable RGs in porcine longissimus dorsi muscle tissue in different developmental stages (Feng et al., 2010). Furthermore, C2C12 myoblast cell line was used herein as an in vitro cell model to study the stability of RGs in skeletal muscle; collectively, PPIAA, ACTB, HPRT, YWHAZ were more stable in C2C12 myoblasts, whereas YWHAZ, TUBA, HPRT, B2M were more stable in skeletal muscle tissue.
Considering all the above findings, ACTB, HPRT, YWHAZ, GAPDH and 18S were selected as RGs to normalize RT-qPCR results of target genes in SOL and GAS, which included TNNI1, TNNC1, TNNI2 and TNNC2. Troponin is the sarcomeric Ca2+ regulator for striated muscle contraction. TNNI1 and TNNC1 are exclusively expressed in slow skeletal muscle fiber types, while TNNI2 and TNNC2 are in fast skeletal muscle fiber types (Gomes, Potter & Szczesna-Cordary, 2002). The normalization of gene expression analysis based on the top ranked RGs (ACTB, HPRT and YWHAZ) yielded more replicable results. When normalizing TNNI1 and TNNC1, a significant difference was observed in SOL and GAS when using GADPH and 18S as RGs, although fold changes were discrepant. When normalizing TNNI2 and TNNC2, no significant differences were found in SOL and GAS using GADPH as RG. John et al. (Hildyard, Finch & Wells, 2019) found that normalization using GADPH yielded conflicting results in a mouse model of Duchenne muscular dystrophy compared with normalization conducted based on other RGs. This suggests that it may be not appropriate to conduct normalization of target genes in RT-qPCR experiments in skeletal muscle tissue using GADPH as the RG.

CONCLUSIONS
In the present study, the expression stability of eight RGs in mouse skeletal muscle tissue was evaluated by five different computational tools. A strong correlation was found among the results obtained with Cq, NormFinder and RefFinder. Thus, a joint analysis of these tools is proposed to enable the proper selection of RGs for similar gene expression experiments. More specifically, based on the data discussed herein, the selection of ACTB as a RG in murine skeletal muscle tissue experiments may be the best choice for more reliable results of gene expression analysis. Moreover, if the experimental conditions permit, the combined use of ACTB, HPRT and YWHAZ could be considered the best option to normalize expression levels of target genes in murine skeletal muscle tissue experiments. The present study provides a useful guide for the selection of RGs for RT-qPCR experiments using murine skeletal muscle tissue.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This research was funded by Sichuan Science and Technology Program, grant number 2020YFN0147; 2021YJ0265; scsztd-2022-08-09, the China Agriculture Research System of MOF and MARA, and the Open Project of Chongqing Municipal Key Laboratory of Pig Science, grant number 21516. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: Sichuan Science and Technology Program: 2020YFN0147, 2021YJ0265, scsztd-2022-08-09. China Agriculture Research System of MOF and MARA. Open Project of Chongqing Municipal Key Laboratory of Pig Science: 21516.