Assessment of suitable reference genes for RT–qPCR studies in chronic rhinosinusitis

Reverse transcription–quantitative polymerase chain reaction is a valuable and reliable method for gene quantification. Target gene expression is usually quantified by normalization using reference genes (RGs), and accurate normalization is critical for producing reliable data. However, stable RGs in nasal polyps and sinonasal tissues from patients with chronic rhinosinusitis (CRS) have not been well investigated. Here, we used a two-stage study design to identify stable RGs. We assessed the stability of 15 commonly used candidate RGs using five programs—geNorm, NormFinder, BestKeeper, ΔCT, and RefFinder. Ribosomal protein lateral stalk subunit P1 (RPLP1) and ribosomal protein lateral stalk subunit P0 (RPLP0) were the two most stable RGs in the first stage of the study, and these results were validated in the second stage. The commonly used RGs β-actin (ACTB) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) were unstable according to all of the algorithms used. The findings were further validated via relative quantification of IL-5, CCL11, IFN-γ, and IL-17A using the stable and unstable RGs. The relative expression levels varied greatly according to normalization with the selected RGs. Appropriate selection of stable RGs will allow more accurate determination of target gene expression levels in patients with CRS.


Results
Reference gene expression profiles. The PCR efficiency of each RG is shown in Supplementary Table S1.
The slopes of the standard curves ranged from −3.280 to −3.107, the efficiencies from 101.8 to 109.8%, the correlation coefficients (R 2 ) from 0.975 to 0.999, and the intercepts from 23.812 to 33.082. The melting curves revealed single peaks and no signal was detected in the negative controls for all primer pairs.

