PPIA and YWHAZ Constitute a Stable Pair of Reference Genes during Electrical Stimulation in Mesenchymal Stem Cells

: Mesenchymal stem cells (MSCs) are multipotent adult stem cells with great potential in regenerative medicine. One method for stimulating proliferation and differentiation of MSCs is via electrical stimulation (ES). A valuable approach for evaluating the response of MSCs to ES is to assess changes in gene expression, relative to one or more reference genes. In a survey of 25 publications that used ES on cells, 70% selected GAPDH as the reference gene. We conducted a study to assess the suitability of six potential reference genes on an immortalized human MSC line following direct current ES at seeding densities of 5000 and 10,000 cells/cm 2 . We employed three methods to validate the most stable reference genes from qRT-PCR data. Our ﬁndings show that GAPDH and ACTB exhibit reduced stability when seeded at 5000 cell/cm 2 . In contrast, we found that the most stable genes across both plating densities and stimulation regimes were PPIA and YWHAZ. Thus, in ES gene expression studies in MSCs, we support the use of PPIA and YWHAZ as an optimal reference gene pair, and discourage the use of ACTB and GAPDH at lower seeding densities. However, it is strongly recommended that similar veriﬁcation studies are carried out based on cell type and different ES conditions.


Introduction
In the field of regenerative medicine, changes in gene expression are often studied to provide evidence that stem cells are undergoing the process of differentiation. Of the techniques used to gather these data, one of the most common is real-time quantitative reverse transcription polymerase chain reaction, qRT-PCR [1]. Analysis of qRT-PCR data using the ∆∆C T method requires a comparison to be made between the detection level of a gene of interest relative to one or more reference genes, also commonly referred to as housekeeping genes, that remain stable during the experimental procedure [2].
Mesenchymal stem cells (MSCs) are of particular interest to the research community due to their immunomodulatory characteristics and their multipotency, which enables differentiation down osteogenic, chondrogenic, adipogenic, and neurogenic lineages [3,4]. Despite the evident potential of MSCs, there are still many hurdles to overcome before their effective use in patients. One main issue is that on implantation, cells often migrate away from the repair site, and therefore efficient integration with the target tissue or organ proves problematic [5]. A promising method for solving this could be the application of electrical stimulation (ES) to cells in vitro or in vivo, which has been shown to stimulate cell proliferation [6][7][8][9], initiate cell migration [10][11][12][13] and accelerate differentiation [14][15][16][17].

Metabolic Activity
The stimulated cells underwent 30 or 60 m/d of ES over 5 days. The cellular metabolic activity was measured on day 5, prior to the final ES period, using the Deep Blue Cell Viability Kit (BioLegend, London, UK). A volume of 1 mL of the solution administered to each well was constituted of 10% of the 'Deep Blue' resazurin solution, and 90% of cell media. The incubation period was 1 h, after which the fluorescence intensity of 200 µL of supernatant was measured at excitation wavelength of 544 nm and an emission wavelength of 590 nm on a plate reader (FLUOstar ® OPTIMA, BMG LabTech, Ortenberg, Germany).

Primer Amplification Efficiency
The primer amplification efficiency (E) was calculated using the following equation: where m is the slope of the standard curve for which C T value (y-axis) is plotted against log 10 of the cDNA mass in ng. The standard curve parameters were calculated by linear regression in GraphPad Prism software.

Data Analysis
NormFinder v0.953 was downloaded from the following website: https://moma.dk/ normfinder-software (accessed on 18 November 2021) and was run through Visual Basic for Applications (VBA) for Microsoft Excel [24]. The ∆∆C T method was used to calculate the fold changes for each gene (Gene 1) versus the gene determined by NormFinder to be the most stable for each seeding density (Gene 2) across unstimulated control and ES groups.
The series of equations used were as follows: The BestKeeper software was downloaded from the following website: https://www. gene-quantification.de/bestkeeper.html, (accessed on 18 November 2021) and was also run on VBA on Microsoft Excel [25].
In the ∆∆C T method results, fold changes for each gene were statistically analyzed using the unpaired t-test, without adjustment for multiple comparisons. Where the control was compared with 30 and 60 m/d of ES, a two-way ANOVA was used Tukey's multiple comparisons test. The metabolic activity assay was also analyzed by a two-way ANOVA.

