Selection of endogenous control genes for normalising gene expression data derived from formalin-fixed paraffin-embedded tumour tissue

Quantitative real time polymerase chain reaction (qPCR) data are normalised using endogenous control genes. We aimed to: (1) demonstrate a pathway to identify endogenous control genes for qPCR analysis of formalin-fixed paraffin-embedded (FFPE) tissue using bladder cancer as an exemplar; and (2) examine the influence of probe length and sample age on PCR amplification and co-expression of candidate genes on apparent expression stability. RNA was extracted from prospective and retrospective samples and subject to qPCR using TaqMan human endogenous control arrays or single tube assays. Gene stability ranking was assessed using coefficient of variation (CoV), GeNorm and NormFinder. Co-expressed genes were identified from The Cancer Genome Atlas (TCGA) using the on-line gene regression analysis tool GRACE. Cycle threshold (Ct) values were lower for prospective (19.49 ± 2.53) vs retrospective (23.8 ± 3.32) tissues (p < 0.001) and shorter vs longer probes. Co-expressed genes ranked as the most stable genes in the TCGA cohort by GeNorm when analysed together but ranked lower when analysed individually omitting co-expressed genes indicating bias. Stability values were < 1.5 for the 20 candidate genes in the prospective cohort. As they consistently ranked in the top ten by CoV, GeNorm and Normfinder, UBC, RPLP0, HMBS, GUSB, and TBP are the most suitable endogenous control genes for bladder cancer qPCR.

Bladder cancer is a major cause of morbidity and mortality in the UK. Data from Cancer Research UK showed that it is the tenth most common cancer in the UK, accounting for 3% of all new cancer cases 1 . There are multiple therapeutic options for bladder cancer, which highlights the importance of developing biomarkers for personalising treatment 2 . Emerging transcriptomic signatures can be progressed towards clinical application using different platforms including quantitative real time polymerase chain reaction (qPCR).
qPCR is a sensitive, affordable and efficient method for measuring gene expression in tissue samples including formalin fixed paraffin embedded tissue (FFPE). RNA from FFPE is generally of poor quality 3 . Formalin fixation results in cross-linking of RNA with other macromolecules including DNA and protein which when dissociated during RNA purification results in fragmentation and reduction in yield of probable material. Pre-amplification of cDNA from these samples is necessary to obtain quantifiable data 4 . Samples of cDNA can then be subject to qPCR. Most qPCR data measure relative gene expression via normalization with endogenous control genes (also known as reference or housekeeping genes). Genes used widely in the past can be affected by tissue type and experimental conditions [5][6][7] , and it is important to identify genes with constitutive and invariant expression for the samples of interest.
Some studies used multiple endogenous control genes that include several co-expressed genes including the ribosomal protein large (RPL) family of proteins 8,9 . Where co-expressed genes are present within a candidate gene Scientific Reports | (2020) 10:17258 | https://doi.org/10.1038/s41598-020-74380-7 www.nature.com/scientificreports/ panel their mutual influence on apparent stability requires consideration. Genes with similar functions tend to exhibit similar gene expression patterns 10 . Gene co-expression can be checked using the on-line tool Genomic Regression Analysis of Coordinated Expression (GRACE) which correlates (Spearman) the expression of a gene with all other genes within TCGA. Vandesompele et al. 11 propose testing a panel of candidate reference genes on a representative number of relevant samples to identify the most stable and optimal number of genes. Test data generated are subject to stability assessment algorithms, the two most commonly used are GeNorm 11 and NormFinder 12 . These algorithms rank genes in order of stability and in the case of GeNorm select the two-gene combination that provides the most stable normalization. GeNorm is considered the optimal algorithm for studies with small sample numbers 12 but over-rates the stability of co-expressed genes in the candidate panel.
Gene expression data are highly dependent on platform 13 so endogenous control gene selection is carried out on the platform of choice. To facilitate selection of control genes, TaqMan endogenous control array cards are available with 16 candidate genes. These genes have been used for normalization in human tissue gene expression studies including bladder 14 , thyroid 15 , hepatocellular 16 , breast 9 , gastric 17 , cervical 18 , endometrial 19 , non-small cell lung 20 , kidney 21 and prostate 22 cancer. Whilst bioinformatic interrogation of TCGA provides a useful verification of gene expression stability it is unsuitable for endogenous control gene selection for TaqMan Array cards as the TCGA database is derived using RNA sequencing.
The primary aim of this work is to facilitate selection of endogenous control genes for the Taqman qPCR gene expression platform for studies of prospective FFPE cancer tissue using bladder cancer as an exemplar. The secondary aim is to evaluate the influence of some covariables including probe length on reverse transcription efficiency and co-expression on stability ranking.

