Optimization of Reference Genes for Normalization of Reverse Transcription Quantitative Real-Time Polymerase Chain Reaction Results in Senescence Study of Mesenchymal Stem Cells

Abstract

Recently, it has been suggested that cellular senescence is associated with stem cell exhaustion, which reduces the regenerative potential of tissues and contributes to aging and age-related diseases. Mesenchymal stem cells (MSCs) attract a large amount of attention in stem cell research and regeneration medicine because they possess multiple advantages and senescent MSCs could be one of the most useful stem cell models in aging studies. It is important to quantitatively evaluate senescence markers to both identify and study the mechanisms involved in MSC senescence. Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) is currently the most widely used tool to quantify the mRNA levels of markers. However, no report has demonstrated the optimal reference genes that should be used to normalize RT-qPCR in senescence studies of MSCs. In this study, we compared 16 commonly used reference genes (GAPDH, ACTB, RPL13A, TBP, B2M, GUSB, RPLPO, YWHAZ, RPS18, EEF1A1, ATP5F1, HPRT1, PGK1, TFRC, UBC, and PPIA) in proliferating or replicative-senescent human adipose-derived MSCs (hAD-MSCs) that were isolated from seven healthy donors aged 29–59 years old. Three algorithms (geNorm, NormFinder, and BestKeeper) were used to determine the most optimal reference gene. The results showed that PPIA exhibited the most stable expression during senescence, while the widely used ACTB exhibited the lowest stability. We also confirmed that different reference genes lead to different evaluations of senescence markers. Our work ensures that results obtained from senescence studies of hAD-MSCs will be appropriately evaluated in both basic research and clinical trials.

Introduction

Recently, cellular senescence has been found to be associated with adult stem cell exhaustion, which causes a decline in the regenerative potential that further leads to aging phenotypes and pathologies in addition to age-related disease [1 –5].

In stem cell research and regeneration medicine, as a type of representative adult stem cell, mesenchymal stem cells (MSCs) have attracted a great deal of attention because of their multiple advantages [6,7]. For example, our and other groups have demonstrated that MSCs can be directed to differentiate into adipocytes, osteoblasts, chondrocytes, neurons, astrocytes, endothelial cells, hepatocytes, and pancreatic cells when cultured under the appropriate conditions [8 –14]. MSCs secrete multiple bioactive molecules that play important roles in maintaining tissue homeostasis and immunomodulation [15,16]. In addition, MSCs appear to be “resident” stem cells in many tissues, where they can be stimulated to proliferate and differentiate when tissue repairs are required [8]. Importantly, MSCs can easily be obtained and cultured in vitro, and they can therefore be used in clinical trials without the limitation of ethics [6,12,17].

Because of these advantages, MSCs are increasingly recognized as a useful stem cell model for aging studies. Zhang et al. reported that ESC-WRT^−/− cells derived MSCs recapitulated the features of premature cellular aging, which included a global loss of H3K9me3 and changes in heterochromatin architecture [18]. Cieslik et al. reported a causal link between deficiencies in MSC differentiation and elevated circulating insulin levels in aging mice and suggested that this might explain the presence of interstitial fibrosis in aging hearts [19].

To advance the study of MSC senescence, it is important to quantitatively evaluate senescence markers to both identify and study the mechanisms of MSC senescence. Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) analysis is currently the most widely used tool for measuring the abundance of gene transcripts. RT-qPCR has many advantages, including sensitivity, speed, and high throughput [20]. As a relatively quantitative method, accurate normalization in this technique is crucial to ensuring reliable evaluations of gene expression [21]. Therefore, reference genes that can be used as normalizers are needed to show no or only minimum variation in the expression levels of genes between different samples and experimental conditions. GAPDH, ACTs (actins), PPIs (cyclophilins), and microglobulin-encoding genes are the most frequently used “traditional” reference genes.

However, the expression levels of many of these genes have been reported to vary considerably between different cells and during different processes [22]. In MSCs, these “traditional” reference genes have also been reported to be unstable, and reference genes viewed as optimal have been shown to vary under different inducing conditions [23]. However, no report has explored which reference genes are optimal for RT-qPCR in senescence studies that use MSCs.

Therefore, in this study, we evaluated the expression levels of 16 commonly used reference genes in proliferating and replicative-senescent MSCs that were derived from human adipose tissues (hAD-MSCs), which were themselves isolated from seven healthy donors aged between 29 and 59 years old. Three algorithms (geNorm, NormFinder, and BestKeeper) were used to determine which the optimal reference gene was.

Materials and Methods

Isolation and expansion of MSCs from adult human adipose tissue

Human adipose tissue was obtained from donors undergoing liposuction according to procedures approved by the Ethics Committee at the Chinese Academy of Medical Sciences and Peking Union Medical College. The isolation and culture procedure was described in our previous study [13,14].

Fluorescence-activated cell sorting analyses

For fluorescence-activated cell sorting (FACS) analysis, the third passage of human adipose-derived MSCs (hAD-MSCs) was used. The cells were detached and washed with PBS containing 0.5% BSA (Sigma) and incubated with primary antibodies for 30 min at 4°C. Working concentrations for primary antibodies against human CD73, CD105, CD90, CD31, CD34, CD44, CD29, and HLA-DR (BD Biosciences) were 10–20 ng/mL, with same-species, same-isotype irrelevant antibody as a negative control. After washing, the cells were incubated with FITC and phycoerythrin (PE)-conjugated secondary antibodies for 30 min at 4°C. Then, washed three times, cells were resuspended in PBS and analyzed by FACSCalibur flow cytometer and the CellQuest Pro software (BD Biosciences).

Senescence-associated β-gal staining

Senescence-associated β-gal (SABG) staining was performed with Senescence β-Galactosidase Staining Kit (Beyotime) according to the manufacturer's instructions. Staining was examined by microscopy.

Western blotting