Primers Are Effective and Selective for Their Respective Genes
A survey was undertaken of 25 published journal articles showing gene expression data in electrically stimulated cells, and 70% of these used GAPDH as the reference gene, with the second most popular (12%) being ACTB (Supplementary Table S1). In addition to these two genes, we also chose to include in our study, 18S, B2M, PPIA and YWHAZ as additional potential candidate reference genes [19,26].
Firstly, we wanted to verify if the primers used for each gene (Table 1) were effective over a range of cDNA concentrations, and that the PCR products were therefore unique and specific. This was achieved by plotting the standard curves for the 6 reference genes, (Figure 1). The values from the linear regression of these standard curves were used to calculate the primer efficiencies ( Table 2). The efficiency values for each of the primers in this study fall between 93.25% and 106.10%, indicating a satisfactory level of amplification with individual gradients between −3.2 and −3.5 (Table 2). We studied the melt curves derived from the StepOne software v.2.3 in order to check that the PCR products were in fact a single double-stranded DNA product. The melt curves confirmed this with the presentation of a single solitary peak, indicating the amplification of a single product for each gene (Supplementary Figure S1).

Cells Remain Viable under ES
The immortalized human Y201 MSC line was selected for use in this study. We next chose to examine whether the ES regimes we decided to use on the MSCs were cytotoxic. The cellular metabolic activity was measured on day 5, prior to the final period of ES, Thus, we assessed the effect of both 30 and 60 min per day (m/d) of ES on cell viability at two seeding densities based on measurements of the cells' metabolic activity ( Figure 2). There was no statistical significance seen between the electrically stimulated cells and the control cells in the 5000 cells/cm 2 group; however, at 10,000 cells/cm 2 , 60 m/d of ES incited a significant increase in the metabolic activity. This effect could be due to oxidative stress resulting from the over-confluence observed in the 10,000 cells/cm 2 group. This chosen seeding density proved to be excessive due to the highly proliferative nature of the Y201 cell line we used for these experiments. Importantly, there were no decreases in metabolic activity that would have been indicative of cell death, and therefore ES cytotoxicity.

Cells Remain Viable under ES
The immortalized human Y201 MSC line was selected for use in this study. We next chose to examine whether the ES regimes we decided to use on the MSCs were cytotoxic. The cellular metabolic activity was measured on day 5, prior to the final period of ES, Thus, we assessed the effect of both 30 and 60 min per day (m/d) of ES on cell viability at two seeding densities based on measurements of the cells' metabolic activity ( Figure 2). There was no statistical significance seen between the electrically stimulated cells and the control cells in the 5000 cells/cm 2 group; however, at 10,000 cells/cm 2 , 60 m/d of ES incited a significant increase in the metabolic activity. This effect could be due to oxidative stress resulting from the over-confluence observed in the 10,000 cells/cm 2 group. This chosen seeding density proved to be excessive due to the highly proliferative nature of the Y201 cell line we used for these experiments. Importantly, there were no decreases in metabolic activity that would have been indicative of cell death, and therefore ES cytotoxicity.

The Effect of the Duration of ES on Reference Gene Activity
Part of this study was to determine if increasing the duration of the ES over five daily applications had an effect on the expression of the six potential reference genes in our study. The ∆∆C T method was used to calculate the fold changes, for which a "gold standard" reference gene should remain as close as possible to 1.0 under all experimental conditions versus the unstimulated control cells and another reference gene. Finding a reference gene pair is particularly advantageous as this corrects for any inaccurate RNA quantification steps and for differences in the efficiency of the reverse transcription process. Furthermore, it was necessary to assess if the fold changes after 60 m/d of ES were statistically different from that of the 30 m/d group or the control group. The fold changes of PPIA, 18S, and GAPDH were normalized to YWHAZ in unstimulated cells or in cells exposed to both durations of ES ( Figure 3). The increases observed in 18S after 30 and 60 m/d were not significant, nor was the rise in expression of GAPDH after 60 m/d. For all three of the reference genes tested, there were no significant differences in expression between the ES durations of 30 m/d and 60 m/d. The clear similarity observed in response for both the 30 m/d and 60 m/d groups provided the rationale to continue focusing solely on the expression changes in cells stimulated for 60 m/d. It was also assumed that a reference gene exhibiting stable, consistent expression after 60 m/d of ES would be likely to also remain stable under a shorter period of ES of the same electric field strength.

