Identification of the best housekeeping gene for RT-qPCR analysis of human pancreatic organoids

In the last few years, there has been a considerable increase in the use of organoids, which is a new three-dimensional culture technology applied in scientific research. The main reasons for their extensive use are their plasticity and multiple applications, including in regenerative medicine and the screening of new drugs. The aim of this study was to better understand these structures by focusing on the choice of the best housekeeping gene (HKG) to perform accurate molecular analysis on such a heterogeneous system. This feature should not be underestimated because the inappropriate use of a HKG can lead to misleading data and incorrect results, especially when the subject of the study is innovative and not totally explored like organoids. We focused our attention on the newly described human pancreatic organoids (hPOs) and compared 12 well-known HKGs (ACTB, B2M, EF1α, GAPDH, GUSB, HPRT, PPIA, RNA18S, RPL13A TBP, UBC and YWHAZ). Four different statistical algorithms (NormFinder, geNorm, BestKeeper and ΔCt) were applied to estimate the expression stability of each HKG, and RefFinder was used to identify the most suitable genes for RT-qPCR data normalization. Our results showed that the intragroup and intergroup comparisons could influence the best choice of the HKG, making clear that the identification of a stable reference gene for accurate and reproducible RT-qPCR data normalization remains a critical issue. In summary, this is the first report on HKGs in human organoids, and this work provides a strong basis to pave the way for further gene analysis in hPOs.


Introduction
In the last few years, there has been a remarkable increase in the use of organoid technology in scientific research. Organoids are three-dimensional cellular bodies that can be generated from pluripotent stem cells (embryonic or induced) or adult tissue-specific ancestor cells [1]. The main reasons for their considerable use are their plasticity and multiple applications, including as new models to study specific diseases [2], in regenerative medicine [3] and the screening of new drugs [4]. Many types of organs can be recapitulated in various types of organoids, including the gut [5], brain [6], kidney [7] and liver [8]. Human organoids have been widely studied for a variety of purposes. In particular, human pancreatic organoids (hPOs) have been characterized and presented to be a potential source of functional cells for the treatment of type 1 diabetes [9,10].
For any type of hPO application such as a three-dimensional platform for tissue regeneration, disease modelling, and drug screening, the gene expression of specific tissue-related genes is a crucial requisite to define their identity. To monitor gene expression, RT-qPCR is often the method of choice due to its sensitive and specific detection, potential for high throughput, rapid and accurate quantification, and high degree of potential automation. In order to obtain gene expression results that are not only accurate but also comparable among different experimental setups, conditions, operators and laboratories, normalization of RT-qPCR data should be performed against one or two housekeeping genes (HKGs). These HKGs must display unchanged cellular expression, irrespective of the experimental conditions, and this is fundamental to achieve reliable results [11].
Moreover, the selection of the HKG for normalization of RT-qPCR data should take into account the expression stability of HKGs among different specimens and experimental conditions [11,12]. This point should not be underestimated because an inappropriate choice of the HKG can lead to misleading data and incorrect results [13,14], especially when the subject of the study is innovative and not totally explored like organoids. Unfortunately, the importance of selecting appropriate reference genes is not adequately emphasized by researchers; therefore, the validation process and the comparison of molecular biology data remain controversial and open topics.
The aim of this study was to compare 12 well-known HKGs (ACTB, B2M, EF1α, GAPDH, GUSB, HPRT, PPIA, RNA18S, RPL13A TBP, UBC, and YWHAZ; S1 Table), belonging to different functional families, on hPOs at different passages (P0, P2, P5, P7, P10) in order to identify the most stable HKGs suitable for hPO analysis. Finally, we quantified the differential expression of a panel of pancreatic markers to assess the effect of HKG variability on their transcriptional patterns. This is the first report on HKGs in human organoids, and this work provides a strong basis to pave the way for further gene analysis in hPOs.

Materials and methods
hPO isolation and culture hPOs were generated starting from adult healthy islet-depleted pancreatic tissue (gently provided by the Diabetes Research Institute, IRCCS Ospedale San Raffaele, Milan, Italy), after approval of the Institutional Review Board (National Transplant Center accredited facility IT000679), as previously described [9]. All experiments were performed according to the amended Declaration of Helsinki. The use of human specimens was approved by the Ethical Committee of Fondazione IRCCS Ca' Granda Ospedale Maggiore Policlinico n˚1982, 14th January 2020.

