Identification of optimal reference genes for quantitative PCR studies on human mesenchymal stem cells

Introduction

Mesenchymal stem cells (MSCs) were initially identified as an adherent, fibroblast-like population obtained from adult bone marrow (BM) by Friedenstein et al (1,2). MSCs are able to differentiate into multiple lineages (3,4) and are increasingly proposed as a therapeutic strategy for tissue regeneration and repair (5). Although BM has been the primary source of MSCs in the past (6–8), the use of BM-derived MSCs (BM-MSCs) is limited due to multiple factors, including the high degree of viral exposure, potential donor morbidity, as well as significant decreases in cell number and proliferation/differentiation capacity associated with age (9), and the highly invasive procedure required in order to obtain BM. Therefore, it was important to find alternative sources to provide MSCs. Further studies have identified additional MSC sources, including adult synovial membranes and the fetal liver and spleen (10–12). However, there has only been a limited number of studies on MSCs isolated from human fetuses (gestational age, 12–16 weeks). Comparison of fetal tissue-derived MSCs (FT-MSCs) (13) and adult-derived MSCs revealed that the biological activity and the differentiative and multiplication capacity of the former were greater than those of the latter (14). Comparison of the in vitro and in vivo characteristics of BM-MSCs and FT-MSCs requires analyses of their respective gene expression profiles in order to elucidate their fundamental mechanisms, including self-renewal during long-term expansion, differentiation into mature cells and tissue-repair properties.

Quantitative polymerase chain reaction (qPCR) is a commonly used technique to determine the relative change in mRNA expression of target genes. Due to the accuracy, ease of use and reproducibility of qPCR analysis, it is frequently used in MSC research. However, qPCR accuracy is influenced by various external and internal factors, including the amount of starting sample, RNA preparation, cDNA synthesis and PCR efficiency. Therefore, it is necessary to normalize gene expression levels by comparison to reference genes (RGs) as internal controls (15). An ideal RG should not be influenced by cell cycle, cell passages or experimental conditions (16); simultaneously, it should be stably expressed in various samples (17,18). However, to the best of our knowledge, no single RG has been reported to be universal and completely constant. Furthermore, increasing evidence indicated that the expression levels of commonly used RGs vary significantly between cell types and experimental conditions (19,20). Thus, the selection of suitable RGs for idiographic study is a prerequisite for any qPCR assay to obtain reliable results. The aim of the present study was to identify and assess the stabilities and reliabilities of ten RGs which are commonly used in BM- and FT-MSCs for qPCR.

Ten common RGs, including 18S, ACTB, B2M, HPRT, GAPDH, TBP, PPIA, RPLP0, PGK1 and RPL13A (Table I)were selected for the present study and their expression stabilities were analyzed using geNorm (21), NormFinder (22) and BestKeeper (23) software. The present study aimed to identify the optimal RGs for further research on BM-MSCs and FT-MSCs.

Table I Summary of reference genes used in the present study.

Table I

Summary of reference genes used in the present study.

Symbol	Name	Function	Accession number
18S	18S ribosomal RNA	Ribosomal subunit	NM_10098.1
GAPDH	Glyceraldehyde-3-phosphate dehydrogenase	Enzyme in glycolysis and nuclear functions	NM_002046
RPLP0	Ribosomal protein, large, P0	Structural component of the 60S subunit of ribosomes	NM_001002.3
ACTB	Beta-actin	Cytoskeletal structural actin	NM_001101
PPIA	Peptidyl-prolylisomerase A	Accelerates the folding of proteins	NM_021130.3
PGK1	Phosphoglycerate kinase 1	Glycolytic enzyme	NM_000291.3
B2M	Beta-2-microglobulin	Component of the MHCI molecules	NM_004048.2
RPL13A	Ribosomal protein L13a	Structural component of the 60S ribosomal subunit	NM_012423.2
HPRT	Hypoxanthine phosphoribosyl transferase 1	Enzyme in purine metabolic pathway	NM_000194
TBP	TATA box binding protein	General transcription factor	NM_003194

[i] MHCI, major histocompatibility complex class I.

Materials and methods

MSCs