The Effect of the Duration of ES on Reference Gene Activity
Part of this study was to determine if increasing the duration of the ES over five daily applications had an effect on the expression of the six potential reference genes in our study. The ∆∆CT method was used to calculate the fold changes, for which a "gold standard" reference gene should remain as close as possible to 1.0 under all experimental conditions versus the unstimulated control cells and another reference gene. Finding a reference gene pair is particularly advantageous as this corrects for any inaccurate RNA quantification steps and for differences in the efficiency of the reverse transcription process. Furthermore, it was necessary to assess if the fold changes after 60 m/d of ES were statistically different from that of the 30 m/d group or the control group. The fold changes of PPIA, 18S, and GAPDH were normalized to YWHAZ in unstimulated cells or in cells exposed to both durations of ES ( Figure 3). The increases observed in 18S after 30 and 60 m/d were not significant, nor was the rise in expression of GAPDH after 60 m/d. For all three of the reference genes tested, there were no significant differences in expression between the ES durations of 30 m/d and 60 m/d. The clear similarity observed in response for both the 30 m/d and 60 m/d groups provided the rationale to continue focusing solely on the expression changes in cells stimulated for 60 m/d. It was also assumed that a reference gene exhibiting stable, consistent expression after 60 m/d of ES would be likely to also remain stable under a shorter period of ES of the same electric field strength.

The Effect of the Duration of ES on Reference Gene Activity
Part of this study was to determine if increasing the duration of the ES over five daily applications had an effect on the expression of the six potential reference genes in our study. The ∆∆CT method was used to calculate the fold changes, for which a "gold standard" reference gene should remain as close as possible to 1.0 under all experimental conditions versus the unstimulated control cells and another reference gene. Finding a reference gene pair is particularly advantageous as this corrects for any inaccurate RNA quantification steps and for differences in the efficiency of the reverse transcription process. Furthermore, it was necessary to assess if the fold changes after 60 m/d of ES were statistically different from that of the 30 m/d group or the control group. The fold changes of PPIA, 18S, and GAPDH were normalized to YWHAZ in unstimulated cells or in cells exposed to both durations of ES ( Figure 3). The increases observed in 18S after 30 and 60 m/d were not significant, nor was the rise in expression of GAPDH after 60 m/d. For all three of the reference genes tested, there were no significant differences in expression between the ES durations of 30 m/d and 60 m/d. The clear similarity observed in response for both the 30 m/d and 60 m/d groups provided the rationale to continue focusing solely on the expression changes in cells stimulated for 60 m/d. It was also assumed that a reference gene exhibiting stable, consistent expression after 60 m/d of ES would be likely to also remain stable under a shorter period of ES of the same electric field strength.

NormFinder Results
To prevent bias and to gain a more detailed analysis of the C T values, we used the NormFinder software v0.953. Stability values for each gene are determined via an algorithm. The lower the stability value, the higher the relative stability of the gene. The NormFinder results are presented for cells stimulated for 60 m/d at each seeding density in Figure 4. The results show that PPIA is the most stable gene for the 5000 cells/cm 2 group and ACTB is the least stable at this seeding density (Figure 4a). For the cells seeded at 10,000 cells/cm 2 , YWHAZ was found to be the most stable, and ACTB was again the least stable ( Figure 4b). GAPDH was ranked as the second most stable gene in both seeding density groups (Figure 4a,b).