Analysis of gene expression stability.
To determine the stability of RGs, the 15 RGs were analysed using geNorm 7 , Normfinder 15 , BestKeeper 16 , and ΔCT 17 . RefFinder 18 , a web based comprehensive evaluation platform, was then used to calculate an overall final ranking based on the results from the above four different algorithms. To evaluate anatomical variations in expression, we divided the 39 patients in the first population into nasal polyp, uncinate process, and control groups. In the second stage of the study, we analysed an independent group of 36 CRS patients to validate the results of the first stage. Additionally, we further classified CRS with nasal polyps (CRSwNP) into eosinophilic CRSwNP (ECRSwNP) and non-eosinophilic CRSwNP (NECRSwNP) according to previously published criteria 19 to evaluate the influence of eosinophilic inflammation. The patient characteristics are shown in Supplementary Table S2. geNorm analysis. According to the analysis of all of the samples in the first stage, RPS18 had the lowest expression stability (M value) that means the most stable gene expression, followed by RPLP0, RPLP2, ATP5F1, and RPLP1 (Fig. 1A). TFRC, ACTB, B2M, GAPDH, and GUSB were the five least stable genes. Of these RGs, RPS18, RPLP0, and RPLP1 were validated as being among the five most stable RGs in the second stage of the study, and TFRC, ACTB, and GAPDH were validated as being among the five least stable RGs ( Fig. 2A). The optimal number of reference genes was also determined using geNorm. The V 2/3 value was below the threshold of 0.15 in all samples, indicating that two genes were sufficient for normalization (Figs 1B and 2B). In the subgroup analyses, RPLP0, RPLP2, and RPS18 showed the highest stability in all tissues in both stages of the study (Figs 1C and 2C). The V 2/3 value was below the threshold of 0.15 in each subgroup.
NormFinder analysis. We also calculated stability values using NormFinder software. The ranking of the five most stable candidate genes is shown in Table 2. Based on analysis of all the samples, TBP was the most stable RG, followed by PGK1, PPIA, RPLP1, and HPRT1, while TFRC, B2M, ACTB, GAPDH, and GUSB were the least stable genes. In the second study stage, PPIA and RPLP1 were validated as being among the five most stable RGs and TFRC, B2M, ACTB, and GAPDH as the least stable RGs. BestKeeper analysis. BestKeeper identified the five most stable RGs listed in Table 3. RPLP2, ATP5F1, RPS18, and RPLP0 were validated as being among the five most stable genes in the second step. TFRC, ACTB, GAPDH, and TBP were validated at the least stable RGs.
ΔCT analysis. The RG stability according to ΔCT analysis is shown in Table 4. Among all of the samples, RPLP1, PPIA, ATP5F1, RPLP0, and HPRT1 were the most stable RGs, while TFRC, B2M, ACTB, GAPDH, and GUSB were the least stable RGs. In the second stage of the study, RPLP1, PPIA, and RPLP0 were validated as having high stability and TFRC, B2M, ACTB, GAPDH as having lower stability.
RefFinder analysis. An overall final ranking was calculated based on the rankings from the previous four different programs, and is shown in Table 5. The comprehensive ranking from the most to the least stable expression is as follows: RPLP1, RPLP0, RPLP2, PPIA, ATP5F1, RPS18, TBP, PGK1, HPRT1, GUSB, YWHAZ, GAPDH, ACTB, B2M, and TFRC. In the second stage of the study, RPLP1, RPLP0, and RPLP2 were validated as being among the five most stable RGs, whereas TFRC, ACTB, and GAPDH were validated as being among the least stable RGs. In the subgroup analyses, RPLP0 was among the two most stable RGs in most of the subgroups, but showed moderate stability in the uncinate process subgroup in the first study stage. RPLP2 showed high stability in all of the subgroups. Conversely, TFRC, ACTB, and GAPDH were the least stable RGs in all of the subgroups. We were not able to detect any apparent trends associated with different anatomical sites or inflammatory differences.
Influence of reference gene choice on the relative expression of target mRNA. To evaluate the influence of the RGs, the expression patterns of IL-5, CCL11, IFN-γ, and IL-17A were evaluated because CRS showed mixed enhanced Th1/Th2/Th17 reactions 14 . We chose the most stable RGs, RPLP0 and RPLP1, and conventional, commonly used, but less stable RGs, ACTB and GAPDH, for normalization. The gene expressions were determined in the population that were analysed in the first part of the study. RG selection affected IL-5 gene expression patterns between the subgroups as follows: RPLP0 (p = 0.0181) and RPLP1 (p = 0.0503), but ACTB (p = 0.0598) and GAPDH (p = 0.1675), assessed by Kruskal-Wallis test (Fig. 3A). CCL11 gene expression patterns were also affected by RGs as follows: RPLP0 (p = 0.0042), RPLP1 (p = 0.0026), and ACTB (p = 0.0029), but GAPDH (p = 0.1390) (Fig. 3B). However, the expression pattern using ACTB was different from using RPLP0 and RPLP1. Although IFN-γ and IL-17A did not show the differences, the genes exhibited similar expression trends when we normalized using RPLP0 and RPLP1. Conversely, when ACTB and GAPDH were used, particularly wide