Briefly, cells were harvested in RIPA lysis buffer (Beyotime). Whole cell protein extracts were qualified by BCA assay (Beyotime), and equal quantities of extracts (20 μg/well total protein) were separated onto NuPAGE 12% polyacrylamide gels (Invitrogen), transferred to nitrocellulose membranes (Bio-Rad), blocked in 5% nonfat milk in Tris-Buffered Saline and Tween 20, and incubated overnight at 4°C with antibodies for P16INK4a (1:500; Abcam), P21^CLP1 (1:800; Abcam), and CyclinA2 (1:150; Abcam), respectively. GAPDH (1:1000; Santa Cruz Biotechnology) on the same membrane was used as a loading control. Signals were revealed after incubation with anti-rabbit IgG secondary antibody (1:3000; NeoBiotechnology) or anti-mouse IgG secondary antibody (1:2000; NeoBiotechnology) coupled to peroxidase by using electrochemiluminescence (Amersham Biosciences). All antibodies were diluted with 3% BSA.

In the analysis of references, GAPDH, ACTB (1:3000; Abcam), and PPIA (1:1000; Cell signaling technology) are on the same membrane. The software Image-Pro Plus 6.0 and Gel-Pro Analyzer 4 were used for gray scale analysis.

Immunofluorescence

The cultured cells were washed with PBS three times, 4% PFA fixed (10 min), and 0.5% Triton X-100 (10 min) to penetrate the cells, then blocked by incubating in 1% BSA for 1 h. The cells were then incubated with antibodies for P16^INK4a (1:100 dilution in 1% BSA; Abcam) and P21^CLP1 (1:100 dilution in 1% BSA; Abcam) overnight at 4°C. The secondary antibody (red) eFluor 660 (eBioscience) was applied at a 1/200 dilution in 1% BSA for 1 h. DAPI was used to stain the cell nuclei (blue) at a concentration of 1.43 μM (1:1000). Then, viewed by OLYMPUS single-laser confocal microscopy to observe P16^INK4a and by fluorescence microscopy to observe P21^CLP1 expression.

RNA isolation and cDNA synthesis

Total RNA was extracted using TRIzol total RNA isolation reagent (Invitrogen) and digested with DNase I and RNase Inhibitor (TAKARA) according to the manufacturer's instructions. Total RNA concentration was verified using a NanoDrop 2000 spectrophotometer (Thermo Scientific). RNA integrities were assessed using electrophoretic techniques. First-strand cDNAs were synthesized using the Quantiscript RT Kit (TaKaRa Bio) according to the manufacturer's instructions. Final cDNA was diluted 1:10 before use in RT-qPCR.

Primer design

All reference candidates' primers for RT-qPCR analysis were designed using the NCBI Primer Designing Tool (www.ncbi.nlm.nih.gov/tools/primer-blast/), selecting only the primers spanning an exon–exon junction, and melting temperature (Tm value) is between 58°C and 62°C (the optimal is 60), PCR product size ranges from 100 to 250 bp. In addition, all designed primer pairs were checked for nonspecific amplification by in-silico PCR (UCSC, http://genome.ucsc.edu) and by performing a BLAST search (NCBI, http://blast.ncbi.nlm.nih.gov). All primer sequences in our study are given in Table S1 (Supplementary Data are available online at www.liebertpub.com/scd) and were synthesized by Sangon.

Primer validation

All primers were examined for their target specificity by endpoint PCR, with hAD-MSCs cDNAs as templates, using the following conditions: the final 25 μL PCR contains 12.5 μL TaqPCR MasterMix (TIANGEN), 1 μL of AD-MSCs cDNA (20 ng), 1 μL of 10 mM of forward and reverse primers, and 10.5 μL of ddH₂O. The PCR program consists of initial denaturation at 94°C for 3 min, 33 cycles of 94°C for 30 s, 60°C for 30 s, and 72°C for 50 s, followed by final extension at 72°C for 5 min. Amplification products were checked on 2% agarose gels for a single band of correct size.

The linearity of target amplification for calculation of primer efficiency was evaluated using quintuplicate serial dilutions (1:1, 1:10, 1:100, 1:1000, and 1:10000), thus, cDNA sample (100–0.01 ng/μL) as templates on a Bio-Rad CFX96 Real-Time System instrument (Bio-Rad Laboratories) as described below. The quantification cycle values versus cDNA concentration input were plotted to determine the slope values and correlation coefficients (R ²). The corresponding RT-qPCR primer efficiencies (E) were calculated according to the following equation: E = 10^−1/slope.

Quantitative real-time PCR

RT-qPCR was performed in a Bio-Rad CFX96 Real-Time PCR Detection System instrument (Bio-Rad Laboratories) with SYBR Green I real-time PCR kit (TaKaRa Bio). Each reaction was performed in triplicates in a reaction volume of 20 μL in 96-well microtiter plates (Applied Biosystems). All reactions contained 1 μL of cDNA (20 ng), 10 μL of SYBR Green Mix, 1 μL of 10 mM of each primer, and 8 μL of DEPC-treated water. The reaction protocol starts with 3-min template denaturation step at 95°C, 40 cycles of 95°C for 5 s, and 60°C for 30 s, and a melt curve at the end of the cycling protocol with measurement from 55°C to 95°C. Blank controls run for each primer master mix.

Baseline and threshold values of RT-qPCR were automatically determined by the Bio-Rad CFX Manager software, and results were imported to Microsoft Excel. The mean values of triplicates were obtained for further calculations. Raw Cq values were transformed to quantities based on 2^−▵▵Cq method. The data obtained have been converted to proper input files according to requirements of the three programs, then analyzed using geNorm (version 3.5), NormFinder (version 0.953), and BestKeeper VBA applets.

Statistical analyses

All statistical analyses were performed with SPSS software version 17.0. The P-values were calculated using a Student t-test, and P-values <0.05 were considered to be statistically significant.

Results

Senescence model of hAD-MSC

The hAD-MSCs were isolated from seven healthy donors, whose ages and genders were 29 (female), 29 (male), 30 (female), 40 (female), 41 (female), 47 (female), and 59 (female). According to these origins, these cells were named A29, A29m, A30, A40, A41, A47, and A59, respectively. Immunophenotype profile showed that they expressed the typical MSC surface antigens and were negative for markers (Supplementary Fig. S1a). The third passage of each hAD-MSC was considered as the proliferating stage and marked with A29P, A29mP, A30P, A40P, A41P, A47P, and A59P, respectively. Cells in the senescent stage were obtained by culturing them through prolonged passages until the doubling time was >1 week. These cells were correspondingly labeled A29S, A29mS, A30S, A40S, A41S, A47S, and A59S.