NormFinder Results
To prevent bias and to gain a more detailed analysis of the CT values, we used the NormFinder software v0.953. Stability values for each gene are determined via an algorithm. The lower the stability value, the higher the relative stability of the gene. The NormFinder results are presented for cells stimulated for 60 m/d at each seeding density in Figure 4. The results show that PPIA is the most stable gene for the 5000 cells/cm 2 group and ACTB is the least stable at this seeding density (Figure 4a). For the cells seeded at 10,000 cells/cm 2 , YWHAZ was found to be the most stable, and ACTB was again the least stable (Figure 4b). GAPDH was ranked as the second most stable gene in both seeding density groups (Figure 4a,b). Interestingly, the stability of the genes studied increased consistently with increasing seeding density. As would be expected, this would imply that cells seeded at the lower density are more greatly affected by the ES potentially due to their heightened levels of exposure.

∆∆. CT Method Results
The ∆∆CT method calculates the fold change of a gene of interest versus a reference gene. Crucially, it is the only way to interpret qRT-PCR data from the raw threshold cycle (CT) values. This is why we considered it important to validate the reference genes using the same method, as this would be used when analysing if a gene of interest was changing under ES conditions. The fold changes for the reference genes were plotted versus the gene determined to be the most stable by NormFinder for each of the seeding densities using the ∆∆CT method ( Figure 5). Interestingly, the stability of the genes studied increased consistently with increasing seeding density. As would be expected, this would imply that cells seeded at the lower density are more greatly affected by the ES potentially due to their heightened levels of exposure.

∆∆C T Method Results
The ∆∆C T method calculates the fold change of a gene of interest versus a reference gene. Crucially, it is the only way to interpret qRT-PCR data from the raw threshold cycle (C T ) values. This is why we considered it important to validate the reference genes using the same method, as this would be used when analysing if a gene of interest was changing under ES conditions. The fold changes for the reference genes were plotted versus the gene determined to be the most stable by NormFinder for each of the seeding densities using the ∆∆C T method ( Figure 5).
The ∆∆C T results showed a significant increase in GAPDH expression while ACTB also exhibited an upwards trend in expression, relative to PPIA at 5000 cells/cm 2 (Figure 5a). The p-value for ACTB was 0.06, therefore it was close to being significant, but failed to reach significance possibly because of the relatively large standard deviation in the values.
In contrast, YWHAZ, 18S, and B2M did not show any significant differences in the ES groups versus PPIA and the control cells (Figure 5a). At the seeding density of 10,000 cells/cm 2 , the expression of none of the potential reference genes were significantly altered following ES (Figure 5b). However, the most consistent pair were PPIA and YWHAZ, with an average fold change in the ES group of 1.08 and 0.90 respectively (Figure 5a,b).
Our results suggest that GAPDH may be a poor choice for a reference gene at seeding densities below 10,000 cells/cm 2 . The significant increase noted in GAPDH contradicts the NormFinder results, which placed GAPDH as the second most stable gene consistently at both seeding densities. This highlights the importance of using more than one method to . Gene expression analysis performed using the ∆∆CT method in unstimulated controls (Control) and cells stimulated for 60 min per day (+ES). Expression is normalized to the control and to the gene determined to be the most stable for the respective seeding density by NormFinder. Each graph corresponds to the two seeding densities used: (a) 5000 cells/cm 2 , and (b) 10,000 cells/cm 2 . * p < 0.05 versus control of the same gene. Each gene was analyzed separately using an unpaired T-test; N = 2 for control, N = 3 for +ES group.
The ∆∆CT results showed a significant increase in GAPDH expression while ACTB also exhibited an upwards trend in expression, relative to PPIA at 5000 cells/cm 2 (Figure 5a). The p-value for ACTB was 0.06, therefore it was close to being significant, but failed to reach significance possibly because of the relatively large standard deviation in the values.
In contrast, YWHAZ, 18S, and B2M did not show any significant differences in the ES groups versus PPIA and the control cells (Figure 5a). At the seeding density of 10,000 cells/cm 2 , the expression of none of the potential reference genes were significantly altered following ES (Figure 5b). However, the most consistent pair were PPIA and YWHAZ, with an average fold change in the ES group of 1.08 and 0.90 respectively (Figure 5a,b).
Our results suggest that GAPDH may be a poor choice for a reference gene at seeding densities below 10,000 cells/cm 2 . The significant increase noted in GAPDH contradicts the NormFinder results, which placed GAPDH as the second most stable gene consistently at both seeding densities. This highlights the importance of using more than one method to determine reference gene stability. The selection of a poor reference gene could easily mask a significant upregulation of a gene of interest.
Despite the lack of significance, ACTB may also not be an appropriate reference gene for studies involving ES due to its observed increase in expression and high variability at 5000 cells/cm 2 . It is important to recognise that this study only focused on an electric field strength of 100 mV/mm and a direct current set-up in a 2D culture environment. Within Figure 5. Gene expression analysis performed using the ∆∆C T method in unstimulated controls (Control) and cells stimulated for 60 min per day (+ES). Expression is normalized to the control and to the gene determined to be the most stable for the respective seeding density by NormFinder. Each graph corresponds to the two seeding densities used: (a) 5000 cells/cm 2 , and (b) 10,000 cells/cm 2 . * p < 0.05 versus control of the same gene. Each gene was analyzed separately using an unpaired T-test; N = 2 for control, N = 3 for +ES group.
Despite the lack of significance, ACTB may also not be an appropriate reference gene for studies involving ES due to its observed increase in expression and high variability at 5000 cells/cm 2 . It is important to recognise that this study only focused on an electric field strength of 100 mV/mm and a direct current set-up in a 2D culture environment. Within the system, there are many variables at play which all could influence the reference gene activity. Based on this precise experimental set-up, ACTB and GAPDH will be avoided in future, especially at lower seeding densities. In contrast, PPIA and YWHAZ would appear to be a stable pair regardless of the number of cells seeded prior to ES.

