Real-time reverse transcription quantitative polymerase chain reaction (RT-qPCR) is a widely used method to measure relative gene expression for biological and medical research, as well as for clinical purposes, such as diagnostics, disease detection and monitoring [1]. To achieve a high level of data reliability and avoid inconsistency in data presentation and interpretation, the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines were developed and published [2]. In RT-qPCR assay the target gene expression should be normalized to the reference gene expression, and the reference gene expression must remain stable across all conditions tested. As housekeeping gene encoded products are involved in maintaining essential cellular functions, these genes are expressed in the majority of cells and tissues and are usually used as reference genes [3]. However, even housekeeping gene expression may dramatically differ depending on tissue type, experimental conditions and treatments. Thus, candidate reference gene stability must be examined for each experimental model to select the most stable genes as suitable references [4]. This selection is particularly important if the study aims to compare the expression levels of target genes in different cell types.

Several mathematical algorithms have been developed to evaluate reference gene stability, and the most frequently used are BestKeeper [5], NormFinder [6], geNorm [7] and the comparative Delta Ct ($\mathrm{\Delta }$Ct) method [8]. GeNorm evaluates the average pairwise variation of each gene, compares it with other analysed genes and calculates the expression stability value M [7]. The lowest M reflects the least pairwise variation that corresponds to the most stable gene in experimental model testing [9]. NormFinder calculates the stability value (SV) for all candidate reference genes tested, with a low SV value indicating stable gene expression [6]. The $\mathrm{\Delta }$Ct method compares the relative expression of ‘pairs of genes’ within each sample and thus identifies the most stable housekeeping genes. The gene with the lowest mean $\mathrm{\Delta }$Ct standard deviation is considered to be the most stable [8]. BestKeeper performs pairwise comparisons of the raw cycle of quantification (Cq) values for each gene, calculating an index using the geometric mean of each candidate gene’s Cq value, as well as the coefficient of variation (CV), the standard deviation (StDev) of the Cq values and the coefficient of correlation (r). The stable reference genes have high r values and low CVs and StDevs [5]. Finally, the RefFinder web-based tool calculates gene stability rankings through geometric averaging ranks obtained with these four methods [9].

When investigating various aspects of signalling and the vesicular transport of transmembrane receptors, researchers have typically used both stable immortalised cancer cell lines and primary cultures, such as mesenchymal stromal cells (MSCs), the latter being model cells with intact regulation of the studied processes [10, 11, 12]. In this study, on the basis of the aforementioned algorithms, we aimed to choose the most stable references among eight commonly used housekeeping genes (GAPDH, POP4, 18S, YWHAZ, HPRT, RPL13A, B2M, and EIF2B1) for a set of human normal, non-malignant cells and pathological cancer cell lines. We sequentially evaluated reference gene stability in the following:

(1) Normal and cancer cell lines.

(2) Cancer cell lines only.

(3) Normal mesenchymal cell lines only.

We found significant differences in housekeeping gene expression between these three groups, as well as less pronounced variability in gene expression within normal mesenchymal cell lines. Reference gene stability analysis identified POP4 and EIF2B1 as the most stable genes within the tested cell panel, while the POP4/GAPDH and 18S/B2M gene pairs appeared to be most suitable as references for cancer cell lines only and normal mesenchymal cell lines only, respectively. We also validated our choice of reference genes for the set of human cell lines by analysing the relative expression of epidermal growth factor receptor (EGFR)—the main component of the EGFR-signalling network, which is commonly dysregulated in cancers and governs cancer progression. Our results are useful for guiding the selection of suitable candidate reference genes for RT-qPCR in further comparative studies of cancer and normal cells.

We analysed the stability of eight housekeeping genes in a wide variety of human cell lines, which can be subdivided into two main subgroups:

(1) Human malignant cells (cervix carcinoma HeLa cell line, breast adenocarcinoma MCF-7, uterine leiomyosarcoma SK-UT-1B, lung carcinoma cell line A549, epidermoid carcinoma A431, breast adenocarcinoma cell line SK-BR-3).

(2) Human non-malignant (normal) mesenchymal cells: endometrial mesenchymal stromal cells (enMSCs, line 2804 and AMO), human umbilical cord MSCs (MSCWJ-1) and human dental pulp MSCs (MSC-DP).