In the senescent stage, the cells underwent irreversible growth arrest and displayed a combination of senescence-related phenotypes, including the morphological changes of cells and nuclei (Supplementary Fig. S1b), the increased secretion of SABG (Supplementary Fig. S1c), the increased expression levels of the senescence-related proteins p21^CLP1 (CDKN1A) and p16^INK4a (Supplementary Fig. S1d–f), and meanwhile the reduced expression levels of the proliferation-related protein cyclinA2 (Supplementary Fig. S1d). MSCs derived from different donors became senescent during different passages displaying individual differences. These results show that our cellular senescence model of hAD-MSCs was successfully achieved.

RNA quality control

To ensure that our results are accurate and reliable, all 14 RNA samples were digested using DNAse I to eliminate contamination by genomic DNA. The A260/A280 ratios of the RNA samples were between 1.82 and 2.08, and the mean A260/A230 ratio was 2.04. We also performed agarose gel electrophoresis to determine the integrity of four randomly selected total RNA samples, and the samples demonstrated a 28S/18S rRNA ratio that was ∼2 (Supplementary Fig. S2d). Furthermore, we tested all of the cDNA samples for residual genomic DNA using PCR with a special pair of primers that was located in the neighboring exons 5 and 6 of the HNRNPA2B1 gene. The results showed that all of our obtained RNA samples were free of genomic DNA contamination because no genomic PCR amplicon was obtained (Supplementary Fig. S2a). Thus, all of the RNA samples were suitable for RT-qPCR analysis.

Reference candidate and primer quality control

We chose 16 commonly used references genes (GAPDH, ACTB, RPL13A, TBP, B2M, GUSB, RPLPO, YWHAZ, RPS18, EEF1A1, ATP5F1, HPRT1, PGK1, TFRC, UBC, and PPIA) as “reference candidates.” The details related to these candidates and their primers are provided in Table 1 and Supplementary Table S1.

Table 1.

Candidate Reference Genes Analyzed in this Study

Symbol	Name	Accession number	Function
GAPDH	Glyceraldehyde-3-phosphate dehydrogenase	NM_001256799	Dehydrogenase, oxidoreductase, catalyzes the reversible oxidative phosphorylation of glyceraldehyde-3-phosphate.
RPS18	Ribosomal protein S18	NM_022551	A ribosomal protein that is a component of the 40S subunit and belongs to the S13P family of ribosomal proteins.
ACTB	Actin, beta	NM_001101	Cytoskeletal actins, involved in cell motility, structure, and integrity.
UBC	Ubiquitin C	NM_021009	A polyubiquitin precursor.
TFRC	Transferrin receptor	NM_001128148	Transferrin receptor activity regulates cellular iron ion homeostasis.
RPLP0	Ribosomal protein, large, P0	NM_001002	A ribosomal protein that is a component of the 60S subunit and belongs to the L10P family of ribosomal proteins.
GUSB	Glucuronidase, beta	NM_000181	A hydrolase that degrades glycosaminoglycans, including heparan sulfate, dermatan sulfate, and chondroitin-4,6 sulfate.
PPIA	Peptidylprolyl isomerase A (cyclophilin A)	NM_001300981	The peptidyl-prolyl cis-trans isomerase (PPIase) family catalyzes the cis–trans isomerization of proline imidic peptide bonds in oligopeptides and accelerates the folding of proteins.
RPL13A	Ribosomal protein L13a	NM_001270491	A member of the L13P family of ribosomal proteins that is a component of the 60S subunit.
PGK1	Phosphoglycerate kinase 1	NM_000291	A glycolytic enzyme that catalyzes the conversion of 1,3-diphosphoglycerate to 3-phosphoglycerate.
TBP	TATA box binding protein	NM_001172085	Modulates the DNA binding activity, and modulation of DNA binding affects the rate of transcription complex formation and initiation of transcription.
B2M	Beta-2-microglobulin	NM_004048	A serum protein, in association with the major histocompatibility complex class I heavy chain.
ATP5F1	ATP synthase, H+ transporting, mitochondrial Fo complex, subunit B1	NM_001688	A subunit of mitochondrial ATP synthase.
HPRT1	Hypoxanthine phosphoribosyl transferase 1	NM_000194	A transferase catalyzes conversion of hypoxanthine to inosine monophosphate and guanine to guanosine monophosphate.
EEF1A1	Eukaryotic translation elongation factor 1 alpha 1	NM_001402	An isoform of the alpha subunit of the elongation factor-1 complex, which is responsible for the enzymatic delivery of aminoacyl tRNAs to the ribosome.
YWHAZ	Tyrosine monooxygenase/tryptoph-an5 monooxygenase activation protein, zeta	NM_001135699	A member of 14-3-3 family of proteins mediates signal transduction by binding to phosphoserine-containing proteins.

All of the primers for the candidate genes were chosen according to the general rules of qPCR primer design. First, each pair of primers was confirmed to identify only one specific PCR amplicon with a known expected length (Supplementary Fig. S2c). Second, the primer efficiency E values ranged between 0.91 and 1.09 (Supplementary Table S1), all of the primers showed correlation coefficients of R ² > 0.99. Third, a melting curve analysis showed a single peak for all of the primer pairs, which further confirmed their specificity (Supplementary Fig. S2b).

Expression levels of reference candidates

We analyzed all of the obtained Cq values for each reference candidate to evaluate their relative expression levels in the 14 hAD-MSCs samples. Figure 1a shows their expression level, ranging from the lowest Cq mean (14.6 for EEF1A1) to the highest Cq mean (26.1 for UBC) and the standard deviations (SDs) of the Cq values for the candidates, which ranged between 0.44 (RPS18) and 1.08 (ACTB). Thus, the mean Cq for every reference candidate was within the detection range of RT-qPCR, and the obtained SD values also indicated that ACTB might be more variable, whereas RPS18, PPIA, and RPLPO might be more stable than the other candidates. However, the SD values were also influenced by RNA concentration, the efficiency of reverse transcription, and other factors. Therefore, we evaluated the expression stability of the reference candidates using three approved algorithm software packages (geNorm, NormFinder, and BestKeeper).

FIG. 1.