BestKeeper Results
In the literature, studies on the most suitable reference genes employ multiple platforms to compare the findings and ultimately verify their conclusions. We, therefore, decided to also input our raw C T values into BestKeeper [25]. The Pearson's correlation coefficient (r) value results generated by BestKeeper were studied, and are shown in Table 3. Each potential gene pair has been ranked in order of the highest average r value over both seeding densities. As an r value of one illustrates a highly positive linear correlation, PPIA and YWHAZ constituted the gene pair with the strongest relationship, and also demonstrated the most consistent pairing with the lowest standard deviation. ACTB consistently appeared in the lowest ranked gene pairings, with the lowest average r values and the highest standard deviations. Next, four out of the possible five pairings including GAPDH were ranked sequentially in the middle of the dataset. However, these r values were still indicative of a strongly positive correlation between the pairs and had a low standard deviation across the seeding densities. An r value above 0.9 would still be determined to be a highly positive correlation, and would indicate that PPIA and GAPDH would make a suitable reference gene pair; however, the significant increases in GAPDH as shown by the ∆∆C T show the innate value in using multiple methods to gain a full and accurate picture of the reference gene expression. The in-depth descriptive statistics of each data set can be found in Supplementary Tables S2 and S3.

PPIA and YWHAZ Expression in Primary Human MSCs
As previously mentioned in Results Section 3.2, an immortalized human MSC line was used in this study. A question we posed was if the findings were also consistent in primary human MSCs, that still retain replicative senescence and are of a different donor. Therefore, we assessed the expression of YWHAZ versus PPIA in primary human bone marrow-derived MSCs ( Figure 6). We found no change in the gene expression after 30 or 60 m/d of ES; therefore, YWHAZ and PPIA were also stable with increasing duration of ES in primary MSCs as well as the Y201 MSC line. This shows that the human immortalized cells used in this study are a highly relevant MSC line, with transferrable results. was used in this study. A question we posed was if the findings were also consistent in primary human MSCs, that still retain replicative senescence and are of a different donor. Therefore, we assessed the expression of YWHAZ versus PPIA in primary human bone marrow-derived MSCs ( Figure 6). We found no change in the gene expression after 30 or 60 m/d of ES; therefore, YWHAZ and PPIA were also stable with increasing duration of ES in primary MSCs as well as the Y201 MSC line. This shows that the human immortalized cells used in this study are a highly relevant MSC line, with transferrable results.