Flow cytometric analysis
hPOs were collected and washed with fresh AdDMEM/F12 medium supplemented with 10 mM HEPES, 1X Glutamax, 1X Penicillin-Streptomycin and centrifuged for 5 min at 300 x g to remove the Matrigel. Then, the cells were dissociated by trypsinization using Tryple (Life Technologies ref. 12604013) for 15 min to obtain a single-cell suspension. We stained the samples using 1:200 FITC mouse Anti-Human Epithelial Cell Adhesion Molecule (EpCAM; BD, cat. #347197), 1:100 Alexa Fluor1 647 Mouse Anti-SRY-box transcription factor 9 (Sox9; BD, cat. #565493) and 1:2000 FITC-conjugated lectin from Ulex europaeus (UEA1) (Sigma-Aldrich, cat. #L9006). For surface marker analysis (EpCAM and UEA1), the cells were incubated with fluorophore-conjugated antibodies for 20 min in the dark at room temperature (RT) and then washed with PBS. For intracellular marker analysis (SOX9), the cells were fixed with 4% paraformaldehyde (Electron Microscopy Sciences, ref. 15710) for 30 min at RT and permeabilized for 1 h on ice with permeabilization buffer (PBS + 5% NGS (Sigma-Aldrich, cat #NS02L) and 0.1% TRITON X (Eurobio, cat #GAUTTR00-01). Cells were then stained with primary antibody resuspended in permeabilization buffer and prepared for analysis, as described for the extracellular markers. Finally, the samples were analysed in a FACSCanto II cytometer with FACSDiva analysis software (Becton Dickinson).

RNA extraction and qPCR analysis
Total RNA extraction was performed using TRIzol reagent (Ambion, cat. #15596-026), according to the manufacturer's instructions. Next, its concentration and quality were measured using a NanoDrop ND-100 spectrophotometer (NanoDrop Technologies), and only samples with 260/280 and 260/230 ratios > 1.8 were accepted. cDNA was synthesized starting from 500 ng of the extracted RNA using iScript Advanced cDNA Synthesis Kit for RT-qPCR (BioRad, cat. #1725038), following the manufacturer's indications. Moreover, to avoid side effects caused by retrotranscription efficiency, we added ERCC RNA Spike-In mix diluted 1:500 to the RNA (Invitrogen, cat. #4456740). The resulting product was diluted 1:10, and 1 μL was used as the template for RT-qPCR. qPCR analysis was completed using SYBR Selected Master Mix (Life Technologies, cat. # 4472908) on a CFX96 thermal cycler (BioRad). All primers used in this study, listed in S2 Table, were designed using Primer 3 software.

Primer efficiency and quality
The efficiency of each pair of primers was tested. Briefly, serial dilutions starting from 100 ng of cDNA (6 dilutions were generated) of cDNA samples were used as templates for qPCR analysis. The standard (STD) function of the CFX Manager™ Software (Biorad, cat. #1845000) was run to generate the linear regression line, and the efficiency of each pair of primers was obtained. To exclude the generation of nonspecific amplicons, qPCR products were run on an 2% agarose gel with dH 2 O. The products were isolated from the gel using the Wizard SV Gel and PCR Clean-Up System (Promega, cat. #A9282) and were sequenced to underline the specificity of the primers used by Sanger sequencing (S1 File).

Best HKG selection
The selection of the best HKG was carried out by focusing on four different methods: ΔCt method (which values the fluctuation of the ΔCt making a comparison between two or more HKGs) [15], geNorm software (www.genorm.cmgg.be; which calculates the stability of each gene through intragroup differences and mean pairwise variation) [16], NormFinder software (www.moma.dk/normfinder-software; which evaluates gene stability using both intragroup and intergroup changes) [17] and BestKeeper software (www.gene-quantification.de/ bestkeeper.html; which ranks the genes in agreement with the standard deviation of their Ct values in correlation with intragroup alterations) [18].
The data were integrated to obtain a final rank, based on the geometric mean, using the RefFinder tool (www.github.com/fulxie/RefFinder).

Statistical analysis
All statistical analyses were performed using Prism7 (GraphPad) on three biological replicates for each experiment. All data passed the Grubb's test excluding the presence of outliers and then were statistically assessed by the non-parametric one-way analysis of variance (ANOVA) test followed by Bonferroni or Dunnett's multiple comparisons post-hoc test, as indicated in the figure legends.

hPO characterization
The hPOs showed a cauliflower-or cyst-like structure when they were embedded in Matrigel, and they could be passaged and maintained in culture for at least 10 passages without showing phenotypic changes (Fig 1A). After culturing for 10 days, the cellular composition was characterized. Flow cytometric analysis showed that the majority of the hPO cells (>95%) were positive for EpCAM, thus recapitulating the epithelial exocrine compartment of the initial pancreatic tissue. Focusing on SOX9 and UEA1, well-known markers for ductal and acinar cells, respectively [9,19,20], we observed a consistent presence of SOX9-positive cells (71.0% ± 10.0%) and a lower percentage of UEA1-positive cells (19.3% ± 5.0%) (Figs 1B and S1).

Expression profile of candidate HKGs
To determine the most reliable reference genes, 12 widely used HKGs, belonging to different functional classes or pathways, were selected to reduce the probability of including co- Representative bright-field images taken after organoid culture of early (< 2 months) and late (> 2 months) passages for 10 days. Scale bar, 500 μm. (B) Immunophenotypic analysis of hPOs at different passages by evaluation of ductal and acinar markers. Data are expressed as the mean ± standard deviation (n = 3/group). One-way analysis of variance followed by the Bonferroni posthoc test for multiple comparisons was used. regulated genes (S1 Table). Agarose gel electrophoresis of the RT-qPCR products showed that all primer pairs amplified products of the predicted size, without evidence of unspecific amplicons in either the cDNA or negative control reaction (Fig 2A). Moreover, all products showed individual dissociation curves at a temperature higher than 75˚C, indicating that they were not primer dimers (Fig 2B). In addition, a standard curve was generated for each HKG using twofold serial dilution of cDNA. The curves showed that the RT-qPCR efficiency of each primer pair ranged between 95% and 105% (S2 Table), which is in agreement with MIQE guidelines [11]. Considering all samples and setting a Ct cut-off > 35 for weak expression, the 12 candidate reference genes exhibited a wide range of expression levels: RNA18S (Ct = 9.73 ± 1) had the highest expression and YWHAZ (Ct = 31.88 ± 1.24) had the lowest expression (Fig 3A). For each reference gene, all Ct values showed that any outliers passed the Grubb's test; thus, no sample was excluded from further analyses. Finally, pairwise comparison between all samples under investigation (all passages analysed) showed a distribution with a correlation coefficient �0.96, indicating high expression similarity (Fig 3B).
The stability ranking for each category was evaluated using four different statistical analyses, NormFinde, geNorm, BestKeeper and ΔCt method, as reported in Table 2. Moreover, since we used various approaches showing differences in the HKG stability, a comprehensive final ranking of the most stable reference genes, based on the geometric mean (geo.mean) of each gene weight generated by the four methods, was obtained using the RefFinder tool. Considering all the passages analysed the best housekeeping gene to normalize a RT-qPCR study on human pancreatic organoids were RPL13A (geo.mean = 1.19), and HPRT (geo. mean = 1.73) whereas the least stable HKG were RNA18S (geo.mean = 10.49) and ACTB (geo. mean = 11.24) (Fig 4).

Impact of HKGs on the expression levels of selected genes
To further evaluate the impact and reliability of the selected reference genes, RT-qPCR analysis was conducted on two specific marker genes for epithelial ductal cells (EpCAM and SOX9), which represent the most enriched cell population in hPO culture. Next, the expression levels were normalized using both the most (RPL13A), the least (ACTB) stable HKGs and exogenous RNA (ERCC RNA Spike-In) as alternative normalization approach, due to their quality and quantity is known [21][22][23]. Moreover, to avoid possible errors related to the use of only one HKG [11], we also normalized the results using the geo.mean of the Ct value of the two best reference genes identified by NormFinder (PPIA and UBC). Finally, we performed the gene expression normalization using GAPDH, the widely used HKG in the literature. After data normalization, we evaluated the modulation of gene expression using hPO P0 as a reference (Fig 5). We observed a stable expression of EpCAM and SOX9 only when our selected HKG (RPL13A) and the geo.mean were used, on the other hand, using ERCC RNA Spike-In or GAPDH to perform the normalization, instability in the expression of the selected ductal genes was observed, and this is in contrast with the protein levels detected in our hPOs (S1 Fig). Moreover, focusing on the less stable HKG (ACTB) an apparent reduction of EpCAM and SOX9 expression was noticed, this behaviour was also in contrast with the detected protein levels (S1 Fig). Taken together, these results demonstrate that the selection and the number of appropriate HKGs to be used are critical parameters to avoid false results and to obtain reliable and significant data by gene expression analysis.

Discussion
Three-dimensional organoid culture systems have emerged as a powerful tool to accurately recapitulate adult stem cell maintenance, differentiation and disease pathophysiology for various organs, including the brain, intestine, liver, pancreas and kidney [24][25][26][27][28]. Recently, our group has described and presented hPOs as a potential source of functional cells for the treatment of type 1 diabetes as they are able to be expanded and cryopreserved without morphological or molecular changes until passage 5 [9]. There is no doubt that organoids have opened remarkable opportunities very rapidly in many fields; however, to envision their use in research and in clinical contexts, from drug efficacy to organ replacement function, human organoids must be fully and extensively characterized. Due to its low cost and ease of operation, RT-qPCR is consistently used as a first step to assess the expression levels of tissue-specific genes [29,30]. Indeed, gene expression quantification is widely considered the best method to confirm or confute the identity of any type of cell or tissue, including the emerging human organoids. The main drawbacks for this "gold standard" method are the imprecision and the high number of technical errors that could affect the reproducibility and the accuracy of the results generated using this technique [11].
To start the molecular characterization, some crucial aspects must be kept in mind. The first one is the appropriate choice of the reference genes to correct for errors in sample quantification as well as sample-to-sample differences in the efficiency of reverse transcription and PCR amplification [31]. The use of non-validated reference genes can lead to erroneous and inconsistent results, which are biologically meaningless. Therefore, it is important that the expression level of the selected HKG is stable and unaffected by the experimental conditions used in the analysis.
Recently, numerous studies focusing on this critical aspect [32,33] have been presented in the literature, but not even one refers to the human organoid context. Indeed, to the best of  our knowledge, there are no previous reports regarding the validation of reference genes for gene expression assays of hPOs. In the current work, we present the first study on hPOs focusing on the identification and validation of a suitable set of HKGs to use in gene expression analysis during short-and long-term culture. To select the best reference gene for hPO transcriptional analysis, we performed a comparison of 12 commonly used HKGs belonging to 5 different functional families: metabolism/cell viability-related genes (GAPDH, HPRT and YWHAZ); a structure/cytoskeleton-related gene (ACTB); protein folding and modification- The relative transcriptional level of the two ductal markers EpCAM and SOX9 in each sample condition considered was normalized with the best and the worst reference gene (RPL13A and ACTB, respectively), the geometric mean (Geo.mean) of two selected HKGs by NormFinder (PPIA and UBC) and one commonly used gene (GAPDH). One-way analysis of variance followed by Dunnett's post-hoc test for multiple comparisons was used. https://doi.org/10.1371/journal.pone.0260902.g005

PLOS ONE
related genes (GUSB, PPIA and UBC); transcription/translation-related genes (EF1A, RNA18S, RPL13A and TBP); and an antigen processing-related gene (B2M). In order to identify the most suitable reference gene, we analysed the expression data using four independent statistical algorithms previously reported: BestKeeper, NormFinder, geNorm and the comparative ΔCt method. All these algorithms gave consistent results, although the most stable genes were not listed in the same order. Next, to provide a comprehensive ranking of the HKGs, we combined all of the data obtained with the four algorithms using the tool RefFinder [34]. This approach can produce precise results by integrating and balancing all of the features of the mentioned methods [13].
To evaluate the transcription levels of the selected HKGs, we first performed a preliminary validation of our primers. Our results showed no dimer formation and an appreciable efficiency (efficiency = 95-105%) combined with an acceptable level of expression (Ct < 35). Altogether, these parameters should be underlined because they play a key role for the validity of this type of study [11].
Based on our statistical analysis, we found the most stable gene to use alone or in combination for normalization of RT-qPCR analysis of hPO isolation and maintenance in culture. Interestingly, to characterize the long-term hPO culture by gene expression analysis, we observed that the most stable HKG was RPL13A, while NormFinder suggested PPIA-UBC as the best couple genes to perform data normalization. This result was in accordance with previous work showing that the mRNA expression levels of RPL13A was stable in a variety of human tissues, including pancreatic cells, making it an excellent reference genes [35][36][37]. On the other hand, ACTB and RNA18S, two of the most-used reference genes, appeared to be less reliable. This is not surprising, since there are many documented reports showing that the mRNA levels of ACTB and RNA18S fluctuate [38]. However, they are still commonly used as reference genes.
Furthermore, our analysis showed that the intragroup and intergroup comparisons could influence the choice of the HKG to be used. Therefore, it is evident that the identification of a stable reference gene for accurate and reproducible RT-qPCR data normalization remains a critical task, particularly when performing gene expression profile analysis using samples from different individuals.

Conclusion
The present study is the first report of a systemic evaluation of potential HKGs suitable for accurate RT-qPCR normalization of hPOs as well as comparative studies with pancreatic tissue. The HKGs were ranked according to their stability in different experimental setups. This work provides a solid basis for future gene expression analysis in the hPO field and can be taken as an example to generally expand the knowledge of this three-dimensional culture system as well as to act as a guide for future studies.