An evaluation of the stability of reference genes. (a) Real-time quantitative polymerase chain reaction quantification cycle (Cq) values are shown for the 16 candidate reference genes in 14 hAD-MSC samples. The median Cq values are shown as lines, with the 25th to 75th Cq percentiles shown as boxes and the range of Cq values obtained from our analysis of the 14 cDNA samples shown as whiskers (n = 3). (b, c) Results of the evaluation of stability using geNorm. (b) Average expression stability values (M) for the 16 reference candidates after stepwise exclusion of the least stable reference candidate. A lower M value indicates more stable expression. (c) A pairwise variation analysis showing the optimal number of reference candidates to use for normalization. When a pairwise variation value is below 0.15, the least number of reference candidates is considered optimal. (d) Evaluation of stability using NormFinder. The lower the variance of the stability value (M value), the more stable the expression of the gene. The data are representatives of three independent experiments. hAD-MSC, human adipose-derived mesenchymal stem cells.

Evaluation of stability using geNorm

GeNorm is a Microsoft Excel VBA applet that was developed by Vandesompele et al. It uses the transformed data represented by the Cq value to determine the average pairwise variation of a chosen candidate gene in comparison to all other reference candidates and then ranks them based on their average expression stability value (M) from most stable to least stable [21,24]. The lowest M value indicates the most stable gene expression, and a value of 1.5 is recommended as the cutoff when selecting RT-qPCR reference genes. In our analysis, when the results from all 14 hAD-MSC samples were combined, the 16 reference genes were ranked as follows: PPIA/RPLPO > ATP5F1 > RPS18 > YWHAZ > HPRT1 > RPL13A > EEF1A1 > GAPDH > UBC > B2M > PGK1 > TBP > TFRC > GUSB > ACTB (Fig. 1b). Therefore, PPIA and RPLPO showed the highest stability, GAPDH was the ninth best candidate, and ACTB, one of the most widely used references genes, appeared to be the most variable among all of the candidates.

In addition, geNorm also determines the optimum number of reference genes that should be used by calculating the pairwise variation (V_n/V_n+1) between the two sequential normalization factors NF_n and NF_n+1. As a general rule, the stepwise inclusion of reference genes is performed until V_n/V_n+1 drops below a theoretical threshold of 0.15, when adding an extra gene (n + 1) becomes unnecessary for accuracy normalization [24,25]. In our study, the pairwise variations for the total analyses were all below 0.15 (Fig. 1c). There is therefore no need to add a third gene for normalization.

Evaluation of stability using NormFinder

Another commonly used gene expression stability analysis software is NormFinder. In contrast to geNorm, NormFinder evaluates the expression stability of each single candidate independently from the stability of the other genes by taking into account variations across subgroups [25,26]. This approach calculates the inter- and intragroup variation of candidate reference genes and then combines all of the results into a stability value, M. Similarly, the lowest value suggests the highest stability [26]. The results of the NormFinder analysis are shown in Fig. 1d. Consistent with the geNorm conclusions, PPIA (0.008) still ranked at the top position and was followed by EEF1A1, ATP5F1, RPLPO, UBC, YWHAZ, GAPDH, RPS18, RPL13A, HPRT1, TFRC, B2M, TBP, PGK1, GUSB, and ACTB. Hence, GAPDH (0.011) ranked seventh, and ACTB (0.026) was still the least stable reference candidate.

Evaluation of stability using BestKeeper

BestKeeper, like geNorm, is a Microsoft Excel-based stability analysis tool. It uses the raw Ct values of the reference genes to calculate the pairwise correlations between each candidate and the geometric mean of the best-suited pair to determine the optimal normalizer among 10 candidates [27]. Gene expression variation is determined using SD, the percent covariance (CV), and the correlation coefficient (r) for candidate reference genes. The gene with the lowest CV ± SD and the highest r is defined as the most stable reference gene.

To select 10 best reference candidates for final BestKeeper analysis, we have evaluated the candidates in different combinations and compared the resulting data. In our study, six of the reference candidates (ACTB, YWHAZ, EEF1A1, HPRT1, PGK1, and TFRC) exhibiting higher SD and CV values were excluded from further analysis, and the summarized results of the final top 10 reference candidates are listed in Table 2. PPIA and RPLPO were the top ranked genes, with small CV ± SD values of 2.562 ± 0.442 and 2.742 ± 0.431, respectively, and the highest r values 0.958 and 0.961, respectively. They were therefore considered to be the most stable reference genes of the candidates. GAPDH was in the middle, while ACTB showed the highest CV ± SD 4.988 ± 0.841 (data not shown), indicating that ACTB was the least stable reference candidate.

Table 2.

Evaluation of Stability Using BestKeeper

Factors	PPIA	RPLPO	RPS18	ATP5F1	RPL13A	TBP	GAPDH	UBC	GUSB	B2M
N	14	14	14	14	14	14	14	14	14	14
Geo Mean [CP]	17.231395	15.716237	15.589816	19.24605	17.9045	24.698826	15.99583	26.260784	24.341644	17.790162
Ar Mean [CP]	17.238469	15.723969	15.595327	19.251817	17.911267	24.707979	16.004414	26.265792	24.355819	17.795336
Min [CP]	16.578956	15.121046	14.952463	18.513637	17.110446	23.915672	15.412757	25.512725	23.226298	17.105979
Max [CP]	17.981466	16.57235	16.323404	20.098326	18.902116	26.075814	17.173652	27.041183	26.284951	18.592362
Std dev [± CP]	0.4416349	0.4311601	0.3816141	0.3989462	0.4384454	0.5825771	0.4348004	0.4748596	0.6852781	0.3611073
CV [% CP]	2.561915	2.7420564	2.4469769	2.0722519	2.4478751	2.3578502	2.716753	1.8079012	2.8136113	2.0292245
Min [x-fold]	−1.525306	−1.488079	−1.530575	−1.721859	−1.727052	−1.693876	−1.528388	−1.641695	−2.154391	−1.579199
Max [x-fold]	1.6247822	1.7713418	1.6321712	1.8820009	1.9867377	2.5259978	2.3558989	1.6772588	3.8085806	1.7086997
Std dev [x-fold]	1.3551393	1.3454065	1.3003089	1.3159103	1.3521684	1.4931553	1.348781	1.3864788	1.6024979	1.2820886
Coeff. of corr. [r]	0.958	0.961	0.927	0.9	0.878	0.833	0.83	0.802	0.757	0.694
P-value	0.001	0.001	0.001	0.001	0.001	0.001	0.001	0.001	0.0010791	0.0040671
Power of HKG [x-fold]	1.9660545	2.0143285	1.7637719	2.0033571	1.9293219	2.3055655	2.020305	1.8248939	2.6113805	1.552447

