Which Reference Genes to Choose for qPCR Normalization? A Comprehensive Analysis in MCF-7 Breast Cancer Cell Line

very little about how these are expressed in MCF-7 sub-clones and when the cell line is cultured over multiple passages (p), given the heterogenic expression behavior often associated with MCF-7 cell line. We investigated the expression dynamics of 12 previously reported and suggested endogenous reference genes using RT-qPCR, available algorithms (NormFinder, geNorm, BestKeeper etc.) and TCGA transcriptomic analysis within identically cultured two sub-clones (culture A1 and A2) of MCF-7 cell line cultured over multiple passages. Further candidate reference genes were used to 4 (2 simulated and 2 backed by qPCR data) to make an evidence-based recommendation of the least variable reference genes that could be used in MCF-7 cell line.

In 2013, an updated approach based on IHC (immuno-histochemistry) markers was introduced by St.
Gallen International Expert Consensus on Primary Therapy of Early Breast Cancer to determine subtypes of breast cancers [2]. Accordingly, luminal A subtype tumors were de ned as estrogen receptor positive (ER+), progesterone receptor (PR) ≥ 20%, HER2 negative, Ki67 proliferation marker < 14% and if available, low recurrence risk tumors based on gene-based assays [2]. MCF-7 neoplastic cells were found to be positive for both Estrogen (ER) and Progesterone (PR) receptors along with having low metastatic activity and hence ful lled the criteria to be classi ed as luminal A molecular subtype tumor cell line [3].
Cell based assays (cell lines) such as MCF-7 represents techniques that can provide more biologically meaningful information than simpli ed biochemical assays [4]. The key reasons for their universal adoption are lower operational costs and the ease of operation in terms of preparing and observing the cells. Further, they represent an unlimited self-replicating source that can be grown in almost in nite quantities [5] yielding unlimited amount of DNA/RNA that enables studies related to validation and downstream functional analysis.
However, MCF-7 cell line like other cell lines is also prone to certain disadvantages. It is vulnerable to genotypic and phenotypic drift during its long-term culturing [5] and is of profound concern since the cell line has been deposited in cell banks for many years now. This has risked and in some certain cases caused arising of subpopulations within the cell line. Subpopulations can cause phenotypic changes over time by the selection of speci c, more rapidly growing clones within a population as demonstrated by Osborne et al. [6] and Resnicoff et al. [7] in 1987.
Over the past decades, extensive evidence that MCF-7 cells showed clonal variations has been reported depicting either differences in phenotypic traits such as estrogen/progesterone responsiveness, epidermal growth factor (EGF) expression or the ability to form tumors in syngeneic mice [8]. Also, genetic variability in the sublines and sub-clones of MCF-7 cell line on karyotypic and chromosomal levels have also been demonstrated by various researchers [9,10,11,12,13,14,15]. Finally, in 2016 Andre et al., illustrated variations among MCF-7 cell line obtained from the same cell bank in the same batch [16].
Presence of such heterogeneity bolsters the need for validating and cross-examining the genetic variability as well as gene expression in the cell line.
A widely accepted technique for validation of gene expression is the Reverse Transcription -quantitative Polymerase Chain Reaction (RT-qPCR). It is highly sensitive, reproducible, simple and high throughput yielding technique that can con rm gene expression differences and measure transcript abundances [17, set or set references to correct for any sampling noise such as differences in the amount of starting material in order to estimate results accurately [19].
These reference genes (previously housekeeping genes) are expressed constitutively and are required for the maintenance of basal cellular functions. In general practice, it is presumed that the endogenous reference gene represents an ideal gene that is su ciently abundant and has stable expression across different tissues and cell lines under different experimental conditions [20]. In addition, it is assumed that the expression levels remain the same amongst biological replicate cell cultures over successive passages. However, some studies suggest that the expression of these reference genes may not be as uniform as previously thought and may also uctuate signi cantly under different experimental conditions [21,22]. Hence, it becomes imperative to validate reference genes before their use in a study as using a non-validated reference gene could lead to misleading interpretations arising from inaccurate results [20,23].
As pointed out in the MIQE guidelines [24], normalization against a single reference gene is not recommended unless a clear evidence of its uniform expression dynamics is described for the speci c experimental conditions. Many studies have been undertaken previously identifying new stable reference genes over the previous two decades for MCF-7 cell line or breast cancer as whole [25,26,27,28,29,30,31,32].
However, these studies didn't undertake a detailed analysis validating the reference gene expression over multiple passages and/or, within sub-clones (biological replicate cultures). This gap in validation leaves us with a cause of serious concern especially if such cell lines are to be regarded as valid models for evaluating the behavior and development of breast cancer and validating their likely response to new drugs and therapies.
The present study, therefore, aims to ll the void by investigating the gene expression of twelve reference genes that were previously identi ed as stable genes by various studies [25,26,27,28,33,34,35,36,37] for MCF-7 cell line but were not accounted and studied together so as to make an evidence supported recommendation of reference genes to be used for normalization of gene of interest in routinely cultured MCF-7 cell line. Further, the study investigated for any differential reference gene expression in two subclones (culture A1 and A2) cultured over multiple passages (p) to further bolster the selection of appropriate reference genes.

Curation of the Dataset and Descriptive Analysis
Three lysates were collected from each passage from both MCF-7 cultures A1 and A2 while two lysates were collected and evaluated for passage 28 (p28) and passage 31 (p31) from culture A1. A detailed selection criteria of candidate reference genes is described in Sect. 5.3. Ampli cation for each of the 12 reference genes produced a dataset with 900 Cq values. As shown in Table 1, RNA18S showed the highest expression in both cultures A1 and A2 (Cq mean = 7.93 and 8.26 respectively), closely followed by RNA28S (Cq mean = 8.27 and 8.35 respectively). Both genes showed high ampli cation levels, appearing close to seven cycles earlier than any other gene in both cultures (ACTB Cq mean = 15.98 and 16.11 respectively). CCSER2 presented the lowest expression levels in both cultures (Cq mean = 26.58 and 26.56 respectively). GAPDH showed the lowest standard deviation (S.D = 0.30) in culture A1 while ACTB showed the lowest standard deviation (S.D = 0.32) in culture A2. The largest variation between Cq values was shown by HNRNPL (S.D = 0.35) in culture A1, while in culture A2, PGK1 (S.D = 0.61) showed the largest difference.  [39]. It has been reported that for a heterogenous sample like breast cancer cells, a CV% lower than 50% is ideal for reference genes to be considered to have stable expression [40].  To determine signi cant relative expression changes, one-way ANOVA was used (Supplementary Table  S1, see Additional le 1). In culture A1, only CCSER2 (Fig. 2C) and PCBP1 (Fig. 2E) showed no signi cant expression changes over successive passages (ANOVA P > 0.05). In culture A2, all the genes selected, including CCSER2 and PCBP1 showed signi cant expression changes over successive passages. Both ACTB (Fig. 1A) and GAPDH (Fig. 1B) showed signi cant expression differences (ANOVA P < 0.01) in both cultures over successive passages, providing evidence that these should be avoided to be used together, if possible, as endogenous reference genes in MCF-7 cell lines.

NormFinder Analysis of MCF-7 Sub-clones
NormFinder [41] estimates the overall variation of gene expression for each candidate reference gene and delivers a stability value, not only identifying the most stable gene but also the best control gene as shown in Fig. 3. It presents the stability of a gene as an estimate of the combined intra-and intergroup variation of the individual gene. It means that the gene with high expression stability is shown by a low standard deviation value and vice-versa.
Further, the algorithm can be used to nd the pair of genes best suited to work together as reference genes as shown in Supplementary Table S2 (see Additional le 1) (only top 10 pairs for both the cultures with their pairwise standard deviations). The algorithm ranked ACTB and PUM1 as the most stable pair (S.D = 0.07) in culture A1 while in culture A2 the rst spot was shared by 2 pairs of genes -RPL13A -SF3A1 and GAPDH -SF3A1 (S.D = 0.07).

geNorm Analysis and Determination of Optimal Number of Genes Needed for Normalization of Dataset
The geNorm algorithm [20] calculates M value for each candidate reference gene based on pairwise comparison. The higher the M value, the less stable the reference gene and vice-versa (Fig. 4). geNorm uses the stepwise exclusion method of the least stable genes by calculating the average M values. It has been reported that for a heterogenous tissue, such as breast cancer, a stability M value of < 1 should be considered for candidature of reference genes. Any gene having stability value M above 1 should not be considered as candidate for reference gene [40]. There were noticeable differences in the results in both cultures. In culture A1, ACTB-HSPCB tied for the rst position with a stability value M = 0.169 (Fig. 4A) while in culture A2, RNA18S-RNA28S tied for the most stable gene (M = 0.177; Fig. 4B).
geNorm can also be used to determine the optimal number of reference genes needed for an accurate estimation of normalization of expression data as shown in Fig. 5. In principle, Vn/Vn + 1 should be less than 0.15, where Vn represents the number of reference genes suitable for normalization. Although it is suggested that 0.15 cutoff be used, it has been reported that this cutoff may not be used strictly to assess the Vn/Vn + 1, especially if Vn/Vn + 1 is 2/3. In such a case, use of three stable control genes should be enough to provide accurate results. For both the cultures A1 and A2, the V2/3 (0.00437 and 0.00468 respectively) was less than the recommended cutoff, indicating addition of a third reference gene would not make a difference in the normalization results.

BestKeeper Analysis
BestKeeper [42] was used to analyze the candidate reference genes expression stability. BestKeeper uses crossing points (CP) to decide, whether the genes are differentially expressed under the applied conditions or not. To arrive at a hierarchy of the best reference gene according to BestKeeper, different criteria have been suggested by Pfa [42], such as considering the i) standard deviation with crossing points (S.D ± CP; recommended cutoff ≤ 1), ii) standard deviation with changes in x-folds (S.D ± x-fold; recommended cutoff ≤ 2) and iii) coe cient of correlation (r; recommended to be as high as possible). Different authors use and report different criteria from the above-mentioned possibilities. For comparison in our study, standard deviation with crossing points was considered rst for both the cultures ( Fig. 6A and B).
As seen in other two previous methods (Sect.

Pairwise Comparative ΔCt Analysis of the MCF-7 Subclones
Comparative ΔCt [43] compares the relative expression of pairs of candidate reference genes within each sample in order to identify and rank the most stable genes. According to Silver et al., if the ΔCq value of the two genes uctuate when analyzed in different samples, it can be concluded that, one or both genes are variably expressed [43]. Conversely, if the ΔCq value remains constant, both genes are expressed stably or are co-regulated among the samples. For a stable gene, it has been suggested that the gene should be strongly expressed, displays minimal uctuations and is independent of expression of other genes. The ranking of the genes was based on the average Standard Deviation (S.D). The higher the S.D, the more variable the expression and less stable the gene and vice-versa (Table 3). In culture A1, GAPDH, CCSER2 and PCBP1 were the top three genes ranked by the algorithm while in culture A2, RNA28S, GAPDH and PCBP1 were the top three genes. It is interesting to note that in culture A1, the three top genes ranked by Comparative ΔCt (Table 3) are ranked in the same order with same rank by BestKeeper (Fig. 6A). In culture A2, only RNA28S has been reported in the top three by both these algorithms (Table 3; Fig. 6A).

RefFinder Analysis
RefFinder [44] is a user friendly, easy to use web based comprehensive tool (https://www.heartcure.com.au/re nder/) that uses NormFinder, geNorm, BestKeeper and Comparative ΔCt to rank and compare candidate reference genes. It measures the geometric mean of attributed weights by these algorithms for generating an overall nal ranking. The rankings from RefFinder are summarized in Table 4. The top three genes in culture A1 were reported to be GAPDH, CCSER2 and PCBP1. For culture A2, RNA28S, GAPDH and PCBP1 were the top three most stable genes. Coe cient of Variation (CV) is the only method where the stability of the gene is not in uenced by any other genes in the study and hence can help in identifying a gene with overall high passage variance [38]. As mentioned above in the CV section, any gene having CV% greater than 50% (threshold) should be excluded from the study. Since none of our genes had CV% greater than the threshold, the next step was to evaluate the results from the ve different approaches as mentioned in Sect. 2.3 to Sect. 2.8.
Based on these approaches, we came to a conclusion that for culture A1, GAPDH and CCSER2 were suitable reference genes as they both were constantly ranked in top three in NormFinder (Fig. 3A), BestKeeper ( Fig. 6A), Comparative ΔCt (Table 3) and RefFinder (Table 4). Further, CCSER2 had no signi cant relative mean expression change over successive passages (Fig. 2C) while GAPDH had shown the least CV% (Table 2) in culture A1. Also, the Pearson Correlation between GAPDH-CCSER2 (r = 0.515) was signi cant (P = 0.004) thereby indicating their positive correlation (Fig. 7A) and corroborating their selection. Based on a similar approach, for culture A2, GAPDH and RNA28S were selected as suitable reference genes. Interestingly, to not exclude a potential candidate for reference gene just because of strict selection of only 2 reference genes as suggested by geNorm (Fig. 5), an experimental candidate with constant high rankings in all the algorithms for both the cultures was also chosen. PCBP1 was identi ed as a strong contender in the race and hence was included in the validation of selected reference genes (Sect. 2.10 to Sect. 2.11). PCBP1 was one of the two genes to have no signi cant relative expression changes in culture A1 (Fig. 2E) and claims the top spots in NormFinder (Fig. 3), BestKeeper ( Fig. 6), Comparative ΔCt (Table 3) and RefFinder (Table 4). It also maintains a constant position in top ve genes per geNorm analysis for both cultures ( Fig. 4A and 4B).

Generation of Data for 2 Genes of Interest (GOIs) & Validation of Reference Genes Using Normalization
To further evaluate the four selected reference genes (GAPDH, RNA28S, CCSER2 and PCBP1) and their utility as reference genes, two simulated datasets were created. The genes simulated were referred to as Gene of Interests (GOI) 1 and 2. Both were assigned some random Cq values (triplicates for 3 lysates over 4 passages for culture A1 and over 5 passages for culture A2) as presented in Supplementary table S3 and S4 (see Additional le 1).
While assigning the Cq values for GOI 1 (simulated to have stable gene expression), it was considered that the difference between Cq values doesn't exceed +/-0.5 Cq. To randomize the Cq values and minimize human intervention, the values were generated using an online random fraction generator tool (https://onlinerandomtool.com/generate-random-fractions). The GOI 1 was then normalized using ΔΔCt method [45] to the selected reference genes for both cultures. Signi cance of fold change was estimated using P < 0.05 (indicating signi cant change) from one-way ANOVA and respective post-hoc tests (with Bonferroni P value correction).
For culture A1, as suggested by different algorithms used in this study, GAPDH -CCSER2 pair did prove its utility when it was used for normalization of the simulated GOI 1 (Supplementary Figure  On similar lines GOI 2 (simulated to have unstable gene expression) was simulated and assigned random Cq values (Supplementary Table S4, see Additional le 1). This time it was considered that the difference between some Cq values should exceed +/-0.5 Cq. The GOI 2 was also normalized using the ΔΔCt method [45] to the selected reference genes for both cultures. For culture A1, GAPDH -CCSER2 pair produced statistically signi cant results (ANOVA P = 0.035) in fold change after normalization of GOI 2, indicating the pair may not after all be able to handle large differences in Cq values (greater than +/ 0.5 Cq) as seen in Fig. 8A and Supplementary Figure S2A (see Additional le 1). Interestingly, all the other reference gene pairs proved their utility and produced insigni cant changes after normalization of GOI 2 (ANOVA P > 0.05) in culture A1 ( Fig. 8A and Supplementary Figure S2A, see Additional le 1). Here as well, for culture A2 GAPDH-RNA28S pair produced statistically signi cant fold changes (ANOVA P = 0.07), indicating it to not be a suitable pair for normalization. All pairs in culture A2 showed only statistically signi cant changes in one passage p25/p29 as seen in Fig. 8B (also see Supplementary Figure S2B, Additional le 1), and hence further normalization with 3 reference genes was done to look for successful normalization.
After performing further normalization of GOI 2 for cultures A1 and A2 in pairs of three reference genes (Supplementary Table S6, see Additional le 1), our analysis showed no suitable pair of reference genes for GOI 2 in culture A2 ( Fig. 8D) but for culture A1, GAPDH-PCBP1-CCSER2 was the only triplet to yield successful normalization (Fig. 8C).

Validation of Selected Reference Genes with Normalization of qPCR Data
To support and continue perusing the selected reference gene pairs, 2 genes of interest namely AURKA (Aurora Kinase A) and KRT19 (Keratin 19) were also used. The qPCR was done, and the dataset was normalized using the reference gene pairs as shown in Supplementary Figure S3 Upon normalization of AURKA, none of the gene pairs in culture A1 were able to normalize AURKA adequately (P < 0.05) as seen in  clones A1 and A2 were combined and the analysis was done. CV% was determined to evaluate the stability of reference genes ( Table 5). As before, comprehensive analysis was done for the combined dataset using NormFinder, geNorm, BestKeeper, Comparative ΔCt and RefFinder. The gene rankings were determined using each of these algorithms as summarized in Table 5. sub-types were removed). In the dataset, selected reference genes were looked for. Normalised Gene expression data was available for 9 of the 12 reference genes (data not available for ACTB, RNA18S and RNA28S) as presented in Table 6. The "normalised_count" is a simple transformation of the "raw_count" (available in the le with extension rsem.genes.results). The "raw_count" values are divided by the 75th percentile and then multiplied by 1000 to obtain "normalised_count" to make the values comparable between experiments. For better data visualization and understanding gene expression of the reference genes, the "normalized_count" were converted to the log scale using log2 (normalized_count + 1) as presented in Supplementary Figure S5. CCSER2 is the most stable gene (from our selected reference genes) identi ed in the database for luminal A sub-type Breast cancer. Further, PCBP1 is also quite stably expressed (Supplementary Figure S5 and Table 6). Interestingly, GAPDH which was identi ed as the most stable gene in both cultures was in contrast identi ed the least stable gene in the transcriptomic analysis.
Going ahead, to evaluate the TPM (transcripts per million), the TCGAbiolinks was used again to retrieve the clinical, morphological and expression data of luminal A sub-type breast cancer patients under the project "TCGA-BRCA" (Platform -Illumina HiSeq; le.type -rsem.genes.results). The "scaled_estimate" (estimated frequency of gene/transcripts among the total number of transcripts that were sequenced) values were obtained. TPM is obtained by multiplying "scaled_estimate" values by a factor of 1 million (10 6 ). In case of TPM, the data was available for 10 of the 12 selected candidate reference genes (except for RNA28S and RNA18S) as shown in Table 7. For better visualization and data analysis, the TPM values were also converted to logarithmic scale using log2(TPM).

Discussion
The human breast adenocarcinoma cell line, MCF-7 has been a standard model among researchers for about ve decades now, serving as a laboratory tool for in vitro studies as well as model for investigation of key cancer driven processes that directly impact patient care and treatment plan [46]. Despite publication of extensive evidence [9,10,11,12,13,14,15], the genetic and phenotypic variance in MCF-7 sub-clones and subpopulations is often not accounted for in the laboratory protocols and is hence overlooked. The main reason for this overlook usually stems from assumptions that by using cells obtained from same batch and same cell bank and by standardizing protocols and limiting the number of passages, laboratories can ensure that their sub-clones will "behave" with su cient stability and reproducibility [16]. Previously as well, the question about accuracy in reproducibility of results with MCF-7 cells has been raised [49], but no substantial proof was reported that dealt with variations in reference gene expression tendencies. The present study addresses this issue and reports a novel evidence that non-tested and vague use of common reference genes among sub-clones for qPCR normalization may lead to inaccurate results and conclusions. A pivotal aspect in any gene expression related study is the selection of the most appropriate and accurate reference gene/s. Given the widespread availability of tens of tools and algorithms including web-based platforms (RefFinder) and stand-alone or R-based algorithms (geNorm, BestKeeper, NormFinder etc.), it becomes di cult for researchers without in-depth mathematical background to choose the appropriate software based on their experimental needs.
Each of these algorithms in the past have been reported to have some advantages and some drawbacks that could render the results based on the speci c experimental needs. Considering the limitations and advantages of all the algorithms, Sundaram et al [39] in their study, suggested to use an approach that involves using CV% analysis seen in Sect. represent the number of cycles that are required for the uorescent signal to cross the threshold. The raw values obtained from and used by the PCR machine aren't necessarily suggestive of the absolute quantity of template in terms of intergenic/inter-run comparisons, but in a simpli ed way, re ects the abundance of the individual template mRNA in wells [50]. As Livak et al., described, the raw Cq values are determined from a log linear plot of PCR signal vs the cycle number, thereby making Cq an exponential term rather than a linear term [45]. Ergo, for any gene expression analysis or comparisons it is necessary to convert the Cq values to a linear scale using the 2 -Cq method [45].
It is interesting to point that on an initial analysis of Our analysis revealed that among the two biological replicate cultures (sub-clones), there are stark differences in the expression patterns and tendencies of the endogenous reference genes. GAPDH-CCSER2 were identi ed as potential genes for culture A1, while GAPDH-RNA28S were identi ed for culture A2 using various algorithms. However, when they were employed for cross-normalization of genes of interests in both cultures, both gene pairs were unable to provide adequate results (Fig. 8). Addition of a third gene PCBP1 to GAPDH-CCSER2 pair helped to yield successful normalization for all 4 genes of interest and in both cultures A1 and A2 (except for GOI 2 in Culture A2; this can be attributed most likely to limitation of GOI 2 as being a simulated gene).
AURKA (Aurora Kinase A) has been known to be associated with playing a key role in centrosome duplication and chromosome segregation during mitosis [52]. Further, it has been reported at many instances to be ampli ed/mutated in several human cancers [53,54,55,56,57,58]. In breast cancer, mixed evidence has been reported with some reporting its association with basal like phenotype [59,60] whilst others suggesting it to be a marker for progression and outcome of luminal like subtype [61].
Keeping this in mind AURKA made a perfect candidate as gene of interest. Although when reference genes were taken in pairs of two, only GAPDH-CCSER2 normalized AURKA expression in culture A2, it was adequately normalized by GAPDH-CCSER2-PCBP1 triplet in both the cultures A1 and A2.
The other gene of interest selected was KRT19 ( Keratin 19). It belongs to the KRT (keratin) family which serves as important markers in RT-qPCR mediated detection of tumors in lymph nodes, peripheral blood and bone marrows of breast cancer patients [62]. KRT19 is one of the smallest intermediate lament KRT protein [63] and has been shown to regulate breast cancer properties [64]. Using Oncomine database and RT-PCR, Saha et al. [65] evaluated expression of KRT19 in MCF-7 and other breast cancer cell lines. They reported that KRT19 was signi cantly overexpressed in MCF-7, MDA-MB-231 and SKBR3, ergo validating its choice in the present study. None of the reference gene pairs in two cultures could normalize KRT19 when they were employed in groups of 2. However, in both the cultures, only GAPDH-PCBP1-CCSER2 yielded successful normalization thereby proving the ability of the triplet pair to handle genes of interest.
The TCGA (The Cancer Genome Atlas Program) represents a joint venture of National Cancer Institute (NCI) and National Human Genome Research Institute (NGGRI) which began in 2006 as a pilot project with three cancer types (lung, ovarian and glioblastoma) which got expanded to present 33 tumor types encompassing a comprehensive dataset describing the molecular changes that occur in cancer [66]. Most samples in TCGA were originally aligned against the Genome Reference Consortium build GRCh37 (hg19) or the "legacy" dataset [67]. However, with advances in technology and drop in sequencing costs GDC (Genomic Data Commons -conceived by NCI) undertook harmonization effort to align the data to GRCh38 (hg38) build ("harmonized" dataset).
The work ow for generating RNA-Seq data in both legacy and harmonized dataset differs substantially [67] leading to introduction of bias between the hg19 and hg38 abundance estimates. However, Gao et al., [67] demonstrated that there exists excellent concordance between the two work ows in relation to BRCA PAM50 subtypes. Further, they reported that relative change between conditions is preserved across all subtypes of BRCA PAM50. Hence, in the present study, the legacy dataset was accessed and analyzed for comparing the expression pro les of the reference genes.
Two different approaches were employed to analyze legacy dataset compare the gene expression, namely normalized_count ( le extension rsem.genes.normalized_results) and TPM obtained from scaled_estimate ( le extension rsem.genes.results). Both approaches revealed CCSER2 to be the most stably expressed. The CV% was estimated at 17.49% (Table 7) from TCGA analysis which was between the range of 15.55% (culture A1) and 25.65% (culture A2) obtained from RT-qPCR (Table 2) in present study. PCBP1 was also expressed quite stably with CV% estimated at 3.96% from TCGA analysis ( Table 7). The RT-qPCR range for PCBP1 obtained was from 14.89% (culture A1) and 24.90% (culture A2). Finally, CV% for the third gene in the triplet (GAPDH-PCBP1-CCSER2), TCGA estimates for GAPDH were at 6.44% while RT-qPCR range was from 12.35% (culture A1) and 26.84% (culture A2) ( Table 2).
Whilst TCGA analysis of CCSER2 and PCBP1 supported our ndings in the present study and bolsters their selection, GAPDH was revealed to unstably expressed, in contrast to our ndings. It is important to point out that the difference in results could be attributed to various underlying limitations. Firstly, TCGA requires all malignancies in its database to be primary, untreated tumors (except cutaneous melanoma) [68]. Further, all specimens deposited were garnered from available frozen materials from different institutes thereby introducing bias in institutional biorepository collections, stemming from institutional research interests, operative protocols, or patient populations [68]. In addition, metastatic diseases or aggressive primary tumors are usually subject to neoadjuvant therapies, which make their inclusion in TCGA database di cult because of limited availability of untreated specimens [68]. However, all these limitations don't contradict the fact that TCGA remains by far the richest source of clinical and research importance especially in developing an integrated picture of commonalities, differences and emergent themes across tumors.
Despite differences in TCGA analysis and our results, use of GAPDH as reference gene in qPCR normalization presents a controversy of its own. GAPDH has been shown to have increased expression in cancers from other body regions specially from cervix, prostate, pancreas, and lungs [69,70,71,72].
Furthermore, it has been reported that GAPDH is overexpressed in MCF-7 cells treated with estradiol [73]. Hence, many studies suggest not to use GAPDH as a control gene to study breast cancer or ranks GAPDH as the least stable gene [25,32,33,73,74,75]. As reported by Liu et al. [25], almost half of the publications in PubMed database used GAPDH as a single reference gene for normalization of gene expression analyses with RT-qPCR. Even with the contradiction regarding consideration of GAPDH as a suitable or non-suitable candidate, attention should be paid to its selection based on the experimental conditions and study design [33].
CCSER2 is described as a "novel housekeeping gene" (nHKG) by Tilli et al. [27], who provided evidence as to its use as a reference gene in breast cancer studies. They demonstrated that expression of CCSER2 is expected to be like the other nHKGs they identi ed that can increase the assessment consistency of normalization. Indeed, CCSER2 was ranked highly in culture A1 in our analysis as well and was con rmed by transcriptomic validation to be least variable in expression. It can be hence, used for future analysis and experiments. In addition, PCBP1, the most trustworthy gene in our analysis with a stable ranking across all platforms shows promising candidature as an endogenous reference gene as reported also by Jo et al [28].
Another gene that came to light in the analysis for culture A2 was RNA28S. The gene has not been reported much in the context for MCF-7 cell line and was added as an in-house suggestion in the present study. However, drawbacks to the use of RNA18S or RNA28S as reference genes have been reported such as absence of puri ed mRNA samples and their high abundance compared to target mRNA transcripts making it di cult to accurately subtract the baseline value in RT-PCR data analysis [20]. Further, the use of these genes as control genes is not suggested due to the imbalance between mRNA and rRNA fractions in these molecules [76]. In addition, it has also been shown that certain drugs and biological factors may also affect the rRNA transcription [77,78].
ACTB, another widely used reference gene, has been previously also veri ed as a candidate stable reference gene for breast cancer tissue and normal tissues [32,79]. ACTB and HSPCB were the best reference genes identi ed by Liu et al. [25] for ER + breast cancer cell lines including MCF-7. They further reported that RNA18S and ACTB was the best pair of genes across all breast cancer cell lines [25].
Despite the widespread acceptance of ACTB for normalization among a set of human breast cancer cell lines of increasing metastatic potential, limitations have been reported as well [29]. Jacob et al. [33] identi ed HSPCB as one of the suitable genes across a variety of cancer cell lines including MCF-7. Our analysis revealed contrasting results. ACTB, HSPCB, RPL13A and RNA18S were not ranked in top three in both cultures. Further, ACTB and HSPCB reported high CV% and showed signi cant fold changes (2 −ΔCq analysis). The difference in our results from the previously reported studies is suggestive of the fact that inter-laboratory replication of results with MCF-7 cell line is not very accurate.
Nonetheless, the authors caution the readers that the results obtained in the present study must be seen in the light of some methodological limitations. An interesting observation was made regarding passage p28 from culture A2 (Figs. 1 and 2), where low Cq values were obtained for all housekeeping genes when compared to other passages in the same culture. A plausible explanation for this observation cannot be provided at this point since Good cell culturing practices were strictly adhered to and the cell culturing, RNA isolation, quality control using PCR and RT-PCR were all done by the same single operator to minimize technical errors. Further, the RNA isolation for all passages (and all lysates) for both the cultures was done in one batch together on the same day and the same is also true for RT-PCR, thereby effectively minimizing any chance of human error.
The study with its results aims to provide the readers and researchers an evidence-based recommendation of the most suitable reference genes in routinely cultured MCF-7 cell line with extensive research and investigations. Further, we report the possibility of existence of a heterogenous differential behavior of the endogenous reference genes within the sub-clones of the MCF-7 cell lines. The authors suggest that more detailed and diverse studies should be undertaken to explore more about the differential expression of the endogenous reference genes. Future studies could be aimed at investigating the reference gene expression amongst other sub-populations of MCF-7 and amongst other breast cancer cell lines.
Lastly, given the widespread use of cancer cell lines not only for basic research but also for drug development and regulatory decision making, ensuring that sub-clones among the cell line are adequately standardized in their expression behavior tends to represent a challenging path forward [16]. A strategy outlined by Tilli et al. [27], whereby experiments are normalized with a panel of reference genes whose expression has been proven to be as minimally variable as possible and as robust as possible under varying conditions gains our support as well. Although geNorm recommended to use two reference genes, we would recommend that three reference genes may be employed to better overcome and handle any reference gene expression variability in the samples that may be present in the sub-clones. Triplet of GAPDH-PCBP1-CCSER2 should provide a potential alternative to traditionally used reference genes for reference gene matrix in MCF-7 cells.

Conclusions
Genetic and phenotypic heterogeneity showcased by MCF-7 cell line poses conundrum with some unanswered questions. Performing reference gene analysis may not be feasible for every sub-clone and hence, based on our results, we suggest using GAPDH-PCBP1-CCSER2 triplet on cultures that have been cultured in conditions that mimics that of our study, whilst seeing the results in the light of the limitations reported in the present study. Going forward, we encourage that studies should be undertaken especially those that can validate MCF-7 behavior in different growth conditions. We are optimistic that MCF-7 cell line will continue to contribute eminently towards improved and novel treatment modalities for breast cancer patients.

Culture and Seeding Conditions of MCF-7 Biological Replicate Cultures (Sub-clones)
Samples were collected from MCF-7 cell line (ATCC, HTB-22) that has been used in our laboratory for previous studies. Samples from different passages were cryopreserved during culturing at different time points. Frozen stocks were taken from MCF-7 cells grown in culturing conditions as described below. For this study, two aliquots from samples that were previously cryopreserved were thawed at the same time.
Culture A1 was cryopreserved after passage 27 (p27) and culture A2 was cryopreserved after passage 24 (p24). They were then cultured in the same conditions simultaneously over multiple successive passages.

Selection of Candidate Reference Gene and Primer Design
In total, 12 candidate reference genes were selected by perusing relevant literature related to breast cancer and/or MCF-7 cell line and literature that reported transcriptomic data based on TCGA (The Cancer Genome Atlas) database [25,26,27,28,33,34,35,36,37]. All selected genes are shown in Table 9. GAPDH and ACTB are two of the most used single control genes reported in more than 90% of the cases in high impact journals [80]. Further these two genes are commonly included in the commercially available cancer pathway kits like those from Qiagen, Life technologies and are also included in Oncotype DX test arrays. Life technologies kit also includes genes like RNA18S and PGK1 along with several others [27]. RNA18S-ACTB has also been identi ed as an appropriate gene pair across all breast cancer cell lines [25]. RPL13A was described by De Jonge et al., as a novel candidate reference gene with enhanced stability among a magnitude of different cell types across varying experimental conditions [37]. HSPCB has been reported previously to be one of the most stable reference gene for ER + breast cancer cells [25,33], thereby in uencing its selection in the present study.
PUM 1 and CCSER2 were identi ed by Tilli et al., using transcriptomic analysis across breast cancer cell lines as less variable and more accurate for research in breast cancer cell lines and tissue samples in comparison to the traditional housekeeping genes [27]. Previously, HNRNPL and PCBP1 have been reported to be among the most highly expressed genes across breast cancer samples and performed better than both ACTB and GAPDH as supported by TCGA transcriptomic validation reported by Jo et al. [28]. ATCB and SF3A1 were identi ed by using high throughput analysis of micro-assay datasets with subsequent validation by RT-qPCR as reported by Maltseva et al. [36]. They further reported these genes as more e cient for analysis of breast cancer samples when compared with the reference gene panel provided by Oncotype DX assay. Finally. RNA28S was added to the selected genes as an experimental inhouse suggestion.
The genes selected, as described above, have been time and again reported as the best or have been described with controversial ndings with studies spanning over almost two decades from 2000 to 2019 and hence, they were never collectively investigated in MCF-7 cell line. The present study, hence, took these ndings from previous studies into account while selecting the genes. The sources of gene primers and primer sequence used for the reference genes are shown in Supplementary Table S9 (see Additional   le 1). The present study employed the primers that were reported before in the literature to maintain inter-reliability and inter-connectivity with the previous studies. Primers for RNA28S and PGK1 were designed using Primer3Plus [81]. The melting curves of all the selected gene primers are presented in Supplementary Figures S8 to S19 (see Additional le 1). Finally, the primer sequences (designed using Primer3Plus) for genes of interest, AURKA and KRT19 are also presented in Supplementary Table S9 (see  Additional le        Determination of optimal number of reference genes needed for normalization (geNorm) for both cultures Legend: Pairwise Vn/n+1 analysis done using geNorm. The recommended cutoff is the lowest Vn/n+1 below the threshold of 0.15. In our case, both the cultures had V2/3 less than 0.15, indicating addition of a third reference gene would not affect normalization results.