Amelioration for an ignored pitfall in reference gene selection by considering the mean expression and standard deviation of target genes

Routine tissue-specific reference genes are often used in expression studies, but target genes are not taken into account. Using the relative RT-qPCR approach, we evaluated the expression of three target genes. At the same time, meta-analyses were conducted in various ethnic groups, genders, and thyroid cancer subtypes. When eight common reference genes were examined, it was discovered that some of them not only lacked consistent expression but also had considerable expression variance. It is worth noting that while choosing a reference gene, the mean gene expression and its standard deviation should be carefully addressed. An equation was developed based on this, and it was used to perform statistical analysis on over 25,000 genes. According to the subtype of thyroid cancer and, of course, the target genes in this investigation, appropriate reference genes were proposed. The intuitive choice of GAPDH as a common reference gene caused a major shift in the quantitative expression data of target genes, inverting the relative expression values. As a result, choosing the appropriate reference gene(s) for quantification of transcription data, and especially for relative studies of the expression of target gene(s), is critical and should be carefully considered during the study design.

www.nature.com/scientificreports/ A showed contradiction for 13 out of 17 PTC samples, whereas only 4 samples (PTC samples 2, 7, 16, and 17) did not show contradiction. In PTC sample 1, Gene A was normalized against GAPDH (red bar) and a negative delta-delta Cq ratio was observed. A positive delta-delta Cq ratio was also observed for gene A just when the gene was normalized against SYMPK (blue bar). The same holds true for gene A in PTC samples 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, and 15. Therefore, only when GAPDH was replaced with SYMPK did gene A show 76.5% inconsistency. Dissimilitude was also observed when gene B (Fig. 2B, 52.9% differences) and gene C (Fig. 2C, 29.4% differences) were normalized against GAPDH and SYMPK. Therefore, when a specific PTC tissue was compared to its adjacent normal tissue, target gene could be reported as overexpressed or downregulated at the same time.
SYMPK and GAPDH expression in normal and PTC tissues. The best way to normalize RT-qPCR data is to pick a reference gene or genes that exhibit the least amount of variation in mRNA expression across all of the samples. SYMPK (Fig. 3, blue bars) showed a narrower range of Cq values than GAPDH (Fig. 3, red bars) in normal tissues, and the same was true in PTC tissues, where SYMPK had less variance.
Statistics on reference genes and target genes. Statistical analyses of target genes (A, B, and C) and reference genes (GAPDH and SYMPK) are presented in Table 1. SYMPK exhibited a lower SD = 1.74 and CqCV% = 5.84 in both the adjacent normal and the PTC tissues than GAPDH (SD = 4.26 and CqCV% = 16.32). To determine the differences in the expression values of GAPDH and SYMPK in the adjacent normal tissues as well as the PTC tissues, separate tissue statistics have been provided. In adjacent normal tissues, the mean Cq values of gene A (29.00) and gene B (28.74) were close to the mean Cq value of SYMPK (29.96) but far from the corresponding Positive and negative deltadelta Cq ratios respectively represent a target gene in PTC tissues that is down-regulated or over-expressed. Y-axis present deltadelta Cq ratios and X-axis show PTC samples.  Figure 3. Expression of SYMPK and GAPDH in normal and PTC tissues. In normal tissues, SYMPK (blue bars) exhibited a smaller range of Cq values than GAPDH (red bars), and this was also true in PTC tissues, where SYMPK had less variance. www.nature.com/scientificreports/ roid carcinoma).Microarray probes were matched to corresponding genes, mean expression values for a probe set were calculated for each gene, and the data was subjected to "removeBatchEffect" (Supplementary Figs. S1 and S2). The expression levels of eight common reference genes were compared in two ways: between normal tissues and each subtype of thyroid cancer, as well as between subtype (Table 2). GAPDH and SYMPK had effect sizes (ES) of 0.235 and 0.151 for the PTC subtype, respectively, when compared to normal tissues; however, the ES of GAPDH was statistically significant (p = 0.0020). GAPDH had statistically significant ES values in both the FTC (p = 0.0012) and ATC (p = 3.19E−17) subtypes. Furthermore, GAPDH had higher ES values than SYMPK in ATC (0.652 vs 0.070 respectively) and FTC (0.389 vs 0.154, respectively) subtypes. Other subtypes, such as FTA, PDTC, and MTC, showed negligible differences between GAPDH and SYMPK expression. GUSB (− 0.024), ACTB (0.032), and HPRT1 (0.037) were the three most ideal reference genes in the PTC subtype, with the lowest insignificant ES (Fig. 4A,B). The best three reference genes for other subtypes were SYMPK (0.070), TBP (− 0.076) and GUSB (− 0.098) in ATC (Fig. 4C,D); ACTB (0.036), HPRT1 (0.063), and GUSB (0.064) in FTC (Fig. 4E,F); GUSB (0.023), HPRT1 (− 0.052), and TBP (− 0.062) in FTA (Fig. 4G,H); GUSB (0.062), HPRT1 (− 0.075), and PYCR1 (− 0.102) in PDTC (Fig. 4I,J); ACTB (0.015), B2M (0.027), and TBP (− 0.070) in MTC (Fig. 4K,L).
The inter-subtype analysis was divided into two parts: the first assessed the differential expression of reference genes between undifferentiated (ATC) subtype and all other subtypes, and the second part was devoted to assessing the differential expression of reference genes between the poorly differentiated (PDTC) subtype and differentiated subtypes (FTA, PTC, FTC, MTC). GAPDH had statistically significant differential expression between ATC and all other subtypes, with the exception of FTC (0.262) and MTC (0.354). When undifferentiated-ATC tissues were compared to differentiated-PTC tissues, the genes GAPDH, ACTB, B2M, HPRT1, and PYCR1 were found to be significantly expressed. The same results were obtained when comparing undifferentiated-ATC tissues to poorly differentiated-PDTC tissues. A gene expression analysis was also performed to compare PDTC to other differentiated subtypes, and none of the reference genes were statistically significant. As a result, only a comparison of PDTC with FTA was reported in Table 2 and the others were omitted.
Intra-sex analyses, as well as sex-subtype interactions. Intra-sex analysis was performed to determine the differentially expressed reference genes in each of the two sexes, and the interaction of sex and subtype was investigated using factorial designs (Table 3). We dealt with 253 samples, including 44 normal, 15 FTA, 119 PTC, 24 FTC, 27 PDTC, and 24 ATC, after 6 out of 14 datasets failed to offer detailed information regarding the sex of the patients. We did not have any FTA-male samples, and no MTC subtype samples were left. Most of the reference genes did not reveal statistically significant differences in expression in intra-sex analysis. The only exceptions were ATC-women, who had statistically different expression of B2M (ES = 0.536, p = 0.0175) and PYCR1 (ES = 0.900, p = 0.0290) genes. The ES value of GAPDH was higher in females than males in PTC subtype (ES = 0.222 vs ES = 0.028 respectively), but the difference was not statistically significant (ES.Female-ES.Male = 0.194, p = 1), according to the interaction analysis. There were also differences in the expression of some other reference genes between females and males (e.g. TBP in ATC and B2M or GUSB in FTC), but using a factorial design to calculate the differences in differential expression revealed no significant differences in the expression of these two genes (p = 1 and 1 or 0.4560, respectively).
The ES of reference genes were depicted in females and males based on their subtypes ( Intra-subtype, inter-sex analysis. With the exception of FTA, inter-sex analysis was performed within subtypes to determine the most appropriate reference gene in different pathological conditions (normal and subtypes, Table 4). TBP, PYCR1, and B2M were the best reference genes in normal tissues (Fig. 6A,B), while ACTB, TBP, and HPRT1, were the best ones in PTC subtypes (Fig. 6C,D). HPRT1, SYMPK, and TBP were the best genes for the FTC subtype (Fig. 6E,F), HPRT1, ACTB and GAPDH for the PDTC subtype (Fig. 6G,H), and B2M, GUSB, and ACTB for the ATC subtype (Fig. 6I,J).
Microarray and RNA-seq data statistics. The TCGA database was used to download raw expression counts of 560 samples, including 502 PTC and 58 normal tissues, and the statistics of this RNAseq data are shown in Table 5. ACTB (2.89), GAPDH (3.08), and SYMPK (3.25) were the top three genes in PTC tissues with the lowest CV% values. In normal tissues adjacent to PTC tissues, SYMPK (CV% = 2.84) was ranked after GAPDH (2.46) and GUSB (2.59). According to the differential expression of the reference genes (Table 6), the top three genes with the lowest ES values were ACTB (− 0.001), TBP (− 0.017), and SYMPK (0.034), respectively. GAPDH had the highest ES value = 0.06 among eight reference genes (Fig. 7). Table 7 shows statistics for microarray pooled data from adjacent normal tissues and each thyroid cancer subtype. While GAPDH was ranked fifth (3.37), the genes with the lowest CV% values in normal tissues were GUSB (2.77), B2M (2.86), and SYMPK (3.10), respectively. GUSB (2.48), GAPDH (2.56), and ACTB (2.86) had the lowest CV% values in PTC tissues, followed by SYMPK (3.38).
To facilitate use, the basic statistic for all 6331 genes in the GEO dataset (Supplementary Table S3) and all 25,705 genes in the TCGA dataset (Supplementary Table S4) were provided. These two tables compare the mean and standard deviation values of prospective target genes with the statistics of candidate reference genes. www.nature.com/scientificreports/

Discussion
In research and clinical detection, RT-qPCR is the gold-standard method for expression evaluation [16][17][18] . The advantageous of RT-qPCR include high sensitivity and specificity, speed of analysis, and real-time monitoring of results 8 . Nature protocols require that appropriate internal reference gene(s), formerly known as housekeeping genes, be validated prior to each study 19,20 . Historically, an ideal reference gene has minimally altered expression under various pathological and physiological conditions such as tumour type and patient sex. It must be free of pseudogene(s) and alternative splicing 15 . We previously investigated eight reference genes and discovered that SYMPK was more stably expressed than conventional reference genes (GAPDH and ACTB) and also lacked pseudogenes 15 . Ribosomal RNA (18S rRNA) is a highly recommended reference gene for RT-qPCR data normalization 21,22 . Unfortunately, 18S rRNA has at least three drawbacks: inhibition by mitomycin C 23 , absence in bulk high-throughput expression platforms, and a clear role in cancer development [24][25][26][27][28] and prognosis 29 . We did not include 18S rRNA in our study due to the aforementioned facts and a previous report about its unstable expression 30 .
GAPDH and SYMPK were used as reference genes to normalize three candidate genes to better understand the consequences of using inappropriate reference genes. GAPDH was chosen because it is the most commonly used reference gene in molecular biology, and we previously reported it as the worst reference gene using Nor-mFinder algorithm 15 . This is in line with a previous study that found GAPDH to be unsuitable for normalizing relative RT-qPCR data from bladder and colon cancer 31 . The gene did not meet the criteria of those authors (e.g. tissues stability, expression level above background, and lack of alternative splicing), so it was eventually ignored despite being ranked in colon cancer. PLA2R1   IPCEF1   DPP6  TFF3  DLG2  DPP4  MPPED2  GGCT  OCA2  RXRG  GPM6A  TUSC3  TBC1D4  BID  GHR  ETV5 PLAU  MYOC  KIT  CTSH  CLUL1  MET  DIO1  SLITRK5  CORO2A  ABCC3  WASF3  CITED1  LGALS1  FMOD  FHL1  CHI3L1  GSTM3  DUSP5  TNFRSF11B  CD55  TGFBI  ZNF536  AGR2  EDN3  PNP  NPY1R  IGFBP6  BBOX1  FCER1G  CCN1  DIRAS3  LTF  SLPI PEBP1   NEBL   PRKCQ  LRP2  KIF14  BCL2  TMEM158  CDH1  MME  HHEX  STEAP1  FCGBP CDS1  MKI67 TFPI2  ARHGAP6  LTBP1  MPPED2 TXNL1  GTSE1  WASF3  CDKN3  CLIP2  SLC6A13  PLAU  AURKA  DLG2  KIT  VCAN  CD14  CXADR  PCP4  ADAMTS2  ADD3  CDH16  SMC4 CXCL8  IRS4  CD300A RRM2  MMRN1 CHKA  KYNU SPP1  SCNN1B  NUDT1  GDF10  SEMA3C  STAC  WFDC2  FOSB  MYO1F  AVPR1A  BGN  ADIRF  MMP9  IL1B  COL9A3 EPHB6  S100A2  ST3GAL5 HLA    www.nature.com/scientificreports/ In this study, the expression of reference genes (GAPDH and SYMPK) was compared between normal and PTC tissues, SYMPK was found to be a better reference than GAPDH because it had less variability. Aside from the lack of alternative splicing, lower CqCV% values for SYMPK gene were obtained from relative RT-qPCR Table 3. Intra-sex analyses, as well as sex-subtype interaction. Using three scenarios, the expression levels of eight common reference genes are compared between normal tissues and each subtype of thyroid cancer. In the first scenario, a linear model was used to calculate differential expression in females: expression ~ tumor, where tumor was a binary variable (tumor vs. normal). In the second scenario, male differential expression was calculated using the same linear model as in the first. In the third scenario, the differential expression of differences was calculated using a complex linear model: expression ~ group, where group was a single factor made up of sex and subtypes. As a result, the binary variable was female vs. male, and the contrast was (female. tumor-female.normal)-(male.tumor-male.normal). For intra-sex and interaction analyses, ES and FWER are presented separately. The interaction could not be calculated because there was no male with FTA subtype. ES effect size, FWER family-wise error rate.  RASSF9  RYR2  GGCT  IPCEF1  MMRN1 MYOC  MPPED2PLA2R1   NPC2  GHR  FN1  UPP1  GPC3  SLC4A4  ETV5 PROS1  EPHB1  LRP4  ENTPD1  TFF3  S100A11  AVPR1A  KIT  WASF3  CORO2A  DIO1  CLUL1  CHI3L1  TGFA  CHRDL1 TTC30A  CLDN10  MDK  APOD PROX1  ADM  SLC26A4  CITED1  CTSC  BMP8A  B3GNT3  ID1   0   5 CASP1  IL7  ALDH5A1  ADAM8 GNA15  NEBL  SEMA3C  MME  MFAP2  BMP7  CKB  NFKBIE  TLE2  PSMB9  MYO1F  NR3C2  BCL2A1  TSHR  RHOG  SASH1  COL7A1 GREM1  PRF1  LRP2  PEBP1  SALL2  FCGR3B  MMP1  CHEK1  EVI2B  C2  TJP3  MAOA  MCM5  RND3  IRS4  GTSE1  RGN  CXCL13  CD3D  TRIM2  CASP3  FAM155B DIO2  ITGB2  www.nature.com/scientificreports/ data in both normal and PTC tissues. The main point of contention is that GAPDH had a significantly higher SD than the target genes, a flaw that makes it decidedly inappropriate for mRNA expression normalization. We performed a meta-analysis on GEO microarray data combined with a comprehensive TCGA RNA-seq data analysis to increase the sample size, include all thyroid cancer subtypes, and involve both sexes. We discovered that GAPDH was significantly upregulated in PTC, FTC, and ATC, and as a result, the gene is unsuitable as a reference gene according to the microarray meta-analysis. GAPDH was found to be significantly upregulated at various stages of tumor differentiation. This idea suggests that GAPDH may be a key promoter of tumor aggressiveness, as previously reported by Chiche et al. in non-Hodgkin's B lymphomas 32 . They proposed that the increased GAPDH levels activated the nuclear factor-κB gene, which in turn increased the activity of hypoxia-inducing factor-1α (HIF-1α). In this study, when FTA and ATC subtypes were compared, the expression of HIF-1α was also upregulated (ES = 0.497, p = 0.0001). www.nature.com/scientificreports/ We provided separate tables to assist researchers in accurately selecting reference genes for their study designs. For example, if researchers want to study different subtypes, Table 2 provides a list of genes, and the gene with an ES closer to zero is the best fit for their research. Researchers could use Table 3 to include the gender of patients in an analysis, and the best genes are those with ES.Female -ES.Male closer to zero. Table 4 is the best reference when a specific subtype is required as well as the gender of the patients, with genes with ES values closer to zero serving as the best reference genes.
Furthermore, a discrepancy was discovered when each target gene was normalized against two different reference genes, SYMPK and GAPDH. (Fig. 2 and Table 8). We hypothesized that the differnce was occurred because of the overlap between the Cq values of target and reference genes. By overlapping, we mean that the Cq values of the reference and target genes are within the same range, and thus samples with positive ddCq mutually neutralize samples with negative ddCq, resulting in a change in the overall expression pattern of a target  abs: absolute value, µ: mean, σ (SD): standard deviation, T: Target gene expression in each subtype, R: Reference gene expression in each subtype. Consider the case where a reference gene has no variation in its expression (σ = 0) and a target gene has σ = 1. If the difference in mean expression between the target and the reference is at least three times the absolute value of the target gene's SD (3σT), the reference gene does not overlap with the target gene ( Supplementary Fig. S4a). However, a reference gene with an SD value of 0.25 necessitates a difference of at least 3.5 units between the reference and target genes' mean expression values (Supplementary Fig. S4b). By doubling the SD value of the reference gene (from 0.25 to 0.5 and from 0.5 to 1), the mean expression values of the reference and the target genes must differ by 4 ( Supplementary Fig. S4c) and 5 units (Supplementary Fig. S4d) respectively. To avoid overlap, we found that twice the absolute value of the SD of the reference gene (2σR) must also be considered for the calculation of the difference between the mean expression of the reference and the target genes. Therefore, it is possible to avoid overlaps between the expression values of reference gene and target gene and stop contradictory gene expression patterns by using Eq. (1).
For all expressed genes in GEO and TCGA, we provided tables with basic statistics, such as mean and SD (Supplementary Tables S3, S4). Our expression data could be a reliable estimate of any population for researchers to compare the mean and SD of desired genes in the above equation because our analyses include large sample sizes representing multiple ethnicities and subtypes in both sexes.
In conclusion, selecting reference gene(s) solely on the basis of specific tissues may result in inaccurate or misleading information. We questioned the common practice of selecting traditional reference genes. In a comprehensive investigation of thyroid cancer subtype, we discovered that GAPDH was significantly influenced by (1) abs (µT − µR) ≥ (3σ T ) + (2σ R) Table 6. TCGA differential expression analysis in PTC samples. PTC samples from the TCGA dataset were analysed for differential expression. ACTB (ES = − 0.002), TBP (ES = -0.017) and SYMPK (ES = 0.035) are the three most stable reference genes in PTC tissues. ES effect size, FWER family-wise error rate.  www.nature.com/scientificreports/ RNA extraction and assessment. Total RNA was extracted from RNAlater-treated samples using a onestep RNA extraction reagent (Bio Basic, Markham, ON, Canada), as directed by the manufacturer. The concentration of isolated RNA was determined using a NanoDrop OneC spectrophotometer (Thermo Scientific, Waltham, MA, USA). A260/A280 and A260/A230 ratios were used to determine RNA purity. The integrity of the RNA was determined using 1.0% agarose gel electrophoresis.
Complementary DNA (cDNA) synthesis. DNase I treatment (Thermo Scientific, Bremen, Germany) was used to remove residual genomic DNA contamination, as directed by the manufacturer. One microgram of total RNA was reverse transcribed in a total reaction volume of 20 μL using the Thermo Scientific RevertAid Reverse Transcriptase kit (Thermo Scientific, Bremen, Germany) according to the manufacturer's instructions.

Design of exon-junction primers.
To avoid amplifying genomic DNA and/or heterogeneous nuclear RNA, all primers were exon junctioned. Beacon Designer 8.1 (Premier Biosoft International, Palo Alto, CA, USA) was used to design primers that span specific exons. Oligo 7 was used to recheck the primers for any unwanted secondary structure (Molecular Biology Insights, Colorado Springs, CO, USA). The NCBI-primer BLAST service was used to confirm the specificity of the designed primers. The melting temperature of the primers was validated using temperature gradient PCR (Sinaclon Bioscience, Tehran, Iran). All of the information on the primer pairs is presented in Supplementary Table S1.
Relative RT-qPCR. In a Bio-Rad Chromo4 device (Bio-Rad, Hercules, CA, USA), a relative RT-qPCR reaction was performed using SYBR Green RealQ Plus 2 × Master Mix (Ampliqon, Odense, Denmark). The RT-qPCR reaction protocol consisted of (i) one cycle of enzyme activation and initial denaturation at 95 °C for 15 min, and (ii) 40 cycles of denaturation at 95 °C for 30 s, annealing for 30 s, and extension at 72 °C for 30 s. After each cycle, the plates were read. All relative RT-qPCR reactions were run in triplicate, with non-template control (NTC) per gene. Table 8. Single PTC sample analysis using ddCq method. Normalization of Gene A against SYMPK reveals that 11 of the 17 PTC samples have negative ddCq, while the remaining six have positive ddCq. When Gene A is normalized against GAPDH, the results are flipped, with six PTC samples exhibiting negative ddCq and eleven exhibiting positive ddCq. Gene B and Gene C also have opposing patterns. PTC samples with positive ddCq mutually neutralize PTC samples with negative ddCq, resulting in a change in the overall expression pattern of a target gene. Each PTC sample is compared with its adjacent normal tissue.ddCq: delta-delta Cq. www.nature.com/scientificreports/ Melt curve analysis. To assess the specificity of relative RT-qPCR, the melt curve was constructed by observing the gradual rise of temperature in 1 °C increments from 55 to 95 °C, followed by plate reading. The temperature (°C, x-axis) was plotted against the derivative of fluorescence change over temperature (y-axis).

Gene expression analysis. Cq values were exported from the Bio-Rad Chromo4 thermocycler into
Microsoft Excel (2013) for further analysis. The average of Cq values for reference and target genes in PTC tissues and adjacent normal tissues was calculated and the Livak method was used for normalization 2 . The delta Cq values were calculated by subtracting the Cq values of a reference gene and a target gene from each sample, and delta-delta Cq was determined by the difference between each PTC tissue and the average of delta Cq in adjacent normal tissues.
Statistical analysis. Microsoft Excel 2013 (Microsoft, Redmond, WA, USA) was used to calculate qPCR fold change, maximum Cq, minimum Cq, standard deviation (SD), mean Cq, and correlation of variation (CqCV%, CqCV% = SD/mean × 100). CV% is a statistical measure that represents the relative dispersion of gene expression values in a dataset, regardless of the mean expression values of the genes. It is used to circumvent the problematic investigation of SD without considering the overall expression.
Data collection. The GEO and The Cancer Genome Atlas (TCGA) databases were used to obtain microarray and RNAseq data, respectively. To scavenge any microarray expression data related to thyroid neoplasm, the GEO database was mined for the keywords "thyroid neoplasm", "thyroid cancer", and "thyroid carcinoma". Exclusion criteria were used, and any data from species other than Homo sapiens was discarded. Cell lines, treatments, therapies, knocked-in and knocked-out models, and any dataset with incomplete phenotype information were excluded from further analysis. To reduce other biases, samples were collected from different countries and from people of various ethnicities. To compensate for the small sample size in different sexes and pathological subtypes, pooled data analyses were performed. As a result, 14 microarray datasets containing 520 samples were used in this study. FTA, PTC, FVPTC, FTC, MTC, PDTC, and ATC were among the thyroid neoplasms represented in the datasets. Microarray datasets are described in detail in supplementary Table S2.