The study was approved by the Ethical Committee of China-Japan Union Hospital, Jilin University (Changchun, China). BM-MSCs were isolated from femur-derived bone marrow samples that were obtained by surgical operation (China-Japan Union Hospital) on otherwise healthy patients (aged between 18 and 43 years) following receipt of their informed consent. FT-MSCs were obtained from Jilin Zhongke Bio-engineering Co., Ltd. (Changchun, China). For all experiments, pools of the various cell types were prepared by mixing equal numbers of cells from five donors of the same passage number. Cells were not cultured for more than four passages.

Identification of MSCs

Flow cytometric characterization of MSCs

MSCs of passage three were labeled with the following anti-human antibodies: CD14-phycoerythrin (PE), CD34-PE, CD45-fluorescein isothiocyanate (FITC), CD73-PE, CD90-FITC, CD105-peridinin chlorophyll, CD44-PE (BD Biosciences, San Jose, CA, USA). A total of 106 labeled cells were evaluated by flow cytometry (Beckman Coulter Fc500, Brea, CA, USA) and the data were analyzed with CXP software (Beckman Coulter Fc500).

MSC differentiation potential

For the differentiation of MSCs into adipocytes and osteoblasts, cells were incubated in adipogenesis differentiation medium (StemPro® Adipogenesis Differentiation kit; Gibco-BRL, Invitrogen Life Technologies, Carlsbad, CA, USA) and osteogenesis differentiation medium (StemPro® Osteogenesis Differentiation kit; Gibco-BRL), respectively according to the manufacturer’s instructions. Adipogenic differentiation was measured by staining cells in wells with Oil Red O (Sigma-Aldrich, St. Louis, MO, USA) on day 21 of culture following fixing the cells with 5% paraformaldehyde (Beyotime, Shanghai, China) for 5 min at room temperature. Following 35 days of incubation, osteogenic differentiation was evaluated by staining the cells with Alizarin Red S (Alizarin S staining kit; Genmed, Shanghai, China).

RNA isolation and cDNA synthesis

Total RNA was isolated using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) according to the manufacturer’s instructions and RNA integrity was electrophoretically verified by ethidium bromide staining. RNA concentrations and A260/A280 nm absorbance ratios were measured spectrophotometrically with a Synergy HT Multi-Mode microplate reader (BioTek Instruments, Inc., Winooski, VT, USA). From 500 ng of RNA, cDNA was synthesized using an RNA PCR kit (avian myeloblastosis virus; TaKaRa, Dalian, China) according to the manufacturer’s instructions. cDNA was stored at −20°C until required.

Selection of candidate RGs and primer design

RGs were selected for analysis based on a literature search on the subject of RG studies in MSCs. In the present study, ten candidate RGs were selected and were as follows: 18S, ACTB (24,25), B2M (26,27), HPRT, GAPDH (28–34), TBP, PPIA, RPLP0, PGK1, RPL13A. The full name, function and accession number of the RGs assessed in the present study are listed in Table I. RGs were selected from varying functional classes, which significantly reduces the chance that the genes may be co-regulated.

