Identifying suitable reference genes for gene expression analysis in developing skeletal muscle in pigs

Sample collection, RNA extraction and next generation sequencing

All animals were sacrificed at a commercial slaughterhouse according to protocols approved by the Institutional Animal Care and Use Committee at the Institute of Animal Science, Chinese Academy of Agricultural Sciences (Approval number: PJ2011-012-03). Longissimus dorsi (LD) muscle samples were collected from Landrace fetuses on the following days post-coitum (dpc): 33, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, and 105 dpc. LD muscle samples were collected from piglets on the following days after birth (dab): 0, 10, 20, 30, 40, 60, 80, 100, 140, 160, and 180 dab. Three biological replicates were collected at each time point, and totally 78 samples were collected. All samples were immediately frozen in liquid nitrogen and stored at –80 °C until further processing. Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA quantity and quality was determined by the Evolution 60 UV-Visible Spectrophotometer (Thermo Scientific). RNA preparations with an A₂₆₀/A₂₈₀ ratio of 1.8–2.1 and an A₂₆₀/A₂₃₀ ratio > 2.0 were selected for this assay. RNA integrity was determined by analyzing the 28S/18S ribosomal RNA ratio on 1.5% agar gels. Only RNA preparations that resolved with three clear bands on these gels were used for the transcriptome sequencing and qPCR analysis.

Selection of candidate reference genes

For the purpose of identifying potential reference genes during skeletal muscle development, candidate reference genes were selected according to previous studies (Martino et al., 2011; Uddin et al., 2011; Wang et al., 2015; Zhou, Liu & Zhuang, 2014). The top 15 stably reference genes in the transcriptome data of skeletal muscle at different developmental stages based on the coefficient of variation (CV) were chosen for further gene-stability evaluation by qPCR method. Lower CV values represent genes with more stable expression in our transcriptome data.

cDNA synthesis

The cDNA synthesis was performed using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific) according to the manufacturer’s instructions for reverse transcription (RT) PCR. A mixture of 2 µg of total RNA and 1 µL of random primer was incubated at 65 °C for 5 min to dissociate the RNA secondary structure. Next, the following reaction was carried out in a total volume of 20 µL: 12 µL of the first reaction mixture, 4 µL of 5× RT buffer, 2 µL of 10 nM dNTP, 1 µL of RevertAid Reverse Transcriptase (200 U/µl) inhibitor, and 1 µL RiboLockRNase Inhibitor (20 U/µl). The reverse transcription reaction was performed at 25 °C for 5 min, followed by 42 °C for 1 h and 5 min at 70 °C. The cDNA was then diluted 7-fold, and stored at −20 °C until use.

qPCR with SYBR green

Each qPCR reaction was performed in a final reaction volume of 20 µL containing 10 µL of SYBR Green Select Master Mix, 7 µL of sterile water, 0.5 µl of gene-specific primers, and 2 µL of template cDNA. The PCR amplifications were performed on a 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) under the following cycling conditions: 95 °C for 5 min, followed by 40 cycles at 95 °C for 15 s and 60 °C for 45 s. Three independent individuals at each time point were used for temporal and spatial analyses. Each qPCR reaction was performed in triplicate for technical repeats. The mean quantification cycle (Cq) value was used for further analysis. The primer sequences were according to or based on those of previous reports as follows (Table 1): B2M, RHOA, RPL32, DRAP1 RNF7, WSB, and H3F3 (Wang et al., 2015); ACTB, AP1S1, API5, and VAPB (Tramontana et al., 2008); GAPDH and RPS18 (Park et al., 2015); and TBP and PPIA (Martino et al., 2011; Uddin et al., 2011).

Analysis of gene expression stability

Gene expression data were transformed to relative quantities using geNorm and NormFinder. GeNorm provided a measure of gene expression stability (M) (Vandesompele et al., 2002; McCulloch et al., 2012)

${\displaystyle Mj=\sum _{k=1}^{n}Vjk∕n-1}$ where:

Mj = gene stability measure,

Vjk = pairwise variation of gene j relative to gene k,

n = total number of number of examined genes.

Lower M values represent genes with more stable expression across specimens being compared and generated a ranking of the putative reference gene expression levels from the most stable (lowest M-values) to the least stable (highest M-values). GeNorm also generated a pairwise stability measure, which can be used to evaluate the suitable number of reference genes for normalization.

NormFinder provided a stability measure (SV), identified the most stable gene, and calculated the best combination of two reference genes. This program focuses on finding the two genes with the least intra- and inter-group expression variation or the most stable reference gene in intra-group expression variation. Genes with lower stability values show a stably expressed pattern, while the higher stability values share the least stably expressed pattern.

The BestKeeper program was used to compute the geometric mean of each candidate gene’s Cq value, to determine the most stably expressed genes based on correlation coefficient (r) analysis for all pairs of candidate reference genes (≤10 genes), and to calculate the percentage coefficient of variation (CV) and standard deviation (SD) using each candidate gene’s crossing point (CP) value (the quantification cycle value; Cq). In the BestKeeper program, genes with higher r values and lower CV and SD values are more stable reference genes.