N, number of samples; Geo Mean [CP], the geometric mean of Cq; Ar Mean [CP], the arithmetic mean of Cq; Min [CP] and Max [CP], the extreme values of Cq; std dev [± CP], standard deviation of Cq; CV [%CP], the coefficient variance expressed as percentage on Cq level; Min [x-fold] and Max [x-fold], the extreme values of expression levels represented as absolute x-fold over or under-regulation coefficient; std dev [±x-fold], standard deviation of the absolute regulation coefficients. Primary estimation of expression stability is obtained based on CV and standard deviation, then the correlation analysis was combinated to determine the best HKG. Cq, quantification cycle.

Comprehensive ranking of candidate reference genes

We compared the ranking results of the above three analyses (geNorm, NormFinder, and BestKeeper) to come to a conclusion regarding the best reference genes for normalization of RT-qPCR (Table 3). While some of the candidates were ranked in different positions because of differences in the methods used between the different algorithms, PPIA was consistently ranked at the top, GAPDH was consistently in the middle (around seven to ninth), and ACTB was always the last. We then calculated the geometric means of the ranking numbers that were produced by the three methods to determine the final ranking. Those with the smallest geometric mean were considered to be the most stable, while bigger represented less stability. As shown in Table 3, PPIA was ranked at the top, GAPDH was ranked in the middle (at seventh), and ACTB was definitively ranked last.

Table 3.

Sixteen Reference Candidates Ranked by Different Methods

Ranking	geNorm	NormFinder	BestKeeper	Final ranking
1	RPLPO& PPIA	PPIA	PPIA	PPIA
2		EEF1A1	RPLPO	RPLPO
3	ATP5F1	ATP5F1	RPS18	ATP5F1
4	RPS18	RPLPO	ATP5F1	RPS18
5	YWHAZ	UBC	RPL13A	RPL13A
6	HRRT1	YWHAZ	TBP	UBC
7	RPL13A	GAPDH	GAPDH	GAPDH
8	EEFIAI	RPS18	UBC	YWHAZ
9	GAPDH	RPL13A	GUSB	EEF1A1
10	UBC	HPRT1	B2M	HRRT1
11	B2M	TFRC	HRRT1	TBP
12	PGK1	B2M	PKG1	B2M
13	TBP	TBP	YWHAZ	PGK1
14	TFRC	PGK1	EEFIAI	GUSB
15	GUSB	GUSB	TFRC	TFRC
16	ACTB	ACTB	ACTB	ACTB

Relative variation in GAPDH and ACTB when normalized using PPIA

To confirm that PPIA is the most stable reference gene, as indicated by the results of the three algorithms, we attempted to evaluate the relative expression levels of two commonly used reference genes, GAPDH and ACTB, using PPIA to normalize their expression in a RT-qPCR analysis. We then calculated the expression levels of GAPDH and ACTB in all of the senescent hAD-MSCs relative to their values in the corresponding proliferating populations of cells. As a result, the relative expression levels of GAPDH were 1.46, 1.08, 1.32, 1.29, 0.77, 1.10, and 0.78 in A29S, A29mS, A30S, A40S, A41S, A47S, and A59S, respectively (Fig. 2a), and the corresponding expression levels of ACTB were 0.77, 0.45, 1.42, 1.04, 0.54, 1.58, and 1.21, respectively (Fig. 2b). Hence, both GAPDH and ACTB showed variation when normalized using PPIA, and the expression of ACTB was clearly more variable than GAPDH, a finding that was consistent with the analysis of the three algorithms.

FIG. 2.

Relative variation in GAPDH and ACTB expression. The levels of GAPDH (a) and ACTB (b) in all seven senescent hAD-MSCs (shown in black) are calculated relative to the values that were obtained from the corresponding proliferating cells (shown in white). All expression levels were normalized to PPIA. All data are shown as the mean ± SD (n = 3). hAD-MSCs, human adipose-derived mesenchymal stem cells; SD, standard deviation.

The protein expression levels of PPIA, GAPDH, and ACTB

So far, we have evaluated reference genes at the mRNA level using RT-qPCR. We sought to determine what happened at the protein level using three of these reference genes. Western blot analysis was used to assess the stability of the PPIA, GAPDH, and ACTB proteins in proliferating and senescent hAD-MSCs, as shown in Fig. 3a. Then, the gray scale analyses were performed to calculate the percentage of each protein's expression level in senescent hAD-MSCs relative to the expression level in the corresponding proliferating cells. Figure 3b shows that the relative level of PPIA expression in A29S, A30S, and A47S was 99.68%, 107.05%, and 112.47%, respectively. The same percentages for GAPDH were 112.82%, 133.71%, and 108.54%, respectively (Fig. 3c). Figure 3d shows that the results for ACTB were 81.64%, 158.81%, and 143.49%, respectively.

FIG. 3.

Western blot analysis of PPIA, GAPDH, and ACTB protein levels. (a) Representative western blots showing the expression levels of PPIA, GAPDH, and ACTB in proliferating hAD-MSCs (A29P, A30P, and A47P) and senescent hAD-MSCs (A29S, A30S, and A47S). The three bar charts represent the results of western blot gray scale analysis and indicate the expression ratios of PPIA (b), GAPDH (c), and ACTB (d) in senescent cells (black bars) to their expression levels in the corresponding proliferating cells (white bars). All data are represented as the mean ± SD (n = 3). hAD-MSCs, human adipose-derived mesenchymal stem cells; SD, standard deviation.

Furthermore, we also calculated the interindividual variability in the expression levels of the three proteins in the different proliferating samples (Supplementary Fig. S3). Thus, these results indicate that PPIA showed the highest stability in protein expression, while the performance of GAPDH was ordinary, and ACTB showed the lowest stability. These results are consistent with the results of the analysis performed using the three algorithms.

Effect of different choice of reference genes

After confirming the stability of the reference genes, we sought to determine whether an evaluation of senescence markers would produce different results if they were normalized to stable versus unstable reference genes in RT-qPCR.