Discussion
In general, most gene expression studies use only one reference gene, rather than a pair of reference genes (Supplementary Table S1). Using pairs of reference genes is the gold standard because any inaccuracies in the relative concentrations of cDNA can be adjusted for in further analysis. Not only is RNA quantification often prone to error, but also the RNA-to-cDNA reverse transcription efficiency can be variable between samples. In addition, there is evidence that the efficiency of reverse transcription may be gene dependent as well [27]. Interestingly, our work has highlighted that ACTB may be problematic in showing high variability between replicates, possibly due to variability in its ability to be reverse transcribed. The input concentration of cDNA is calculated by assuming a 'perfect' reverse transcription efficiency of 100% for all genes [27]; however, the actual efficiency is unknown and will differ, as already mentioned, between genes and between samples.
Knowing therefore that two reference genes are strongly correlated, and using them both as a baseline, gives further confidence that any differences assessed between the reference genes and the genes of interest will be genuine. It is proposed that the CT values of PPIA and YWHAZ be averaged and then compared with that of a gene of interest. A further advantage of using a pair of reference genes is the ability to quickly identify a scenario in which either gene does change. A similar analysis, as carried out in this study, could then be conducted to find a new, more suitable pair.
Many other published reference or 'housekeeping' gene articles use BestKeeper, NormFinder, and also GeNorm [2,19,22]. It is acknowledged that a further platform could add greater accuracy in the determination of the most stable reference genes; however, this was achieved using the ∆∆CT method and stable genes were still able to be confidently identified. Furthermore, it has previously been shown that evaluating the analysis from multiple approaches provides a more accurate estimate of the best reference genes; GeNorm was also determined to be less rigorous than NormFinder, particularly where some genes had a high variability across samples [28].