Expression analysis of candidate reference genes in developing skeletal muscle

We performed qPCR assays to measure the expression levels of 15 candidate reference genes (ACTB, API5, B2M, GAPDH, RNF7, H3F3, PPIA, AP1S1, PPIA, RHOA, RPS18, RPL32, TBP, WSB, and VAPB) in the LD muscle samples at 15 embryonic stages (33, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, and 105 dpc) and 11 postnatal stages (0, 10, 20, 30, 40, 60, 80, 100, 140, 160, and 180 dab) of Landrace pigs. To minimize experimental error, triplicate amplifications were performed for individual experiments. The Cq values were computed to quantify the candidate reference gene expression levels. A higher Cq value means lower gene expression levels. Analysis of gene expression stability was based on the Cq values generated by qPCR (Fig. 1). Among all the tested genes, GAPDH had the lowest mean Cq value (15.94) and AP1S1 had the highest mean Cq value (26.36). All candidate reference genes were abundantly expressed in skeletal muscle and showed wide variations in expression levels at different developmental stages. Therefore, it was necessary to evaluate gene expression stability and determine the suitable number of reference genes for accurate gene expression profiling in developing skeletal muscle.

GeNorm analysis of candidate reference gene expression stability

We calculated the gene expression stability values (M value) for the 15 candidate reference genes using the geNorm program. Genes with lower M values have more consistent expression levels. The M values of the candidate genes are presented in Fig. 1. When all developmental stages were analyzed as one data set. The results revealed that API5 and H3F3 had the lowest M values, whereas GAPDH had the highest M value. This indicated that API5 and H3F3 were the most stably expressed gene pair of the 15 candidate reference genes, whereas GAPDH had the most variable expression (Fig. 2) in developing skeletal muscle across prenatal and postnatal periods. In the prenatal muscle samples, API5 and RPS18 expression was the most stable, whereas GAPDH and DRAP1 expression was the most variable (Fig. 3). When only postnatal muscle samples were analyzed, API5 and RPS18 expression was the most stable, whereas B2M expression was the most variable (Fig. 4). The geNorm analysis demonstrated that GAPDH was the most variably expressed gene in all developmental periods, suggesting that GAPDH was not a suitable reference gene for gene expression analysis in developing skeletal muscle. By contrast, the stability of API5, RPS18, and H3F3 expression suggested that they were suitable reference genes to use as internal controls. When gene expression was analyzed in developing skeletal muscle across all tested prenatal and postnatal periods, API5 expression was the most suitable to use as a reference gene for normalization analysis in expression profiling studies.

One single reference gene might not provide sufficient control for gene expression analyses in developing skeletal muscle. Therefore, we used geNorm to analyze the optimal number of reference genes required to obtain reliable results from RT-qPCR studies. GeNorm was used to calculate the average pairwise variation (V) value between two sequential normalization factors; it has a cut-off value of 0.15 for the pairwise variation according to the previous study (Wang et al., 2015), below which the inclusion of an additional reference gene is not required for reliable normalization of qPCR analyses. When all developmental stages were analyzed together, the V_n/V_n+1 value ranged from 0.059 to 0.111, which were all lower than the cut-off value of 0.15 (Fig. 5). These results indicated that two reference genes were optimal for gene expression analysis of all tested developmental periods. The results were similar for gene expression analysis of the embryonic data set (Fig. 6), and two reference genes were sufficient for analysis. By contrast, the V value decreased significantly with the addition of reference genes in the postnatal data set, although all values were lower than 0.15 (Fig. 7). These results indicated that the three most stable reference genes were required for accurate normalization of gene expression data for the postnatal period.

NormFinder analysis of candidate reference gene expression stability

Next, we used NormFinder to rank the most stable and the least stable genes by calculating the gene expression stability value and standard error. NormFinder analyses showed that API5 was the most stable reference gene with the lowest stability value (SV = 0.088) in all tested developmental periods (Table 2). API5 was the most stable reference gene in postnatal periods, whereas RPS18 was the most stable reference gene in prenatal periods (Table 2).

BestKeeper analysis of candidate reference gene expression stability

Then, BestKeeper program was used to evaluate the reference gene expression stability. We used BestKeeper to identify the optimal reference genes on the basis of the correlation coefficient (r), CV, and SD values (Table 3). The program can calculate r values for up to 10 genes. Therefore, we selected the top 10 candidate genes based on the previous results. In the BestKeeper program, genes with higher r values (≥0.900) and lower CV and SD values are considered as stable and suitable reference genes. In all tested developmental periods, API5 expression had the lowest CV value (2.09) and almost the lowest SD value (it was slightly larger than that of RPL32 expression). Therefore, we propose that API5 is the most suitable reference gene for expression analysis of developing skeletal muscle during the tested prenatal and postnatal stages. API5 also was selected as the most stable gene during the postnatal period, whereas VAPB was the most stable gene for the analysis of developing skeletal muscle during the embryonic period.