Primer pairs used for qPCR were designed using primer three input (http://flypush.imgen.bcm.tmc.edu/primer/primer3_www.cgi). The qPCR primers were synthesized by Sigma Genesys (Sigma-Aldrich) with melting temperature (Tm) at 60±1°C. All primers were purified by Ultrapage (Sangon, Biotech, Shanghai, China). Primer efficiencies were determined using a 10-fold dilution series of cDNA as templates for qPCR reactions. An approximation of PCR efficiency was calculated using the slope of the calibration curve according to the following equation: E=10−1/slop, where ‘slop’ represented the linear regression slope (18). Reactions were performed in triplicate and data were analyzed by using the 2−ΔΔCt method (10). The primer sequences and corresponding amplicon sizes are listed in Table II.

Table II Primer sequences, product sizes and PCR efficiency.

Table II

Primer sequences, product sizes and PCR efficiency.

Gene	Primer sequences (5′-3′)	Product size (bp)	PCR efficiency
18S	F-GTGGAGCGATTTGTCTGGTT R-AACGCCACTTGTCCCTCTAA	115	1.90
GAPDH	F-ATGGGGAAGGTGAAGGTCG R-GGGGTCATTGATGGCAACAATA	108	1.99
RPLP0	F-CTGGAAGTCCAACTACTTCCT R-CATCATGGTGTTCTTGCCCAT	160	2.74
ACTB	F-GAAGATCAAGATCATTGCTCCT R-TACTCCTGCTTGCTGATCCA	111	1.89
PPIA	F-TCCTGGCATCTTGTCCAT R-TGCTGGTCTTGCCATTCCT	179	2.17
PGK1	F-GCCACTTGCTGTGCCAAATG R-CCCAGGAAGGACTTTACCTT	102	2.62
B2M	F-CTATCCAGCGTACTCCAAAG R-GAAAGACCAGTCCTTGCTGA	188	2.08
RPL13A	F-CGAGGTTGGCTGGAAGTACC R-CTTCTCGGCCTGTTTCCGTAG	121	2.00
HPRT	F-CCTGGCGTCGTGATTAGTGAT R-AGACGTTCAGTCCTGTCCATAA	131	1.78
TBP	F-GCACAGGAGCCAAGAGTGA R-GTTGGTGGGTGAGCACAAG	174	2.10

[i] PCR, polymerase chain reaction; F, forward; R, reverse.

qPCR

qPCR was performed in 96-well plates with the ABI PRISM 7500 Sequence Detection system (Perkin Elmer, Inc., Waltham, MA, USA). PCR conditions were as follows: 50°C for two minutes, 95°C for 10 min followed by 40 cycles of denaturation at 95°C for 15 sec, annealing at 58°C for 15 sec and extension at 72°C for 30 sec during which fluorescence was measured. Expression levels were recorded as cycle threshold (Ct). Data were acquired using the 7500 Software (Applied Biosystems Life Technologies, Foster City, CA, USA). The mean Ct values of the triplicate reactions were used for data analysis.

Data analysis

Expression stabilities of the ten RGs were assessed via the three commonly used software programs geNorm, NormFinder and BestKeeper. In geNorm and NormFinder, Ct values were converted into relative quantities via the formula 2−(Ct−lowest Ct). The raw Ct values were used directly for BestKeeper analysis. These three programs are based on Microsoft Excel (Microsoft Corp., Redmond, WA, USA) using various algorithms to determine the expression stability of RGs. GeNorm calculated a gene expression stability measure (M) and pairwise variation (V) parameter. M is the mean pairwise variation for a given gene compared to other tested genes and stepwise exclusion of the gene with the highest M value followed by recalculation was performed until the two most stable genes were left. Lower M values represent higher expression stability. M=1.5 was used as an experimental parameter. Above this value, the gene was considered to be unreliable as an RG. V was calculated to determine the minimal number of RGs required to normalize the expression of genes of interest. V=0.15 was also used as an experimental parameter; below this value, the number of RGs was sufficient for valid normalization.

NormFinder computed RG stability values via an analysis of variance-based model. Lower values indicated higher stabilities. NormFinder was also able to compare inter- and intra-group variations in gene stability.

BestKeeper analyzed RG stability based on the standard deviation (SD) and coefficient of correlation (r) of all RGs. SD values were obtained from the Ct values of each RG, and r values were the correlation coefficient calculated using Pearson’s pair-wise correlation analyses between each RG and the geometric mean of the Ct values. Those genes with an SD >1.0 were considered to be unreliable as a stable RG and the remaining genes are ranked according to their r values.

Results

Isolation and characteristics of MSCs

Human MSCs were obtained from human BM and FT. The adherent cells had fibroblastic morphologies (Fig. 1). The cell-surface antigen profiles of these cells at three passages in culture were analyzed by flow cytometry. These cells were strongly positive for MSC-specific surface markers, including CD44, CD73, CD90 and CD105, but were negative for CD14, CD34 and CD45. These cells also exhibited mesenchymal differentiation potential, as assessed by culturing them in adipogenic and osteogenic medium (data not shown).

Figure 1 Morphologies of MSCs from two human tissue samples. Cells were isolated and cultured in 100-mm Petri dishes. Cells displayed a fibroblast-like shape. Cells displayed were passage three. Phase contrast magnification, ×100. MSC, mesenchymal stem cell; BM, bone marrow; FT, fetal tissue.

Amplification specificity and efficiency of primers

The primer sequences, corresponding amplicon sizes and PCR efficiencies of primers, are listed in Table II. The amplification performance of each primer pair was tested by qPCR. Primer specificities for qPCR were verified by dissociation curve analysis and 2% agarose gel electrophoresis. No primer-dimers or multiple bands/peaks were detected, confirming a single amplified band of a predicted size (Fig. 2). All primer pairs amplified a single PCR product of expected size, produced a slope >-4.0, exhibited correlation coefficients (r2) >0.97 and the corresponding qPCR efficiencies were in the range of 1.8–2.7 obtained from a 10-fold dilution series of the template cDNA. These data indicated that the amplification efficiencies of the primers analyzed reached the standard requirements of conventional qPCR.

Figure 2 Specificity of primers and amplicon lengths. PCR amplification products were analyzed by agarose gel electrophoresis and dissociation curves. (A) PCR products were run through 2% agarose gel. A single band of anticipated size indicated that the PCR product was specific. The difference in band brightness displays the different expression levels of reference genes. (B) PCR products were analyzed by their respective dissociation curves. The y-axis represents fluorescence intensity while the x-axis shows temperature (60–95°C). A single peak represented a specific PCR product. M, size markers; N, negative control; PCR, polymerase chain reaction.

Expression levels of RGs

The expression levels of all ten RGs were analyzed by comparison of Ct values using two biological samples with three technical replicates. Fig. 3 exhibits the mean Ct values for each gene in two samples. All the RGs displayed Ct values in MSCs, ranging from 12.39 for ACTB to 20.16 for TBP. Among these RGs, ACTB (Ct 12.39–14.19) and PPIA (Ct 13.55–13.62) displayed the highest RNA transcription levels. The lowest RNA transcription levels were observed for TBP (Ct 18.37–20.16), followed by 18S (Ct 13.36–16.02). The individual RGs had varying expression ranges across the samples. Among the ten RGs in the present study, 18S had the greatest variation in its transcript expression levels (2.7 cycles), while PPIA and B2M had significantly lower expression variations (0.06 and 0.14 cycles, respectively). These Ct standard deviations indicated an initial hypothesis concerning the variability in expression. The RGs were ranked from greatest to lowest variability in expression: 18S>TBP>ACTB>PGK1>GAPDH>RPLP0>RPL13A>HPRT>B2M>PPIA.

Figure 3 Expression levels of the reference genes in the two cell groups. Values are presented as the mean ± standard deviation (n=5). Ct, cycle threshold; BM-MSC, bone marrow derived mesenchymal stem cell; FT-MSC; fetus tissue derived MSC.

GeNorm analysis

The average stability values (M) of the ten RGs in all samples analyzed are indicated in Fig. 4A. PPIA and B2M (M=0.154) were identified as the most stably expressed RGs, while ACTB (M=1.031) was the least stably expressed. These ten RGs demonstrated inconspicuous variation between BM-MSCs and FT-MSCs as the M-value of each gene was below the default limit of M≤1.5. The stability ranking of the selected RGs was B2M/PPIA>RPL13A>HPRT>RPLP0>GAPDH>PGK1>TBP>18S>ACTB. GeNorm defines a pairwise variation of 0.15 as the cutoff value, below which the inclusion of an additional RG is unnecessary (28). Here, the V2/3 value (the pairwise variation upon increasing the number of normalization factors from two to three) was 0.127, which was below the cutoff value; hence, the use if a combination of two RGs was sufficiently stable for the analysis of BM-MSCs and FT-MSCs (Fig. 4B).

Figure 4 Selection of the most suitable RGs for normalization among samples using geNorm analysis. (A) M of the ten RGs analyzed. The x-axis from left to right indicates the ranking of the genes according to their expression stability, lower M values indicate higher expression stability. (B) Determination of the optimal number of RGs for normalization was conducted. Software calculated the normalization factors from at least two genes at which variable V defines the pair-wise variation between two sequential normalization factors. V2/3 exhibited the value below the cut-off value of 0.15, indicating that use of two RGs for normalization is necessary, whereas the addition of a third RG is optional. M, expression stability measures; RGs, reference genes.

NormFinder analysis

According to the outputs from the NormFinder analysis, the most stable gene was RPL13A (stability value, 0.082), followed by PPIA and B2M. The most unstable genes were ACTB, 18S and TBP (stability values, 0.909, 0.885 and 0.707, respectively; Fig. 5). Among the RGs analyzed in the present study, the ranking of the stabilities according to NormFinder software was RPL13A>PPIA>B2M>RPLP0>GAPDH>HPRT>PGK1>TBP>18S>ACTB. The most stable combination was PPIA and RPLP0, which produced a stability value of 0.131.

Figure 5 Expression stability values of the reference genes calculated by NormFinder. Lower stability values indicated more stably expressed genes. The x-axis from left to right indicates the ranking of the genes according to expression stability.

BestKeeper analysis

The results of BestKeeper analysis are displayed in Table III. According to BestKeeper, genes with an SD>1.0 are considered to have an unacceptable range of variation. The results of the present study indicated that all of the RGs analyzed, excluding 18S, had an SD level under 1.0 and may therefore be used as credible RGs (Fig. 6A). 18S had an SD>1.0 and was therefore eliminated from further analysis and the remaining genes were ranked according to their coefficient of correlation (r). HPRT was the most unstable RG and ACTB was the most stable RG according to the results of BestKeeper analysis (Fig. 6B). The ranking of the stabilities was ACTB>TBP>RPL13A>B2M>RPLP0>GAPDH> PGK1>PPIA>HPRT.

Figure 6 Evaluation of RGs using BestKeeper analysis software. (A) Std dev values of the reference genes evaluated by BestKeeper. (B) r values of RGs evaluated by BestKeeper. Higher stability values implied greater stability. std dev, standard deviation; r, coefficient of correlation; RG, reference gene.

Table III Stability of ten reference genes evaluated by BestKeeper.

Table III

Stability of ten reference genes evaluated by BestKeeper.

	ACTB	PPIA	RPLP0	RPL13A	HPRT	GAPDH	B2M	PGK1	18S	TBP
geo Mean [CP]	13.25789	13.59014	13.76474	13.84548	14.15068	14.15517	14.18604	14.48262	14.6256	19.23467
ar Mean [CP]	13.289	13.5905	13.777	13.84875	14.157	14.1765	14.18625	14.50275	14.688	19.26175
min [CP]	12.279	13.492	13.005	13.541	13.446	13.063	14.109	13.648	13.048	18.14
max [CP]	14.343	13.736	14.452	14.207	14.487	15.228	14.3	15.706	16.127	20.841
std dev [±CP]	0.899	0.091	0.5515	0.29775	0.3555	0.6255	0.06975	0.64675	1.329	0.89475
CV [%CP]	6.764994	0.669585	4.003049	2.150014	2.511125	4.412232	0.491673	4.459499	9.048203	4.645217
coeff. of corr. [r]	−0.951	0.468	0.903	0.939	−0.341	0.863	−0.92	0.718	0.953	0.948

[i] geo Mean [CP], geometric mean of CP; ar Mean [CP], arithmetic mean of CP; min [CP] and max [CP], the extreme values of CP; stddev [±CP], standard deviation of the CP; CV [%CP], coefficient of variance expressed as a percentage on the CP level; coeff. of corr. [r], correlation coefficients of each gene. CP, Cross Point.

Comparison of software analyses

The results obtained from geNorm, NormFinder and BestKeeper were analyzed comprehensively. GeNorm and NormFinder indicated that PPIA, B2M and RPL13A were the most stable RGs, whereas ACTB, 18S and TBP were the least stable RGs. However, ACTB was ranked as the most stable gene based on r values obtained from BestKeeper analysis. A summary of the rankings produced by the three software programs is presented in Table IV.

Table IV Ranking of reference gene stability.

Table IV

Ranking of reference gene stability.

			BestKeeper
Rank	geNorm	NormFinder	r	SD
1	B2M,PPIA	RPL13A	ACTB	B2M
2	*	PPIA	TBP	PPIA
3	RPL13A	B2M	RPL13A	RPL13A
4	HPRT	RPLP0	B2M	HPRT
5	RPLP0	GAPDH	RPLP0	RPLP0
6	GAPDH	HPRT	GAPDH	GAPDH
7	PGK1	PGK1	PGK1	PGK1
8	TBP	TBP	PPIA	TBP
9	18S	18S	HPRT	ACTB
10	ACTB	ACTB		18S

* Rank 2 is blank since B2M and PPIA were the the two most stable genes identified by geNorm analysis which was unable to distinguish which was first and which was second.

{ label (or @symbol) needed for fn[@id='tfn5-mmr-11-02-1304'] } SD, standard deviation; r, coefficient of correlation.

Discussion

RGs are required to provide a scale in qPCR experiments. If the scale varies significantly between samples in an experiment, it loses the function of measuring gene expression in those samples. An ideal RG should be neither influenced nor regulated by experimental conditions or treatments. Increasing evidence indicated that there is no single RG that is able to be used for multiple experiments; however, an increasing number of studies suggested that a group of putative RGs for certain specific experimental setups may be recommended for future studies (35–37).

MSCs, primarily derived from BM, have been examined widely for their capacities in repairing damaged tissues (38–40). The prospective clinical applications of BM-MSCs are varied. Their ability to differentiate into desired mature cell types and additional indirect mechanisms have been recognized to have important roles in the treatment of autoimmune diseases (40). FT-MSCs, which are self-renewing and pluripotent, are a valuable research tool and have potential for use in regenerative medicine. Thus, fetal tissue provides a potential substitute to BM as a source of MSCs (41). In order to make comparisons between the mRNA expression features of two sources of MSCs, it is important to identify a standardized, reproducible set of RGs to normalize qPCR analysis.

To the best of our knowledge, the present study was the first to validate the stability of RGs used for the comparison of BM-MSCs and FT-MSCs. Expression stabilities of ten candidate RGs, including 18S, GAPDH, RPLP0, ACTB, PPIA, PGK1, B2M, RPL13A, HRPT1 and TBP, were estimated using three statistical algorithm software packages: geNorm, NormFinder and BestKeeper. Evaluation with geNorm as well as NormFinder indicated that B2M, PPIA and RPL13A were the most stably expressed RGs. In addition, the two programs rendered ACTB, 18S and TBP the least stably expressed RGs and the rank orders were also coincident. However, the results of the analysis using BestKeeper significantly differed from those obtained using geNorm and NormFinder. BestKeeper indicated that ACTB, TBP and RPL13A were the three most stable RGs, while 18S, HPRT and PPIA were least stable. Of note, ACTB was ranked as most stable by BestKeeper, but least stable by geNorm and NormFinder. Rank orders obtained were expected to differ between the three software programs given their distinct statistical algorithms. For example, unlike geNorm and BestKeeper, NormFinder used a mathematical model based on assessment of intra- and intergroup variations to determine the optimal RGs and therefore showed decreased sensitivity toward co-regulation of RGs. There was no consensus regarding which was the best method. In the present study, the most stable RGs (RPL13A, PPIA and B2M) determined by geNorm and NormFinder were identical, as were the least stable genes (ACTB, 18S and TBP).

ACTB is widely used as an RG to evaluate target gene expression levels in MSCs (42–46). However, in the present study, ACTB was found to have the lowest stability among the selected RGs in BM-MSCs and FT-MSCs. GAPDH was previously considered as the gold standard RG, despite a lack of experimental evidence (37,47–49). However, GAPDH was ranked sixth by geNorm and fifth by NormFinder and BestKeeper in the present study. It was therefore concluded that GAPDH and ACTB were not reliable RGs for the normalization of qPCR data in BM- and FT-MSC research and they are not recommended for use throughout this field of research.

In conclusion, the present study identified RPL13A, PPIA and B2M as suitable genes for the normalization of qPCR data and indicated that 18S, ACTB, and TBP were unsuitable for normalization of mRNA expression levels of BM- and FT-MSCs.

Acknowledgements

The present study was supported by the National Natural Science Foundation of Jilin Province of China (no. 20130101167JC).

References