Discussion
Here, we investigated the stability of candidate reference genes for RT-qPCR in nasal polyps and sinonasal tissues using a two-stage study design. Biological and experimental errors introduced throughout the RT-qPCR process need to be accounted for 5,6 , and normalization using an RG is a simple and popular method of internally controlling for such errors, as well as controlling for different input RNA amounts during the reverse transcription step 5 . High RNA integrity and purity are critical for obtaining meaningful and reliable gene expression data and ensuring reproducibility of results 2,20 . Poor RNA integrity may generate misleading differences in gene expression   measurements 1,20 . Stable RGs in nasal polyps and sinonasal tissues differ between intact and degraded RNA samples 21 . In this study, we chose samples with RNA integrity number (RIN) values ≥7 to avoid false results. We used four different programs, geNorm, NormFinder, BestKeeper, and ΔCT to identify stable RGs. Unfortunately, we did not identify any RGs that were universally stable across these four programs. The discrepancies in the results occurred because the different programs use different algorithms, for example a pairwise comparison or a model-based approach 7,15-17 . To overcome the discrepancies and obtain a final ranking, we used RefFinder software. Two RGs, RPLP1 and RPLP0, were identified as the two most stable RGs in the first stage of the study, and this was validated in the second stage. Our findings indicate that RPLP0 and RPLP1, either singly or in combination, are suitable for normalizing gene expression in nasal polyp and sinonasal tissues. RPLP0 has been used as an RG in previous CRS studies [22][23][24] . However, this study is the first to demonstrate its stability in nasal polyp and sinonasal tissues. Conversely, TFRC, ACTB, and GAPDH were among the five least stable genes throughout all of the algorithms, including RefFinder. Based on these results, these conventionally used RGs should not be used.
Narrower standard deviations and the different p-values in the gene expressions were revealed when we used the most stable RGs, RPLP0 and RPLP1, for normalization, compared with when we used the commonly used, less stable genes GAPDH and ACTB. The selection of RGs could shift results from indicating significant differences to being non-significant, and vice-versa. Normalizing using GAPDH, especially, produced a different expression trend. GAPDH mRNA expression levels are known to differ between different tissues and between the same tissues in different individuals 25 . In this study, GAPDH showed lower stability in all of the anatomical sites and inflammation patterns studied compared to other RGs. It has been reported that GAPDH and ACTB were not suitable as RGs for quantitative analysis of gene expression in asthma, which has similar pathophysiology to CRS 26 . These results emphasise the importance of choosing stable RGs for normalization.
A large number of studies have investigated the validation of reference genes in many different tissues and cell types. However, to the best of our knowledge, few have examined the suitability of reference genes in CRS. Perez-Novo et al. 21 investigated 16 samples, including degraded samples, from ethmoid and maxillary sinuses from patients with nasal polyps and CRS. In intact RNA, they found that the genes for hydroxymethyl-bilane synthase (HMBS) and succinate dehydrogenase complex, subunit A (SDHA) in CRS, and ACTB and TBP in nasal polyps were the most stable among nine candidate reference genes analysed using geNorm. GUSB and ATPase

Total
Nasal   plasma membrane Ca 2+ transporting 4 (ATP2B4) were identified as reliable genes for normalization of cystic fibrosis transmembrane conductance regulator (CFTR) gene expression in the nasal mucosa and nasal polyps in patients with cystic fibrosis 27 . These RGs differ from those identified in the current study. These differences may be attributed to differences in sample sizes, sampling location, differences in pathophysiology within CRS subgroups, or ethnic differences between Western and Japanese populations [28][29][30] .

Conclusion
This study identified suitable RGs for normalizing target gene expression levels in nasal and sinus tissues using a two-stage study design. RPLP0 and RPLP1, either singly or in combination, are suitable for normalizing gene expression in nasal and sinus tissues, whereas TFRC, ACTB, and GAPDH were less stable RGs according to all of the algorithms used. Use of appropriate reference genes will facilitate the generation of accurate, robust, and reproducible gene expression studies in CRS.  and clinical characteristics were obtained from the patients prior to surgery, including age, sex, asthma status, and history of sinus surgery. Preoperative computed tomography scans were assessed according to the classification described by Lund and Mackay 31 . The preoperative polyp-grading system used a 5-point scale (score of 0-4) according to the recommended guidelines 32 . Blood samples were taken before surgery, and complete blood counts and serum total IgE levels were determined.

