The human breast adenocarcinoma cell line, MCF-7 has been a standard model among researchers for about five 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 sufficient stability and reproducibility [16].
Further, laboratories can argue that by adhering strictly to the recommendations from Good Cell Culture Practice (GCCP) [47] and employing SNP/STR cell authentication techniques, one can reproduce their results with MCF-7 sub-clones. However, this may not be necessarily enough in connection with MCF-7 cells. In fact, MCF-7 estrogen disruptor assay failed to get international validation in 2016 by US NICEATM (US National Toxicology Program Interagency Center for the Evaluation of Alternative Toxicological Methods). The main reason cited for failure was centered on concerns regarding the inter-laboratory reproducibility of results [48].
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 difficult 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 specific 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. 2.2 (calculated from 2-Cq linearized values), 2-ΔCq analysis (with calibrator) using ANOVA (or Kruskal Wallis test) seen in Sect. 2.3, to roughly estimate the stable reference genes. Further, they suggest using NormFinder (Sect. 2.4) to analyze the stability of the genes and conclude the best pair/s of reference genes (Sect. 2.9). Lastly, it was recommended following the selected pair/s further by analyzing the normalized profile of the target gene/s (gene of interest) to validate the selection (Sect. 2.10 and 2.11).
The Cq (quantification cycle) values (previously referred to as Ct – threshold cycle or Cp – crossing point) represent the number of cycles that are required for the fluorescent 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 simplified way, reflects 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 Table 1, one could argue that the mean Cq values of all genes in both cultures are very comparable to each other and hence the gene expression is stable for all genes in both cultures. However, equal raw Cq values doesn’t necessarily imply that the gene expression is also uniform, it merely implies presence of nearly equal template mRNA in the wells. A 2x-3x fold significant difference between gene expression may still exist once the raw Cq values are converted to a linear scale. As shown by CV% analysis (Table 2), the variance between gene expression in both cultures is evidently present, that could otherwise have been overlooked in case raw Cq values were only considered. Further, the 2-ΔCq analysis (Figs. 1 and 2) also shows significant fold change for many genes over multiple passages (Supplementary table S1, see Additional file 1). The raw Cq values can, however, give information regarding the intragroup and intergroup variation, as the gene with low variations both inter- and intra-group would have less Cq dispersion [51]. Also, the use of Cq values can be done on semi-quantitative basis whereby a comparison simply indicates that there were fewer/more copies of gene A than gene B in a well.
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 identified as potential genes for culture A1, while GAPDH-RNA28S were identified 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 amplified/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 filament 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 significantly 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 workflow 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 workflows 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 profiles of the reference genes.
Two different approaches were employed to analyze legacy dataset compare the gene expression, namely normalized_count (file extension rsem.genes.normalized_results) and TPM obtained from scaled_estimate (file 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 findings in the present study and bolsters their selection, GAPDH was revealed to unstably expressed, in contrast to our findings. 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 difficult 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 identified that can increase the assessment consistency of normalization. Indeed, CCSER2 was ranked highly in culture A1 in our analysis as well and was confirmed 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 purified mRNA samples and their high abundance compared to target mRNA transcripts making it difficult 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 verified as a candidate stable reference gene for breast cancer tissue and normal tissues [32, 79]. ACTB and HSPCB were the best reference genes identified 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] identified 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 significant 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.