Therefore, we selected five widely used senescence-associated markers: (1) P16^INK4a , which is encoded by the INK4a/ARF locus that is correlated with the chronological age of nearly all tissues in both mice and humans [28 –30]; (2) P21^CLP1 , which is, like P16, a CDK inhibitor and has also been extensively reported to increase with cellular senescence in aging organisms in both mice and humans [31,32]; (3) LMNB1, a protein component of nuclear lamina; this is a biomarker of senescence processes that declines in expression during cell senescence [33,34]; (4) poly-ADP-ribose polymerase 1 (PARP1), a factor that has been reported to be upregulated with age and can impair DNA repair signaling [35 –37]; (5) and ISG15, which has been shown to be significantly lowered in expression as telomere position changes: as telomeres become shorter during aging, chromosome looping and gene expression patterns change and related genes, such as ISG15, are correspondingly regulated [38].

We analyzed the expression levels of the five senescence markers in all of the 14 samples (proliferating and senescent MSCs) with RT-qPCR, and the results were normalized to PPIA, GAPDH, and ACTB, respectively. As shown in Fig. 4 and Supplementary Fig. S4, normalization to different reference genes resulted in different relative fold expression of the target genes. Using unsuitable reference genes led to the false evaluation of target genes, an effect that was particularly exemplified in A47 cells, as shown in Fig. 4.

FIG. 4.

Effect of choice of different reference genes. The mRNA levels of P16^INK4a (a), p21^CLP1 (b), LMNB1 (c), PARP1 (d), and ISG15 (e) in senescent hAD-MSCs (A47P) are presented relative to their expression levels in corresponding proliferating cells (A47S). The values were normalized to PPIA, GAPDH, or ACTB, as indicated. All data are shown as the mean ± SD (n = 3). P-values were obtained using unpaired t-tests. ***P ≤ 0.001, **P ≤ 0.01, and *P < 0.05, and “ns” indicates P ≥ 0.05. hAD-MSCs, human adipose-derived mesenchymal stem cells; SD, standard deviation.

When normalized to PPIA, GAPDH, or ACTB, the expression level of P16^INK4a was higher in A47S (senescent stage) cells, by 9.04-, 8.26-, and 4.65-fold, respectively, than the expression levels observed in A47P (proliferating stage) cells (Fig. 4a). When normalized using PPIA or GAPDH, the expression of P21^CLP1 was 1.63- or 1.49-fold higher in A47S cells, whereas its level was reduced to 0.83-fold when normalized using ACTB (Fig. 4b). In A47S cells, the expression of LMNB1 was reduced to 0.15-, 0.13-, or 0.08-fold, respectively (Fig. 4c); the expression of PARP1 was reduced to 0.56-, 0.52-, or 0.29-fold, respectively (Fig. 4d); and the expression of ISG15 was increased to 2.08-, 1.90-, or 1.07-fold, respectively (Fig. 4e).

Notably, as shown in Fig. 4b and e, the senescent state of the cells was not recognized (ISG15) or even refuted (P21^CLP1 ) when gene expression was normalized using ACTB. Therefore, choosing a stable reference gene is very important when evaluating senescence using RT-qPCR analyses.

Discussion

The properties, potencies, and effects of senescent MSCs need to be explored to determine the potential of applications that involve this cell type in regeneration medicine and aging studies. Before initiating a study of senescence mechanisms, it is essential to screen for senescence-associated genes. RT-qPCR is considered the most efficient quantitative method for assessing RNA abundance. However, it requires assays that are carefully designed and that include strategies to optimize reactions to maximize their sensitivity and accuracy [39,40]. In earlier studies, the normal strategy for RT-qPCR normalization was to use a single housekeeping gene, usually GAPDH or ACTB, without further validation.

However, many reports have suggested that the expression of these genes can vary up to 10-fold across different samples [41,42]. In 2002, Vandesompele et al. suggested that experiments to determine the ideal reference genes should be carefully performed for each individual system by evaluating several genes [24].

In this study, we aimed to identify the optimal references genes to use in studies of senescence in hAD-MSCs. Cells were isolated from seven human donors with different age and both genders. To our knowledge, this is the first study to systematically determine the optimal references genes for a MSC-related senescence study. To obtain comparable results, we designed a rigorous set of parallel experiments that used the same conditions. We obtained data that displayed a high level of consistency and then analyzed the data using three independent expression stability analyses (geNorm, NormFinder, and BestKeeper). Despite the differences between these distinct statistical algorithms, PPIA was consistently ranked as the most stable reference gene, while ACTB, one of the most widely used reference genes in aging studies, was ranked last. GAPDH, another traditional reference gene, ranked in the middle.

Meanwhile, these results are consistent with results that were reported in 2013 and showed that “ACTB was the most unstable gene when evaluating adipose tissue-derived mesenchymal stromal cells” [43] and also the results reported in 2015 showing that “beta-2-microglobulin and peptidyl-prolyl isomerase A were the optimal reference genes for normalizing RT-qPCR data obtained from four different human MSCs” [44].

We considered the fact that ACTB encodes a cytoskeletal protein that is essential for cytoskeletal structure and kinetics. It should therefore not be surprising that its expression might be influenced by cellular senescence because senescence is a process that is characterized by morphological changes that result from changes in cellular structure and motility. A correlation between the actin cytoskeleton and senescence in human cells has been noted in some reports [45 –47], but it has not yet been adequately studied. GAPDH is a glyceraldehyde-3-phosphate that is involved in dehydrogenase activity and energy metabolism that has been reported to be associated with senescence in cancer cells [48 –50]. Finally, PPIA encodes a member of the peptidyl-prolyl cis-trans isomerase (PPIase) family and is involved in protein folding, but its role in senescence has not yet been explored.

To further validate the effect of reference genes, we assessed the expression of five senescence-associated genes (P16^INK4a , P21^CLP1 , LMNB1, PARP1, and ISG15) in both proliferating and senescent hAD-MSCs with normalization to PPIA, GAPDH, and ACTB, respectively. The results of assays in which expression was normalized to ACTB were inconsistent and showed statistically significant variation.