Patients
Sampling and total RNA extraction. At the surgery, we removed nasal polyps from patients with CRS, and uncinate processes from CRS patients with/without nasal polyps and control patients. The control group consisted of 10 patients with pituitary tumours or anatomical variants, without endoscopic and radiological evidence of sinus disease, in the first population. Under 0-degree endoscope, the samples were obtained with a scalpel and Grunwald's forceps without local anaesthesia instillation. The tissues were immediately immersed in RNA later (Ambion, Austin, TX, USA) and stored at 4 °C for 1-2 days, then at −80 °C until analysis. We extracted RNA using NucleoSpin RNA (Macherey-Nagel, Düren, Germany) in the first stage of the study and an miRNeasy Mini Kit (Qiagen, Hilden, Germany) for the second stage, both according to the manufacturer's instructions. The quality and quantity of the extracted RNA were determined by measuring the ratios of absorbance at 260/230 nm and 260/280 nm using a DeNovix DS-11 spectrophotometer (DeNovix, Wilmington, DE, USA) and a NanoDrop 2000 spectrophotometer (Thermo, Wilmington, DE, USA). RNA integrity was confirmed with an RNA 6000 Nano Chip using an Agilent 2100 Electrophoresis Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). An RIN ≥7 was considered adequate for analysis, on a scale where 1 indicated the most degraded and 10 the most intact profile. qPCR. A sample of total RNA was reverse transcribed into cDNA using PrimeScript RT Master Mix (Takara, Shiga, Japan) in the first stage of the study and an iScript cDNA Synthesis kit (Bio-Rad, Hercules, USA) for the second stage, according to the manufacturers' instructions. Reverse transcription was performed in a TaKaRa PCR Thermal Cycler Dice Gradient (Takara) and iCycler Thermal Cycle system (Bio-Rad). cDNA was stored at −20 °C until qPCR experiments were performed.
We examined the expression levels of the 15 reference genes provided in the Human Housekeeping Gene Primer Set (Takara) ( Table 1). All of the PCR products ranged from 75 to 200 bases. qPCR amplification reactions were performed using a Light Cycler 96 (Roche, Mannheim, Germany) for the first stage of the study and a CFX96 Touch Real-Time PCR Detection System (Bio-Rad) for the second. Amplifications were performed with 30 s enzyme activation at 95 °C, followed by 40 cycles of denaturation at 95 °C for 5 s, and then annealing/extension at 60 °C for 20 s. At the end of each run, melting curve analysis was performed from 65 °C to 95 °C. Briefly, 2 μl of cDNA equivalent to 10 ng of total RNA was used as a template in a total reaction volume of 10 μl containing 5 μl of SYBR Premix Ex Taq II (Takara), 200 nM of each primer, and RNase/DNase-free water.
The slope, efficiency, correlation coefficient (R 2 ), and intercept of each primer pair were determined from the standard curve created using 5-point serial dilutions of cDNA template mixture, and analysed by qbase plus (Biogazelle, Ghent, Belgium). The absorbance ratios (mean ± standard deviation) at 260/230 nm (2.07 ± 0.14) and 260/280 nm (2.13 ± 0.08) indicated that the RNA samples were pure and free of protein. The mean RIN values were 8.04 ± 0.58.
Statistical analysis. The stability of gene expression was analysed using geNorm 7 , Normfinder 15 , BestKeeper 16 , and ΔCT 17 . RefFinder 18 was then used to calculate an overall final ranking based on the results from the above four different algorithms. Raw Cq values were analysed in qbase plus and the expression stability (M value) of each reference gene was calculated as the average pairwise variation for the reference gene with all of the other genes. A low M value represented stable gene expression 7 . Additionally, the pairwise variation (V n/n+1 ) between the two sequential normalization factors (NF n and NF n+1 ) was calculated to define the optimal number of reference genes. When the V value was below the cut-off of 0.15, it was not necessary to include additional reference genes 7 . The NormFinder add-in for Microsoft Excel was used to estimate both intragroup and intergroup variation for each candidate reference gene. These two variation values were combined to give a stability value, with the lowest value indicating the most stable expression 15 . BestKeeper calculated RG stability based on the standard deviation (SD) of Cq values; RGs with SD > 1 were excluded 16 . ΔCT generated 'pair of genes' comparisons between each RG and the other RGs within each sample and calculated the average SD against the other RGs 17 . Finally, based on the rankings of these four algorithms, RefFinder, a web based comprehensive evaluation platform, was used (http://www.leonxie.com/referencegene.php).
Statistical analysis was performed using IBM SPSS Statistics version 23 (IBM Corporation, Armonk, NY, USA) and graphs were produced using GraphPad Prism 6 for Mac OS X (GraphPad Software Inc., San Diego, CA, USA). Continuous and categorical data were compared among subgroups using Kruskal-Wallis and χ 2 tests. The Kruskal-Wallis with post-hoc Dunn's multiple comparison test was used for comparing more than two groups. Nasal polyp score was compared between two groups using Mann-Whitney U-tests. A value of p < 0.05 was considered to indicate a statistically significant difference. Data Availability Statement. All data generated or analysed during this study are included in this published article (and its supplementary information files).