We compared the distribution of the Cq values of the reference genes for all the cells tested (Fig. 2). The 18S gene showed the highest expression (Cq = 13.1–18.0, median of 14.6), while EIF2B1 had the lowest expression (Cq = 24.7–31.1, median of 27.5). The Cq values (median meanings) of other housekeeping genes ranged from 21 to 26.

We first evaluated the stability of the candidate reference genes in the entire panel of cell lines. The reference gene stability for the panel was obtained using four algorithms (the comparative $\mathrm{\Delta }$Ct method, BestKeeper, NormFinder, and geNorm) and ranked with the RefFinder online tool (https://www.ciidirsinaloa.com.mx/RefFinder-master/) (Fig. 3).

According to the RefFinder comprehensive stability ranking, POP4, 18S and YWHAZ were the three most stable housekeeping genes in this model, while the three least stable genes were HPRT, RPL13A and B2M.

Based on this ranking, we excluded the least stable reference genes to avoid possible bias in calculations from unstably expressed genes. The most stable reference genes as determined by RefFinder were POP4 and EIF2B1 (Fig. 3). The five most successful candidate reference genes were individually analysed in detail by each of the four algorithms. According to geNorm, all five genes tested had stability values lower than the cutoff value of M = 1.5, indicating that they were not unreliable [7]. The lowest M value was 0.667 for EIF2B1 and pOP4, so these genes can be interpreted by geNorm as the most stable ones. According to the BestKeeper algorithm, only the 18S gene had an StDev value less than 1 (0.723) and thus could be considered stable [5]. Both the NormFinder and $\mathrm{\Delta }$Ct methods revealed POP4, GAPDH and YWHAZ as the most stable housekeeping genes.

Next, we analysed the stability of housekeeping genes in the subgroups of cancer and normal mesenchymal cell cultures separately.

First, we compared the RefFinder comprehensive ranking results for these two subgroups (Figs. 4,5). In the cancer cell lines, POP4, GAPDH and EIF2B1 were the most stable genes, and B2M, RPL13A and HPRT1 were the least stable. This differed from the ranking of housekeeping genes in the subgroup of mesenchymal stromal cell cultures (MSCs), where none of the POP4, GAPDH or EIF2B1 genes were among the most stable ones. Moreover, GAPDH in MSCs had the highest ranking value. In MSCs, RefFinder indicated that 18S, B2M and RPL13A were the most stable genes (Fig. 5), while GAPDH, HPRT1 and YWHAZ were identified as the most variable genes.

The geNorm M value varied from 0.674 for EIF2B1 and POP4 to 1.469 for RPL13A in the cancer cells (Fig. 4) and from 0.429 for 18S and RPL13A to 0.989 for GAPDH in MSCs (Fig. 5). Thus, the M values in both subgroups were lower than the cutoff stability value M = 1.5, and these values were much lower for MSCs. According to BestKeeper, 18S and RPL13A (Fig. 4) in cancer cells and 18S, B2M, RPL13A and POP4 in MSCs had StDev values 1 (Fig. 5).

Both the NormFinder and $\mathrm{\Delta }$Ct methods revealed GAPDH and POP4 as the least variable housekeeping genes in cancer cell lines, while RPL13A was indicated as the most variable (Fig. 4).

The NormFinder and $\mathrm{\Delta }$Ct algorithms showed similar rankings when housekeeping gene expression was analysed in MSCs (Fig. 5). Namely, B2M, 18S and RPL13A were found to be the most stable, whileYWHAZ, HPRT1 and GAPDH were the least stable. The averaged ranking presented by RefFinder revealed 18S and B2M as the least variable genes, and YWHAZ and GAPDH showed the greatest variability (Fig. 5). For all four algorithms used, the stability values in the MSC subgroup were lower than those in the cancer cell subgroup, elucidating the extent of metabolic pathway variability within these subgroups of cell lines.

Finally, to demonstrate that the results concerning the target gene expression can differ depending on the reference genes chosen, we conducted an analysis of the expression of epidermal growth factor receptor (EGFR), a transmembrane protein involved in the regulation of cell survival, growth, proliferation and differentiation. As the dysregulation of EGFR expression and signalling is closely linked to cancer progression [25], investigations of EGFR-mediated signal transduction and intracellular transport peculiarities in cancer versus normal (non-malignant) cells are of considerable interest to many researchers.

First, the expression data were normalised by the geometric average of the three most stable reference genes selected according to RefFinder ranking (POP4, EIF2B1, YWHAZ). We found that EGFR mRNA was decreased in the MCF-7 and SK-UT-1B cell lines and increased in the A431 line compared to HeLa (Fig. 6A). We did not detect any differences in EGFR mRNA levels between HeLa and the other cell lines tested. However, when the results were normalised by the three least stable genes identified (HPRT1, RPL13A, B2M), the aforementioned changes were not detected. In this case, we did find a false increase in EGFR mRNA content in SK-BR-3 cells compared to HeLa cells (Fig. 6B). Thus, inaccurate normalisation leads to significant misinterpretations of EGFR expression data.

In this study, we performed an RT-qPCR analysis of the expression stability of eight candidate reference genes from the variety of cancer and normal cell lines frequently used in cell studies. We selected the most widely used human cancer cell lines (HeLa, MCF-7, SK-UT-1B, A549, A431 and SK-BR-3) as well as several primary mesenchymal cell lines as the model of normal cells (two lines of human endometrial MSCs, one line of human umbilical cord MSCs [MSCWJ-1] and one line of dental pulp MSCs [MSC-DP]). We deliberately avoided the use of in vitro immortalised cell lines as a model of normal cells in our study, as the gene expression profile of such cells may significantly differ from that of primary cells of the same origin [26].

Our results provide technical information about the reference gene expression stability, and we don’t aim to distinguish the mechanisms affecting variability in gene expression in the experiments. The observed changes between the cell line panels could have originated from the changes in transcriptional and posttranscriptional regulation at the subcellular level, as well as from the differences in cell subtype composition within the samples. However, in this section we attempt to compare our data with those from previous research to discuss possible factors affecting observed differences in the reference gene expression stability.

RefFinder identified the POP4 gene as the most stable in the group of cancer cells and in the entire panel of cells. The POP4 gene (processing of precursor 4, ribonuclease P/MRP subunit) encodes RPP29, a protein subunit of the two complexes: the endoribonuclease for the mitochondrial RNA processing complex involved in mitochondrial RNA editing and the ribonuclease P complex. In the latter complex, RPP29 acts as the central protein that interacts with other proteins of the complex and with the RNA included in the ribonuclease P (RNase P) complex [27]. Ribonuclease P is an ancient enzyme found in all kingdoms of life that is involved in the maturation of the 5^′ end of transfer RNA (tRNA) [28]. RNase P and its components are well known to have non-canonical functions in regulating chromatin assembly and DNA damage responses, and they process basic biomolecules, such as pre-tRNA and other RNA substrates [29]. Thus, stable expression of the POP4 gene is likely critical to normal and especially malignant cells, at least in vitro, due to its multiple physiological functions in RNA processing and genome regulation.

According to our data, the EIF2B1 gene was also among the most stable housekeeping genes in a wide panel of cell lines of different origins. This gene encodes one of the subunits of eukaryotic translation initiation factor 2B (eIF2B). eIF2B is the guanine nucleotide exchange factor for the eIF2 complex. eIF2 delivers an initiator, methionyl-tRNA, to the ribosome to promote translation initiation, and eIF2B regulates the activity of eIF2 by exchanging GDP for GTP on the $\gamma$-subunit of eIF2, thus activating eIF2 [30]. As eIF2 expression appears to be stable within the cancer cells, the EIF2B1 functional role is likely essential and strongly regulated even in malignant cell lines, despite their genome instability. Indeed, eIF2 is known to be involved in tumorigenesis, while the inhibition of the eIF2B subunit in vitro has an anticancerogenic effect. The regulation of mRNA translation initiation is a simpler and quicker step for controlling protein content in cells than a transcriptional step. Cancer cells require effective working machinery to regulate the mRNA translation initiation, as cancer progression is dependent on cancer cells’ opportunity to adapt rapidly to the changes of a microenvironment [31].

Different cellular investigations have often identified the EIF2B1 and POP4 genes as the most stable. EIF2B1 and POP4 were revealed to be one of the most stable genes among the 24 candidate genes tested in total RNA samples isolated from surgically removed tumor and non-malignant tissues of 10 patients with pancreatic carcinoma [32]. EIF2B1 and POP4 were also found to be the most stable of the 32 housekeeping genes when mesenchymal stem cells derived from the menstrual blood of women with endometriosis were compared with those of healthy women [16]. However, in cancerous and non-malignant tissues of the human uterine cervix, POP4 and EIF2B1 are among the most variable genes [33]. According to our data concerning only MSCs of different origins (including endometrium-derived stromal cells), POP4 and EIF2B1 show average stability (the fourth and fifth places according to the RefFinder stability ranking, respectively).

The most stable candidate reference genes for MSCs were 18S rRNA and B2M (Fig. 5). 18S ribosomal RNA, which is encoded by the RNA18SN5 gene, is the structural component of the eukaryotic ribosomal small subunit and is involved in the catalysis of protein synthesis at ribosomes [34]. The 18S ribosomal RNA gene is highly conserved among eukaryotes due to its function and is often considered to be stable [35, 36, 37]. However, 18S is an unsuitable reference in experiments where only mRNA is isolated from cells or only oligo(dT) primers are used during reverse transcription. This is because rRNA contains no poly(A) tail [38]. In our experiments, we used both oligo(dT) and random primers during cDNA synthesis, so we suggest that 18S rRNA is a suitable reference gene for MSCs.

The B2M was also found to be of low variability in MSCs (Fig. 5). In contrast, our data concerning cancer cells only and cancer cells+MSCs revealed the high variability of B2M between samples (Figs. 3,4) [39]. B2M encodes $\beta$2-microglobulin, an essential component of MHC class I molecules involved in immune modulation and surveillance. MHC class I complexes are present on the surface of most cells in the body, indicating their belonging to this organism. B2M has been previously shown to be an unstable reference gene not only within the malignant cells originating from different tissues but also within a set of lung cancer cell lines, where B2M demonstrated the lowest expression stability out of ten commonly used housekeeping genes [40]. Unstable expression of this gene in cancer cells may be explained by the fact that B2M expression often changes during malignisation [41]. Downregulation of MHC-I has been described in 40–90% of human tumors, including melanomas, squamous cell carcinoma and neuroblastoma, often correlating with worse prognosis because loss or mutations in MHC-I components allow malignant cells to escape from their recognition by cytotoxic CD8+T-cells [42]. However, in some cancer types, such as gliomas, B2M is overexpressed [43]. B2M product functions outside the cell but not intracellularly as the products of other tested reference genes do. Under in vitro conditions, there is no interaction with cells of the immune system, and accordingly, there is no natural selection for the ability of cancer cells to avoid this interaction. Combined with genomic dysregulation in malignant cells, this might lead to the high variability in vitro in B2M expression in cultured cells. Indeed, B2M is often viewed as an unsuitable reference in cancer cell studies [44].

We found the most variable gene for MSCs to be GAPDH (Fig. 5), while in cancer cell lines, GAPDH demonstrated high stability (Fig. 4). GAPDH is one of the most popular and widely used reference genes. The primary but not the only role of GAPDH in a cell is the enzymatic conversion of glyceraldehyde-3-phosphate to 1,3-biphosphoglycerate in the glycolysis cascade, a step required for glucose metabolism [45]. Barber et al. [46] conducted a comprehensive study concerning GAPDH stability in 72 human non-malignant tissues. They found that the expression of GAPDH was highly variable between tissues, which is consistent with our data. According to Barber et al. [46], the expression of GAPDH differed by a maximum of 15-fold between the highest and lowest expression in skeletal muscle and breast tissue types, respectively. However, in some tissues, such as the nervous system, no significant variation was observed [46].

In cancer cell lines, the most unsuitable gene for normalisation was HPRT1 (Fig. 4). The HPRT1 gene encodes hypoxanthine phosphoribosyltransferase 1, the enzyme involved in the synthesis of both guanine and inosine in a salvage pathway [47]. The unsuitability of HPRT1 as a normalisation control in cancer studies has been demonstrated. HPRT1 was found to be highly variable in cancer cells and malignant tissues originating from patients with lung, colon, prostate and pancreatic cancer at both the mRNA and protein levels [48]. However, in some experimental models, HPRT1 may be chosen as a suitable reference [49].

We did not find any relation between the expression level (observed Cqs, Fig. 2) and the reference gene stability rankings. For example, 18S was the most stable, and GAPDH was the least stable in the MSCs, although both of them were highly expressed. In the cancer cell lines, POP4 and GAPDH are stably expressed, but POP4 is one of the lowest-expressed genes. Thus, the expression level should not be considered when selecting reference genes for RT-qPCR.

If we compare housekeeping gene stability values obtained by all four algorithms in cancer versus non-malignant mesenchymal cells, the values are significantly higher in the cancer cell group. The differences in gene expression in cancer cells are difficult to predict, as in some cases, gene expression can be noticeably different even in the same type of cancer due to differences in the amount of transcripts of certain genes [4]. In contrast, the housekeeping gene expression of MSCs from different tissues was highly homogeneous (Fig. 5).

To visualise the effect of the normalisation strategy on the RT-qPCR results, we normalised the relative expression of the EGFR gene by the three most-stable or three least-stable reference genes in the entire panel of cell lines (Fig. 6). Using optimal normalisation strategy, we found that EGFR mRNA was increased in A431 and decreased in the MCF-7 and SK-UT1-B cell lines compared to HeLa. However, in the case of normalisation to three unstable reference genes (HPRT1, RPL13A, B2M), no differences between these cell lines in the EGFR mRNA content were detected. Instead, a false elevated EGFR gene expression was found in the SK-BR-3 cell line, cells that are known to express a low level of EGFR [50].

Our data are in accordance with previous published data on EGFR protein expression in these cell lines [51, 52]. Epidermal growth factor receptor, encoded by the EGFR gene, is one of the key growth factor receptors, involved in the control of cell growth, maintenance, proliferation and migration. EGFR disturbances, such as mutations and EGFR overexpression, are often observed in cancer cells and appear to be one of the reasons for cancer progression. One of the well-known examples of cancer cell lines with overexpressed EGFR is the human epidermoid carcinoma cell line A431. A431 cells are characterized by the amplification (up to 110 times) of the EGFR gene, elevated mRNA and protein level. EGFR protein levels have been found to be 2–100 times elevated above that present in normal fibroblasts [53]. Our data reflected these alterations and large differences in EGFR expression within cancer cell panel. Moreover, our data are in accordance with previously obtained results concerning EGFR protein content in MSCs, which is similar within different MSCs and comparable to EGFR levels in HeLa cells [52]. Thus, using geometric averaging of the expression data for several reference genes without checking the stability of the wide panel of candidate genes could lead to a misinterpretation of RT-qPCR results.

It should be noticed that our data are limited to analysed cell lines, as reference gene expression might differ in another cell types or experimental conditions. According to different estimations, about 3000–7000 genes have been defined as housekeeping genes (i.e., those that are ubiquitously expressed in all tissue/cell types; involved in maintenance of basal cellular functions) [54]. At present, it is obvious that no universal reference genes exist for RT-qPCR analysis. However, the accumulated information from many studies, including this one, may narrow the list of potential candidate genes and ease the choice of suitable reference genes for certain experimental tasks, conditions and biological materials (cells/tissues) used. Nevertheless, the reference gene stability should be validated for each experimental model.

Our results emphasise that different reference genes should be used in RT-qPCR studies depending on the cell types. We showed that the POP4 and EIF2B1 genes are the most suitable among housekeeping genes tested for the entire panel of cells, while POP4 and GAPDH could be the reference genes of the best choice for cancer cells, and 18S and B2M for non-malignant MSCs. We showed that housekeeping gene expression in MSCs is more homogeneous than in cancer cell lines. Additionally, we validated our choice of reference genes by comparing EGFR gene expression normalised to the three most-stable (POP4, EIF2B1, YWHAZ) versus the three most-variable (HPRT1, RPL13A, B2M) housekeeping genes according to RefFinder ranking. The results of the present study could be helpful for the design of the experiments using the aforementioned cell lines. However, it should also be noted that even when the same cell lines are analysed under new experimental conditions, such as treatment with different pharmacological agents, hypoxia, oxidative stress, epigenetic modifications, etc., the optimal reference genes should be independently verified. Our research aims to facilitate the selection of a subset of suitable reference genes for RT-qPCR analysis from the wide panel of cell lines frequently used as pathological or normal models in various fundamental studies.

All data reported in this paper will be shared by the Vera Kosheverova upon reasonable request.

VK and AS designed the research study. VK, AS and RK performed the research. VK, AS, EK and MK analyzed the data. VK wrote the first draft of manuscript. EK reviewed and edited the first draft of the manuscript. EK managed the research project and acquired funding. All authors contributed to editorial changes in the manuscript. All authors read and approved the final manuscript. All authors have participated sufficiently in the work to take public responsibility for appropriate portions of the content and agreed to be accountable for all aspects of the work in ensuring that questions related to its accuracy or integrity.

Not applicable.

This research was funded by the Russian Science Foundation grant #23–14–00335.

The authors declare no conflict of interest.