These results indicate that choosing an unsuitable reference gene can yield misleading results, and this is especially true when the expression levels of target genes, such as P21^CLP1 , PARP1, and ISG15 in our study, are not markedly altered under different conditions. In addition, as has been shown, in stem cells, slight changes in crucial protein factors can have effects that are strong enough to alter a cell's stage. Therefore, using a stable reference gene is of substantial importance for obtaining reliable experimental results in studies of stem cell aging.

In conclusion, our findings suggest that in senescence studies using hAD-MSCs replicative-senescence cell model, PPIA is the optimal reference gene to use for normalization in RT-qPCR assays. Our findings contribute to basic research in this field in addition to translational research into clinical cell therapies for aging-related diseases.

Footnotes

Acknowledgments

This study was supported by grants from the National Key Scientific Program of China (no. 2011CB964901), the National Collaborative Innovation Program, and the National Science and Technology Major Projects for “Drug Research and Development” (no. 2014ZX09101042). The authors also thank Stephen Dalton (University of Georgia) for reading and recommending this article.

Author Disclosure Statement

No competing financial interests exist.

References

Sousa-Victor

, Gutarra

, Garcia-Prat

, Rodriguez-Ubreva

, Ortet

, Ruiz-Bonilla

, Jardi

, Ballestar

, Gonzalez

, et al. (2014). Geriatric muscle stem cells switch reversible quiescence into senescence. Nature, 506:316–321.

Lopez-Otin

, Blasco

, Partridge

, Serrano

and Kroemer

. (2013). The hallmarks of aging. Cell, 153:1194–1217.

Bernet

, Doles

, Hall

, Kelly Tanaka

, Carter

and Olwin

. (2014). p38 MAPK signaling underlies a cell-autonomous loss of stem cell self-renewal in skeletal muscle of aged mice. Nat Med, 20:265–271.

Cosgrove

, Gilbert

, Porpiglia

, Mourkioti

, Lee

, Corbel

, Llewellyn

, Delp

and Blau

. (2014). Rejuvenation of the muscle stem cell population restores strength to injured aged muscles. Nat Med, 20:255–264.

Signer

and Morrison

. (2013). Mechanisms that regulate stem cell aging and life span. Cell Stem Cell, 12:152–165.

Trounson

and McDonald

. (2015). Stem cell therapies in clinical trials: progress and challenges. Cell Stem Cell, 17:11–22.

Yao

, Zhou

, Wan

, Su

, Sun

, Mizuguchi

, Yoshioka

, Nakagawa

, Zhao

and Gao

. (2014). Polyethyleneimine-coating enhances adenoviral transduction of mesenchymal stem cells. Biochem Biophys Res Commun, 447:383–387.

Gazit

, Pelled

, Sheyn

, Kimelman

and Gazit

. (2014). In: Essentials of Stem Cell Biology. Lanza

and Atala

, eds. Academic Press, Cambridge, MA, pp. 255–266.

, Zhu

, Qu

, Li

, Lin

, Liao

, Wang

, Xu

and Zhao

. (2013). Stepwise differentiation of human adipose-derived mesenchymal stem cells toward definitive endoderm and pancreatic progenitor cells by mimicking pancreatic development in vivo. Stem Cells Dev, 22:1576–1587.

10.

Feng

, Han

, Li

, Wang

, Li

, Yao

and Zhao

. (2014). Generation of highly purified neural stem cells from human adipose-derived mesenchymal stem cells by Sox1 activation. Stem Cells Dev, 23:515–529.

11.

DelaRosa

and Lombardo

. (2010). Modulation of adult mesenchymal stem cells activity by toll-like receptors: implications on therapeutic potential. Mediators Inflamm, 2010:865601.

12.

Volarevic

, Nurkovic

, Arsenijevic

and Stojkovic

. (2014). Concise review: therapeutic potential of mesenchymal stem cells for the treatment of acute liver failure and cirrhosis. Stem Cells, 32:2818–2823.

13.

Cao

, Sun

, Liao

, Meng

, Han

and Zhao

. (2005). Human adipose tissue-derived stem cells differentiate into endothelial cells in vitro and improve postnatal neovascularization in vivo. Biochem Biophys Res Commun, 332:370–379.

14.

Fang

, Liao

, Shi

, Yang

and Zhao

. (2004). Multipotency of Flk1CD34 progenitors derived from human fetal bone marrow. J Lab Clin Med, 143:230–240.

15.

Liu

, Ren

, Ge

, Cheng

, Zenke

, Keating

and Zhao

. (2015). Sca-1+Lin-CD117- mesenchymal stem/stromal cells induce the generation of novel IRF8-controlled regulatory dendritic cells through Notch-RBP-J signaling. J Immunol, 194:4298–4308.

16.

Kfoury

and Scadden

. (2015). Mesenchymal cell contributions to the stem cell niche. Cell Stem Cell, 16:239–253.

17.

Farini

, Sitzia

, Erratico

, Meregalli

and Torrente

. (2014). Clinical applications of mesenchymal stem cells in chronic diseases. Stem Cells Int, 2014:306573.

18.

Zhang

, Li

, Suzuki

, Qu

, Wang

, Zhou

, Liu

, Ren

, Xu

, et al. (2015). Aging stem cells. A Werner syndrome stem cell model unveils heterochromatin alterations as a driver of human aging. Science, 348:1160–1163.

19.

Cieslik

, Trial

, Carlson

, Taffet

and Entman

. (2013). Aberrant differentiation of fibroblast progenitors contributes to fibrosis in the aged murine heart: role of elevated circulating insulin levels. FASEB J, 27:1761–1771.

20.

Walker

. (2002). Tech.Sight. A technique whose time has come. Science, 296:557–559.

21.

Galiveti

, Rozhdestvensky

, Brosius

, Lehrach

and Konthur

. (2010). Application of housekeeping npcRNAs for quantitative expression analysis of human transcriptome by real-time PCR. RNA, 16:450–461.

22.

Suzuki

, Higgins

and Crawford

. (2000). Control selection for RNA quantitation. Biotechniques, 29:332–337.

23.

Ragni

, Vigano

, Rebulla

, Giordano

and Lazzari

. (2013). What is beyond a qRT-PCR study on mesenchymal stem cell differentiation properties: how to choose the most reliable housekeeping genes. J Cell Mol Med, 17:168–180.

24.

Vandesompele