Discussion
In general, most gene expression studies use only one reference gene, rather than a pair of reference genes (Supplementary Table S1). Using pairs of reference genes is the gold standard because any inaccuracies in the relative concentrations of cDNA can be adjusted for in further analysis. Not only is RNA quantification often prone to error, but also the RNA-to-cDNA reverse transcription efficiency can be variable between samples. In addition, there is evidence that the efficiency of reverse transcription may be gene dependent as well [27]. Interestingly, our work has highlighted that ACTB may be problematic in showing high variability between replicates, possibly due to variability in its ability to be reverse transcribed. The input concentration of cDNA is calculated by assuming a 'perfect' reverse transcription efficiency of 100% for all genes [27]; however, the actual efficiency is unknown and will differ, as already mentioned, between genes and between samples.
Knowing therefore that two reference genes are strongly correlated, and using them both as a baseline, gives further confidence that any differences assessed between the reference genes and the genes of interest will be genuine. It is proposed that the C T values of PPIA and YWHAZ be averaged and then compared with that of a gene of interest. A further advantage of using a pair of reference genes is the ability to quickly identify a scenario in which either gene does change. A similar analysis, as carried out in this study, could then be conducted to find a new, more suitable pair.
Many other published reference or 'housekeeping' gene articles use BestKeeper, NormFinder, and also GeNorm [2,19,22]. It is acknowledged that a further platform could add greater accuracy in the determination of the most stable reference genes; however, this was achieved using the ∆∆C T method and stable genes were still able to be confidently identified. Furthermore, it has previously been shown that evaluating the analysis from multiple approaches provides a more accurate estimate of the best reference genes; GeNorm was also determined to be less rigorous than NormFinder, particularly where some genes had a high variability across samples [28].
Of the primers used, three out of the six had efficiencies of over 100%; this means that the rate of doubling was slightly above two per cycle [29]. The ∆∆C T method calculations rely on the assumption that each cycle is accompanied by a doubling of the number of target sequence molecules; therefore efficiencies of less than 90% are more problematic as the molecules are not achieving said threshold [30,31].
The seeding densities of 5000 and 10,000 cells/cm 2 have both commonly been used in previous ES studies [13,14,[32][33][34]. The electric field strength was selected primarily because the same ES device had also shown pro-osteogenic effects at 100 mV/mm in the literature, and it falls within the range of physiological EFs that are established in vivo for both development and regeneration [14,35]. The metabolic activity differences at the higher seeding density do not constitute an unusual finding; increases were observed in rat adipose-derived MSCs in the same ES set-up and in cells seeded on a substrate formed of hydroxyapatite, respectively [14,34]. Evidently, the seeding density plays a role in the cells' response to ES. However, a full analysis of this effect falls out of the scope of this study and will be appropriately investigated in future work.
In a survey taken of published literature using ES, 70% of the papers cited GAPDH as the reference gene in the study. Based on our findings, the use of GAPDH as the only reference gene in ES studies should be verified before continued use, especially at lower seeding densities. However, it is equally acknowledged that the field of ES research uses a vast array of parameters, different cell types, as well as varying modalities for delivering the stimulation to the cells such as DC ES, AC ES, a capacitive set-up, or even electromagnetic stimulation. Even a factor as seemingly small as the electrode material can have a large impact on the cellular response due to differences in faradic by-products [36]. Therefore, the results in this study only serve as an indication of how MSCs react in the system described. Any decisions to use PPIA and YWHAZ as a reference gene pair should be first validated in each specific ES set-up.
Hydrogen peroxide (H 2 O 2 ) is produced at the cathode during DC ES, and therefore, the levels of exposure per cell to the H 2 O 2 concentration would be elevated when fewer cells are present in the well. The reason for GAPDH's increase in expression at the lower seeding density of 5000 cells/cm 2 could be linked to findings that GAPDH is involved in regulating intracellular redox levels [37].
A possible explanation for the lack of stability of ACTB could be due to its role in cell migration [38]. MSCs in a direct current electric field migrate towards the anode at electric field strengths as low as 25 mV/mm [13]. Consequently, this could provide rationale for ACTB's involvement in the cellular response.
All in all, 18S, while stable, has drawbacks due to the fact that it is expressed at much higher levels than the average gene of interest, appearing after 10-12 cycles. Due to the difference between the gene of interest and the reference gene being an exponent of 2 in fold change calculations, small differences are highly sensitive, and errors can be reduced by selecting a reference gene with a C T value more similar to that of the genes of interest in the study. Moreover, due its single exon any primers targeting 18S cDNA could also amplify the corresponding genomic region. Thus, contaminating DNA might skew the findings unless removed prior to reverse transcription.
Overall, comparing the findings from NormFinder, the ∆∆C T method, and BestKeeper has proven the value of using multiple platforms. As much as NormFinder identified ACTB as the least stable gene at both seeding densities, it consistently placed GAPDH as the second most stable gene. This exposes a flaw in solely relying upon NormFinder, as when the fold changes were calculated, GAPDH significantly increased in the 5000 cells/cm 2 group. The Pearson correlation coefficient analysis from BestKeeper ranked GAPDH mostly in the middle of the dataset. This further highlights the importance of selecting a stable reference gene pair. A poor reference gene has the potential of giving inaccurate gene expression results, which could confound and delay the advancement of knowledge within the research community. Thus, the value in carefully choosing the reference gene pairs used in gene expression studies cannot be underestimated.
In summary, our study identified PPIA and YWHAZ as the most stable reference genes following ES. These two genes were similarly found to be the most stable following long-term MSC culture [19]. The ∆∆C T method showed an undeviating similarity in expression of both genes. This was further reinforced by the BestKeeper Pearson's correlation coefficient data, which ranked PPIA and YWHAZ as the best potential gene pair overall with a very highly positive and consistent relationship at both seeding densities. This finding was also replicated in primary MSCs, which provides confidence that the conclusions from this study are not only applicable to the immortalized MSCs, but to primary cells from a different donor.
Due to the wide variation in parameters, conditions, electrical and temporal regimes applied by numerous research groups, this paper cannot state that PPIA and YWHAZ are unilaterally the most stable reference genes under ES. However, it is hoped that this study will encourage researchers to thoroughly assess the transcriptional response of commonlyused reference genes across different time points, varying electric field strengths, and duration of stimulation prior to investigating pathways of interest.