Materials and methods
Patient samples. Pre-treatment FFPE grade 3 MIBC samples were available from a prospective (n = 12) and retrospective (n = 16) patient cohort. Samples were obtained via the Manchester Cancer Research Centre Biobank under research tissue bank ethics (Ref: 18/NW/0092). The cases were graded by an experienced subspecialist Uropathologist (HD). Mean (range) block age was 6 (3-8) months for the prospective cohort and 15 (14)(15)(16)(17) years for the retrospective cohort. Two 10 µm sections for RNA extraction and a 4 µm section for histological verification of tumour cellularity (> 30%) were obtained from each block. RNA was extracted from the two 10 µm sections using the Roche High Pure FFPET RNA isolation kit. TCGA bladder cancer RNA-seq data (n = 401) were also analysed.
TLDA cards and single tube assays. Table 1 lists the endogenous control genes tested along with the probes and their amplicon length. Sixteen genes were on the endogenous control TLDA cards and single tube assays were set up for succinate dehydrogenase complex flavoprotein subunit A (SDHA) a gene demonstrating particularly low variability in bladder cancer cells 9 . Single tube assays were also run for RPL11, RPL24 and GNB2L1 gene to examine the biasing effect of co-expression on gene stability. To investigate the effect of probe size on Ct values two different probes were selected for RPL11, RPL24 and GNB2L1. RNA extraction, quantification and reverse transcription. RNA was extracted using the Roche High Pure FFPET RNA isolation kit from two 10 µm sections. The detailed steps of extraction were performed according to manufacturer's recommendations. RNA quantification and purity were determined on a NanoDrop UV-Vis Spectrophotometer (Thermo Fisher Scientific Poole UK) and a Qubit fluorometer (Invitrogen Paisley UK). Reverse transcription and pre-amplification steps were carried out on a 2720 thermal cycler (Applied Biosystems UK). qPCR was carried out on a Quantstudio 12K (Applied Biosystems UK). Complementary DNA (cDNA) was generated using a high capacity RNA-to-cDNA kit (Life Technologies Ltd UK). One sample of cDNA was subject to pre-amplification using a custom preamp pool mix consisting of gene expression assay corresponding with genes present on the TaqMan human endogenous control card array (Applied Biosystems®). A further sample of cDNA underwent pre-amplification using a preamp pool mix prepared by mixing single tube assay (Thermo Fisher Scientific UK) components for the panel of candidate genes not present on the endogenous control card array. A preamp TaqMan Fast Advanced Master Mix (Thermo Fisher Scientific UK) was used for both samples. The pre-amplification step was carried out for 14 cycles on a 2720 thermal cycler.
qPCR. Samples pre-amplified using the control array primer pool were mixed with Fast Mastermix (2X) and loaded onto the endogenous control plate and subject to qPCR on a Quantstudio12 (Applied Biosystems). Samples pre-amplified using the pooled single assay primer pool were loaded into 96 well plates preloaded with individual gene probes and Fast Mastermix (2x) and subject to qPCR on a Quantstudio12 using the fast cycle option. Data analysis. TCGA was accessed using the Firebrowse portal selecting RSEM normalised RNAseq bladder cancer. Each gene was examined for co-expression with other genes within the candidate panel using the on-line tool GRACE. Coefficient of variation (CoV) values for the expression of each gene were calculated and significant differences between mean values determined using the Student t-test.
GeNorm (https ://genor m.cmgg.be/) and Normfinder (NormFinder software-moma.dk) algorithms were used to determine the most stable genes from the list of candidates. The software packages were used as excel addons. To determine the biasing effect of co-expressed genes on apparent gene expression stability using GeNorm, the analyses were carried out with all candidate genes and repeated after removing 3 of the 4 co-expressed genes.
Statistical significance between means was determined using the Mann Whitney U test calculator (Wilcoxon rank-sum) non-parametric test. . Gene expression was determined in triplicate and the intra-sample standard deviation (SD) of the cycle threshold (Ct) values (number of cycles required for the fluorescent signal to cross a threshold) calculated for the 28 samples. Table 1 lists the mean SD for each of the 23 candidate endogenous control gene probes. The mean ± SD of the SD of the Ct values for triplicate assays for 16 endogenous control genes assayed in 28 samples on TaqMan arrays was 0.153 ± 0.071 (n = 448 gene-probe/sample combinations; range 0.079 to 0.37). For the 7 single tube assays the mean SD for the Ct values were 0.082 ± 0.022 (n = 196 gene-probe/sample combinations; range 0.046 to 0.129). To investigate inter-assay reliability three samples were assayed by qPCR on two separate TaqMan arrays/96 well plates run on two different days. The mean ± SD of the SD of the Ct values for each gene/sample (n = 48) was 0.21 ± 0.04 for samples loaded onto TLDA cards and 0.07 ± 0.039 for the single tube assays loaded into 96 well plates. Table 2 lists the mean, SD and CoV of the Ct values for each of the 23 (16 on the TLDA array card and 7 single assay) gene probes assayed in the prospective and retrospective patient cohorts. The mean Ct values for all 23 probes were significantly lower when assayed in the prospective compared with the retrospective cohort (p < 0.0001). Figure 1 shows the mean Ct values for two probes of different lengths for RPL11, RPL24 and GNB2L1. Ct values were significantly higher (significance levels indicated in Fig. 1) with the longer probes, except for GNB2L1 in the prospective samples. Shorter length gene probes for these genes were selected for subsequent analyses.

Effects of FFPE block age and gene probe length.
Gene stability. Table 3 lists the candidate endogenous genes by stability as determined by CoV, GeNorm and NormFinder in the prospective (n = 12) and TCGA (n = 401) bladder cohorts. Figures 2 and 3a plot the candidate endogenous gene stability ranking by GeNorm in the prospective and TCGA cohorts respectively. GeNorm also defines the stability value for the combination of the two most stable genes. All candidate genes were expressed stably with values below the recommended M = 1.5 cut-off. In the prospective cohort GeNorm identified the combination of SDHA and IPO8 as the most stable. Five genes (UBC, RPLP0, HMBS, GUSB, and TBP) were present in all the ten most stable genes ranked by CoV, GeNorm and NormFinder. In the TCGA cohort PPIA and TBP had the greatest stability by both CoV and NormFinder. However, GeNorm ranked the four co-expressed genes RPLP0, RPL11, RPL24 and GNB2L1 (Fig. 3) as exhibiting the most stable expression. Interestingly the next two genes were PPIA and TBP. www.nature.com/scientificreports/  www.nature.com/scientificreports/ To explore the possibility of bias due to co-expression GeNorm analysis was carried out omitting 3 of the 4 co-expressed genes. Figure 3 shows that when 3 of the 4 co-expressed genes were excluded from the analysis, expression levels of PPIA and HPRT1 were the most stable with TBP and HMBS in third and fourth place. Co-expression accounts for some of the apparent high stability of these four genes when analysed collectively alongside all candidate genes by GeNorm. However, when analysed in the absence of co-expressing genes, RPL11 and RPL24 rank sixth, RPLP0 seventh and GNB2L1 tenth suggesting that their expression is sufficiently stable to use as endogenous controls.
Assessing the performance of the selected endogenous control genes. Genes in the bladder cancer cohort extracted from the TCGA were ordered by CoV and the candidate endogenous genes highlighted (Fig. 4). All the candidate endogenous control genes fell within the lower 50% of CoV values. The most stable Table 3. Gene stability ranking: The ten most stable genes selected on the basis of lowest CoVs, or by inputting Ct values into GeNorm and Normfinder algorithms from the prospective bladder cancer cohort and TCGA bladder cancer cohort. (*' common' refers to genes which appear in the top ten most stable genes in all three measures of stability (CoV, GeNorm and Normfinder)).    TFRC  B2M  PGK1  YWHAZ  ACTB  UBC  IPO8  SDHA  GAPDH  POLR2A  GUSB  HMBS  HPTR  TBP  PPIA  RPLP0  GNB2L1  RPL11

Discussion
Measuring gene expression is increasingly important for a diverse range of clinical applications [23][24][25][26] . Purification of RNA, an essential prerequisite for qPCR, removes other cellular components and data must be normalised based on the stable expression of endogenous control genes. Analysis of gene expression studies showed using a single endogenous control gene 11 can produce gene expression error values of over 20-fold suggesting that multiple endogenous control genes are required for normalisation. Selection of endogenous control genes needs to be rigorous and take account of potential confounding factors which may be study specific. Commonly used endogenous control genes shown to be stable in one tissue type and set of conditions may be unsuited for others. GAPDH and β-actin and suitable for qPCR normalisation in some situations because they are expressed at high and constant levels in many cells and tissues 27,28 . However, a study using GAPDH and ACTB as endogenous control genes demonstrated aberrations in qPCR results due to the regulatory effects of microRNAs on the expression of these genes 7 . Further, the expression of GAPDH correlates highly with CA9, a marker of hypoxia, precluding the use of GAPDH as an endogenous control gene for studies involving solid tumours.
In this study we have described a workflow that uses a combination of laboratory and software tools to select a set of endogenous control genes for qPCR studies. The protocol is summarised in Fig. 6.
Both GeNorm and Normfinder can identify the most stable from the least stable endogenous control genes 29 . However, in common with other studies 30,31 the order of gene stability ranking by the two algorithms differed for both the prospective and TGCA cohorts. GeNorm uses pairwise comparison of candidate endogenous control genes to test for gene expression stability 11 to stepwise eliminate the least stable genes. NormFinder is a mathematical model based on ANOVA which calculates an overall average expression level for all genes to which it compares the expression of each individual gene and ranking according to stability 12 . For small sample size studies GeNorm is recognised as the more reliable algorithm for determining gene expression stability 29 . On the other hand, Normfinder is considered more robust than GeNorm for studies with larger sample numbers. GeNorm can preferentially select genes that are coregulated which mutually reinforce and so bias the apparent expression stability of co-expressed genes.
Using GeNorm, POLR2A and HMBS were identified as the most stable gene combination in the prospective sample cohort. GeNorm but not Normfinder ranked the four co-expressed genes as the most stable of the candidate group when analysed in the TCGA cohort. To test the possibility that GeNorm was selecting these genes due to bias through co-expression the analysis was repeated excluding 3 of the 4 co-expressed genes in our endogenous control gene panel. In each case the remaining gene remained within the top ten genes but ranked lower. This finding suggests that these genes can still be used for normalisation but when analysed together using GeNorm multiple co-expressed genes can provide mutual reinforcement in stability score which overstates their actual stability, at least in part, explaining the difference of the overall endogenous control ranking by GeNorm and Normfinder. Final selecƟon based on genes that consistently perform in the top 10 most stable genes by each analysis Confirm stability by analysing the candidate genes in the relevant tumour cohort in TCGA to assess stability and produce a visual display of the co-ordinaƟon between gene expression in each sample www.nature.com/scientificreports/ Interestingly GeNorm did not rank the co-expressed genes highly in the prospective muscle invasive bladder cancer cohort. Overall ranking of the endogenous control genes by both GeNorm and Normfinder differ between the two cohorts. These different rankings are not surprising as gene expression data in TGCA is acquired using RNAseq and normalised. Gene expression data in the TCGA is also acquired using RNA extracted from freshfrozen tissue which will be less modified than that from FFPE. However, it has been shown that RNA expression acquired using FFPE maintains the fidelity of patterns in biological signals and relationships with patient outcomes consistent with studies using fresh-frozen tissue 32 .
Taqman PCR gene expression methodology requires complete hybridisation of gene probes to cDNA sequences to register a hit which would suggest that shorter gene probe sequences will improve gene expression detection especially in degraded samples. Consistently Ct values were found to be significantly lower compared with long probes when using shorter gene probes for RPL11, RPL254 and GNB2L1 demonstrating the importance of choosing shorter length probes to reduce the risk of sample gene dropout especially when using archived FFPE samples.
In summary, our work highlights the importance of probe length and the need to account for co-expression when selecting a panel of endogenous control genes for downstream application in clinical samples. We identified a set of six genes that are stably expressed in FFPE bladder cancer samples and are suitable for use as endogenous control genes. We also recommend use of our workflow to harmonise the process of endogenous control selection qPCR-based studies.