, De Preter

, Pattyn

, Poppe

, Van Roy

, De Paepe

and Speleman

. (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol, 3:RESEARCH0034.

25.

Zhou

, Zhang

, Xia

, Wang

, Chen

, Fang

and Fan

. (2014). Selection of suitable reference genes for normalization of quantitative real-time polymerase chain reaction in human cartilage endplate of the lumbar spine. PLoS One, 9:e88892.

26.

Andersen

, Jensen

and Orntoft

. (2004). Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res, 64:5245–5250.

27.

Pfaffl

, Tichopad

, Prgomet

and Neuvians

. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper—Excel-based tool using pair-wise correlations. Biotechnol Lett, 26:509–515.

28.

Jeck

, Siebold

and Sharpless

. (2012). Review: a meta-analysis of GWAS and age-associated diseases. Aging Cell, 11:727–731.

29.

Krishnamurthy

, Torrice

, Ramsey

, Kovalev

, Al-Regaiey

, Su

and Sharpless

. (2004). Ink4a/Arf expression is a biomarker of aging. J Clin Invest, 114:1299–1307.

30.

Ressler

, Bartkova

, Niederegger

, Bartek

, Scharffetter-Kochanek

, Jansen-Durr

and Wlaschek

. (2006). p16INK4A is a robust in vivo biomarker of cellular aging in human skin. Aging Cell, 5:379–389.

31.

Capparelli

, Chiavarina

, Whitaker-Menezes

, Pestell

, Hulit

, Ando

, Howell

, Martinez-Outschoorn

, Sotgia

and Lisanti

. (2012). CDK inhibitors (p16/p19/p21) induce senescence and autophagy in cancer-associated fibroblasts, “fueling” tumor growth via paracrine interactions, without an increase in neo-angiogenesis. Cell Cycle, 11:3599–3610.

32.

Baker

, Weaver

and van Deursen

. (2013). p21 both attenuates and drives senescence and aging in BubR1 progeroid mice. Cell Rep, 3:1164–1174.

33.

Freund

, Laberge

, Demaria

and Campisi

. (2012). Lamin B1 loss is a senescence-associated biomarker. Mol Biol Cell, 23:2066–2075.

34.

Shimi

, Butin-Israeli

, Adam

, Hamanaka

, Goldman

, Lucas

, Shumaker

, Kosak

, Chandel

and Goldman

. (2011). The role of nuclear lamin B1 in cell proliferation and senescence. Genes Dev, 25:2579–2593.

35.

Altmeyer

and Hottiger

. (2009). Poly(ADP-ribose) polymerase 1 at the crossroad of metabolic stress and inflammation in aging. Aging (Albany NY), 1:458–469.

36.

Noren Hooten

, Fitzpatrick

, Kompaniez

, Jacob

, Moore

, Nagle

, Barnes

, Lohani

and Evans

. (2012). Coordination of DNA repair by NEIL1 and PARP-1: a possible link to aging. Aging (Albany NY), 4:674–685.

37.

Zhang

, Xiong

and Cao

. (2014). Mechanisms controlling the smooth muscle cell death in progeria via down-regulation of poly(ADP-ribose) polymerase 1. Proc Natl Acad Sci U S A, 111:E2261–E2270.

38.

Robin

, Ludlow

, Batten

, Magdinier

, Stadler

, Wagner

, Shay

and Wright

. (2014). Telomere position effect: regulation of gene expression with progressive telomere shortening over long distances. Genes Dev, 28:2464–2476.

39.

Peters

, Helps

, Hall

and Day

. (2004). Real-time RT-PCR: considerations for efficient and sensitive assay design. J Immunol Methods, 286:203–217.

40.

Nolan

, Hands

and Bustin

. (2006). Quantification of mRNA using real-time RT-PCR. Nat Protoc, 1:1559–1582.

41.

Warrington

, Nair

, Mahadevappa

and Tsyganskaya

. (2000). Comparison of human adult and fetal expression and identification of 535 housekeeping/maintenance genes. Physiol Genomics, 2:143–147.

42.

Tricarico

, Pinzani

, Bianchi

, Paglierani

, Distante

, Pazzagli

, Bustin

and Orlando

. (2002). Quantitative real-time reverse transcription polymerase chain reaction: normalization to rRNA or single housekeeping genes is inappropriate for human tissue biopsies. Anal Biochem, 309:293–300.

43.

Amable

, Teixeira

, Carias

, Granjeiro

and Borojevic

. (2013). Identification of appropriate reference genes for human mesenchymal cells during expansion and differentiation. PLoS One, 8:e73792.

44.

, Yang

, Bai

, Xuan

and Wang

. (2015). Identification of appropriate reference genes for human mesenchymal stem cell analysis by quantitative real-time PCR. Biotechnol Lett, 37:67–73.

45.

Cho

, Ryu

, Oh

, Park

, Lee

, Kim

, Jang

and Park

. (2004). Morphological adjustment of senescent cells by modulating caveolin-1 status. J Biol Chem, 279:42270–42278.

46.

Chen

, Tu

, Catania

, Burton

, Toussaint

and Dilley

. (2000). Involvement of Rb family proteins, focal adhesion proteins and protein synthesis in senescent morphogenesis induced by hydrogen peroxide. J Cell Sci, 113 (Pt 22):4087–4097.

47.

Gourlay

, Carpp

, Timpson

, Winder

and Ayscough

. (2004). A role for the actin cytoskeleton in cell death and aging in yeast. J Cell Biol, 164:803–809.

48.

Phadke

, Krynetskaia

, Mishra

and Krynetskiy

. (2011). Accelerated cellular senescence phenotype of GAPDH-depleted human lung carcinoma cells. Biochem Biophys Res Commun, 411:409–415.

49.

Liao

, Hsu

, Chuah

, Lee

, Hu

, Huang

, Yang

and Chiu

. (2014). Radiation induces senescence and a bystander effect through metabolic alterations. Cell Death Dis, 5:e1255.

50.

Nicholls

, Pinto

, Li

, Wang

, Simpson

and Liu

. (2012). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) induces cancer cell senescence by interacting with telomerase RNA component. Proc Natl Acad Sci U S A, 109:13308–13313.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.21 MB

0.09 MB

0.07 MB

0.10 MB

0.02 MB