Comprehensive selection of reference genes for quantitative RT-PCR analysis of murine extramedullary hematopoiesis during development

The purpose of this study was to perform a comprehensive evaluation and selection of reference genes for the study of extramedullary hematopoiesis during development and the early post-natal period. A total of six candidate reference genes (ACTB, GAPDH, HPRT1, PPID, TBP, TUBB3) in four organs (heart, liver, spleen, and thymus) over five perinatal time points (Embryonic days 14.5, 16.5, 18.5, Post-natal days 0, 21) were evaluated by quantitative real-time PCR. The expression stability of the candidate reference genes were analyzed using geNorm, NormFinder, Bestkeeper, Delta CT method, and RefFinder software packages. Detailed methodology for isolation of high quality/purity RNA and analysis is presented. Detailed analysis demonstrated that TBP is the best single reference gene for embryonic samples and HPRT1 is the best single reference gene for post-natal and pooled embryonic and post-natal samples. Organ-level analysis demonstrated that HPRT1 was the most suitable reference gene for heart, liver and thymus samples, while TBP was the best candidate for spleen samples. In general, TUBB3 was consistently the least stable gene for normalization. This is the first study to describe a systematic comprehensive selection of reference genes for murine extramedullary hematopoietic tissues over a developmental time course. We provide suggested reference genes for individual tissues and developmental stages and propose that a combination of reference genes affords flexibility in experimental design and analysis.


Introduction
Extramedullary hematopoiesis (EH), the process of blood cell differentiation from hematopoietic stem cells occurring in organs outside of the bone marrow (BM), is a normal component of pre-natal development. However, postnatally and through adulthood, physiologic hematopoiesis takes place primarily in the BM while EH is considered pathologic. Post-natal EH Animals were time mated and the day the vaginal plug was identified defined as embryonic day (E) 0.5. Embryos were isolated from pregnant dams that had been anesthetized with isoflurane and humanely euthanized by cervical dislocation. Postnatal (P) animals were humanely euthanized using isoflurane and cervical dislocation. Tissues (liver, spleen, thymus, heart) were collected, snap frozen in liquid nitrogen and stored at -80˚C. Three biological replicates of each tissue type were collected at five developmental stages: E14.5, E16.5, E18.5, P0, and P21.

RNA isolation, quantification, purity and integrity
Total RNA was isolated from frozen tissue using silica-membrane based kits from Qiagen. RNeasy Micro kit (Qiagen, Hilden, Germany) was used for all embryonic ages and P0, except for liver samples for which we used the RNeasy Mini kit in all ages. P21 samples were isolated using RNeasy Mini kit (Qiagen, Hilden, Germany). Beta-mercaptoethanol (BME) was added to the lysis buffer to prevent degradation by RNases. Fully automated RNA isolation (QIAcube, Qiagen, Hilden, Germany) was used for samples at E18.5, P0, and P21. E14.5 and E16.5 samples were processed manually. On-column DNAse digestion was also performed for all samples. RNA concentration and purity (260/280nm ratio) was determined using Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE). RNA integrity was determined by capillary electrophoresis (Bioanalyzer 2100, Agilent Technologies, Waldbronn, Germany).

cDNA synthesis and pre-amplification
The cDNA synthesis was carried out in the presence of ProtoScript II First Strand cDNA synthesis kit using 150 ng RNA, 1 μl of oligo d(T)23 VN (50 μM) and 1 μl of random primer mix (60 μM) in a final volume of 20 μl according manufacturer's instructions (New England Bio-Labs, Ipswich, MA). To pre-amplify small amounts of cDNA we used Taqman PreAmp Master Mix (Applied Biosystems, Foster City, CA). The PreAmp mix was performed using a primer pool consisting of 50 nM final concentration of each primer, cDNA at a final concentration of 22.5 ng in a final volume of 10 μl. The pre-amplification reaction was conducted as follows: 1 cycle of 95˚C for 10 minutes and 14 cycles of 95˚C for 15 seconds and 60˚C for 4 minutes. After pre-amplification samples were stored at -20˚C.

Candidate reference genes
Six reference genes with distinct functions in live cells were selected for evaluation of their expression profile (Table 1): Beta-actin (ACTB), Glyceraldehyde 3-phosphate dehydrogenase (GAPDH), Hypoxanthine-guanine phosphoribosyltransferase (HPRT1), Peptidylprolyl isomerase D (PPID), TATA box binding protein (TBP), and Tubulin beta 5 (TUBB3). The NCBI reference number of the gene used to design the primers, primer sequence, Universal ProbeLibrary (UPL) (Roche) probe number are detailed in Table 1. The amplicon size was kept in a narrow size range (61-123 bp) to ensure similar PCR efficiencies for comparison between assays. The UPL probes were labeled at the 5' end with fluorescein (FAM) and at the 3' end with a dark quencher dye. Locked Nucleic Acids (LNA) incorporated into the sequence of each UPL probe to increase specificity.

Primer/Probe design and real-time PCR
We used ProbeFinder (Universal ProbeLibrary Assay Design Center, Roche, Indianapolis, IN) for mouse species to design primers and UPL probes for six reference genes (Table 1). Real-time PCR reactions were prepared in a final volume of 10 μl containing 0.1 μl of UPL probe (10 μM), 0.2 μl of left (10 μM) and right (10 μM) primers, 5ul of 2X Kapa probe fast qPCR master mix (Kapa Biosystems, Wilmington, MA), 2 μl of cDNA and 2.7 μl molecular grade water. Three biological replicates for each sample were run in a 384-well white plate optimized for Roche 480 Light Cycler (Phenix Research Products, Candler, NC). Non-template controls were run in triplicate. Interplate calibrators (IPC) were also included in triplicates in each plate as a method for inter-run variation compensation (Tataa Biocenter AB, Goteborg, Sweden). All plates were run in Roche LightCycler 480 machine using the Mono Color Hydrolysis Probe/UPL program: 1 cycle of pre-incubation step at 95˚C for 5 minutes, 45 cycles of denaturation step at 95˚C for 10 seconds, 60˚C for 30 seconds, and 70˚C for 10 seconds) with acquisition mode "single" at the last step of amplification. qPCR efficiency RNA isolation and cDNA synthesis were processed independently for each sample. Three independent cDNAs were prepared for each sample. cDNA from each organ was pooled and three pools of cDNA were obtained for determining qPCR efficiency. All the pooled specimens were tested using five serial dilutions: 50, 5, 0.5, 0.05, and 0.005 ng (final concentration of cDNA). Efficiency of amplification for each gene was determined by analysis of the slope of the standard curve (Light Cycler 480 SW 1.5.1, Roche Applied Science, Penzberg, Germany). The highest possible quality PCR reaction produces a standard curve with an efficiency of 2, which refers to the number of target nuclei acid doubles with each amplification cycle. The error is the mean squared error of the single data points fit to the regression line. An error < 0.2 was considered acceptable.
In the case of geNorm analysis, stability of reference genes is determined by calculating the M value (average pairwise variation analysis). The most stable gene yields the lowest M value and the two most stable reference genes are determined by stepwise exclusion of the least stable gene. In NormFinder software, the stability value of the candidate genes is calculated according to their intra-and intergroup variations. Stably expressed genes will have low stability values.
BestKeeper utilizes a two-way ranking: Pearson's correlation coefficient and BestKeeper computed SD values. Gene ranking by Delta Ct method is calculated by pair-wise comparisons.
Using raw Ct values, the average SD of each gene is inversely proportional to gene stability. The more stable gene has the lowest CV and SD. RefFinder integrates all four programs (geNorm, NormFinder, BestKeeper, and Delta Ct method) and calculates a weighted geometric mean for a ranking of the best reference genes. Means of different groups were compared and analyzed using the Student's t-test. Standard error (SEM) was used. GraphPad Prism version 7.0b (GraphPad Software) was used for the box plot, statistical analysis and graphs plotting.

Quality of RNA samples
Isolation of quality of RNA including integrity and purity is crucial for comparative gene expression experiments. We used a micro-capillary electrophoretic RNA separation to estimate the RNA integrity number (RIN) and the ratio of 28S:18S ribosomal RNA by Agilent 2100 Bioanalyzer [30]. RIN values range from 7.8 (spleen) to 10 (heart, liver and thymus) ( Fig  1A). Because splenic samples have high nuclease and nucleic acid content, BME was added to the lysis buffer during RNA isolation to minimize RNA degradation. Splenic sample RIN without BME ranged from 2.6-3.

PCR performance and efficiency
Experiments involving embryonic organs often rely on small amounts of source tissue. To overcome this limitation, we utilized a pre-amplification kit (PreAmp, Applied Biosystems) to amplify small amounts of cDNA into many real-time PCR reactions, thereby allowing testing of all selected reference genes. An initial test was done in triplicates using P0 samples (heart, liver, spleen, and thymus) and HPRT1 reference. We compared ten serial dilutions (original cDNA, PreAmp pure cDNA, and PreAmp cDNA dilutions 1:5, 1:20, 1:50, 1:100, 1:200 1:400, 1:800, and 1:1,600) to confirm that there is no variation in the expression results introduced by the pre-amplification step and to determine which is the best dilution to use for further experiments (data not shown). We selected the PreAmp 1:100 cDNA dilution for further analysis. Using a standard curve generated by serial dilutions of cDNA from a pool of embryonic tissue and a pool of post-natal tissue, efficiency was calculated for each of the 6 candidate genes. PCR efficiency ranged from 1.9 to 2.2 ( Table 2), indicating good to excellent efficiency for all 6 primers.

Reference gene expression stability over time
In order to identify the reference gene(s) with the lowest variation in expression among different tissues and different developmental time points, raw Cp values of the 6 reference genes were plotted for embryonic samples (Fig 2A), post-natal samples ( Fig 2B) and pooled embryonic and post-natal samples ( Fig 2C). All three plots demonstrate a similar trend, with the highest variation in TUBB3 compared with the other genes.
In order to detail the temporal variability in our selection of reference genes, we analyzed the genes based on the ages E14.5, E16.5, E18.5, P0, and P21 using pooled cDNA from the four organs of interest (heart, liver, spleen, thymus) at those time points (S1 Fig). We found that expression of each individual reference gene is near-constant across the developmental stages, validating the grouping of time points into embryonic, post-natal, and pooled. Additionally, ACTB, GAPDH, HPRT1, and TBP demonstrate very similar levels of expression to one another over time. However, TUBB3 and PPID differ from the other candidate reference genes, demonstrating consistently higher gene expression over the time course.
In order to determine the appropriate number of reference genes to utilize for accurate normalization across tissues and time points, we used the geNorm and NormFinder modules of   GenEx (Fig 3). geNorm ranks reference genes based on their expression stability (M), where a low value is the more stable than a high value (Fig 3A-3C). The algorithm uses paired expression values to sequentially eliminate the gene that shows the highest variation relative to the other genes, resulting in a pair of recommended reference genes. As shown in the Fig 3, the order of gene stability ranking is similar for embryonic and post-natal samples. HPRT1 and GAPDH (M = 0.81 for embryonic; M = 0.84 for post-natal) were the most stable pair of reference genes (Fig 3A and 3B). However, in the case of pooled samples the most stable pair of genes were HPRT1 and PPID (M = 0.93, Fig 3C). NormFinder was used to calculate the expression stability value (standard deviation) for the six candidate reference genes. NormFinder identified TBP as the single best candidate gene for embryonic (SD = 0.52), post-natal (0.43), and pooled specimens (0.10) (Fig 3D-3F). NormFinder also calculates the Accumulated Standard Deviation (Acc SD) in order to determine the optimal number of reference genes (Fig 3). When we analyzed embryonic or postnatal samples, the Acc SD indicated the lowest values for a combination of three reference genes for embryonic (TBP, GAPDH, HPRT1, Acc SD = 0.4477) and post-natal (TBP, HPRT1, A comprehensive ranking of the six candidate reference genes for embryonic, post-natal and pooled samples were calculated using geNorm, NormFinder, BestKeeper, Delta Ct, and RefFinder ( Table 3). The results indicate that HPRT1 and/or TBP were the most frequently identified single best gene for embryonic samples. HPRT1 and/or PPID were the most frequent single best gene identified in post-natal samples. HPRT1 and/or TBP were the most frequently identified best reference genes for pooled samples. RefFinder analysis involves an integrated calculation including geNorm, NormFinder, BestKeeper and Delta Ct methods. This analysis demonstrated that TBP is the best reference gene for embryonic samples and HPRT1 is the best reference gene for post-natal and pooled samples.

Gene expression analysis by organ
In order to understand tissue level variability in our selection of reference genes, we analyzed the genes based on the four organs of interest (heart, liver, spleen, thymus) (Fig 4). By using a pool of cDNA from each of the developmental stages for each organ, we found significant differences in expression within organs in all of the reference genes except PPID, for which there were no significant differences between organs. TUBB3 demonstrated high degrees of variation between organs.
The results of NormFinder, geNorm, DeltaCT, and RefFinder demonstrate consensus, indicating that HPRT1 is the most stable reference gene in liver and thymus (Table 4). However, BestKeeper indicates HPRT1 was the second most stable gene for liver and thymus. In heart, HPRT1 is the most stable reference gene according geNorm, BestKeeper, and RefFinder. While NormFinder results show ACTB as the most stable gene in heart, Delta CT methodology demonstrates that TBP is the most stable gene. The spleen has either GAPDH (NormFinder, geNorm) or TBP (geNorm, BestKeeper, Delta CT, and RefFinder) as the best reference gene. TUBB3 was consistently the least stable or the second least stable in all organs. By using RefFinder we were able to identify consensus reference genes between the 4 different algorithms. Based on these analyses, HPRT1 was the most suitable reference gene for heart, liver and thymus samples, while TBP was the best candidate for spleen samples (Table 4).

Discussion
For reliable and reproducible qPCR analyses, selection of the correct reference gene(s) for a given experimental setup is of paramount importance. EH occurs under both physiologic and pathological conditions. EH plays an essential role during normal fetal development, is a normal response to infection and inflammation, can accompany malignant conditions such as myelofibrosis, and can also occur during chronic inflammation [17]. No previous studies have identified potential reference genes in hematopoietic organs at different stage. In this study, we have described a step-by-step procedure for selection of reference genes to study EH. Genes were selected based on stability analysis in multiple tissues and over multiple developmental time points. We found that the combination of HPRT1 and TBP were the best pair of reference genes for normalization of gene expression in EH studies involving only embryonic or both embryonic and post-natal time points. For studies involving only post-natal time points, the combination of HPRT1 and PPID appear to be the optimal combination. RNA quality can vary depending on the complexity of the sample, making the choice of RNA extraction methodology a critical component of experimental design. To establish an optimal protocol for the diverse collection of EH organs, we initially compared TRIzol (Invitrogen) (data not shown) and a column-based RNeasy micro kit (Qiagen). TRIzol is a monophasic solution of phenol and guanidinium isothiocyanate that simultaneously solubilizes qPCR reference gene selection for extramedullary hematopoiesis biological material and denatures protein. The subsequent addition of chloroform causes phase separation where protein RNA is extracted. Other have reported that RNA isolated with TRIzol is often contaminated with genomic DNA, significantly impairing interpretation of RT-PCR results [31]. We also noted poor RNA quality (RIN as low as 1) when using TRIzol (data not shown). The RNeasy micro kit (Qiagen) technology uses selective binding properties of a silica-based membrane to maximize quality of isolated RNA. Using this approach, we noted maximal, or near-maximal RIN with a clean electropherogram for all four organs ( Fig  1A and 1B). Other authors have reported difficulties in isolating RNA from certain tissues, including spleen, because of high levels of endogenous nucleases [32]. To overcome this limitation, we utilized BME to prevent RNA degradation in the lysis buffer in all samples. With these approaches, we were able to obtain highly pure RNA (A 260 /A 280 1.85-2.25) for further analysis. Finally, because of our interest in embryonic as well as post-natal time points, we were faced with the challenge of small sample amounts. To overcome this limitation, we utilized preamplification of cDNA with target specific primer pairs.
Because the embryonic period spans major organogenesis and maturation of tissues and organs during rapid growth of the body, embryonic samples represent potential sources of high variation. Mamo et al [33] showed moderate changes for the genes ACTB, GAPDH,  HPRT1 and others in an in vitro and in vivo experiment at different stages using whole embryos. The cells comprising the embryo have a vastly heterogeneous nature, which leads to greater variation in the endogenous biological processes, and greater variation in the sensitivity of cells to their environment [34]. Because these factors can introduce intra-and inter assay variations, we elected to isolate the specific EH organs of interest (heart, liver, spleen, thymus) and analyze gene expression in these organs independently instead of analyzing whole embryo homogenates, as has been previously described [35]. Accordingly, when analyzing variation in Cp, we noted little variability in post-natal samples as compared to embryonic samples (Fig 2). Further, individual genes showed varying degrees of variability, with TUBB3 highest, and ACTB and GAPDH demonstrating the least variability. The analysis of reference genes for studying EH was aided by the use of free and commercially available software packages (NormFinder, geNorm, BestKeeper, DeltaCT, RefFinder) designed for reference gene validation. These software packages differed slightly in the selection of best reference gene or combination of reference genes, emphasizing the importance of understanding the software package's approach and limitations. RefFinder is a package that integrates the other algorithms and provides a comprenhensive analysis by assigning a weight to each individual gene and calculating the geometric mean of these weights in order to generate a final overall ranking. RefFinder identified TBP as the most stable reference gene for embryonic tissue, which echoed the results from NormFinder and the Delta CT. For postnatal and pooled samples (embryonic and post-natal altogether) the most stable reference gene identified by RefFinder was HPRT1. Similar results were obtained with geNorm, BestKeeper, and Delta CT method for post-natal, and NormFinder and BestKeeper for pooled samples. Based on these analyses, we propose that the maximal flexibility for studying EH in all of these tissues/times is afforded by using a set of three reference genes. We have identified the best trio of reference genes to allow flexibility of analyses between embryonic and post-natal tissues. HPRT1/TBP/GAPDH were the best genes in combination for embryonic samples, HPRT1/ PPID/TBP were the best combination for post-natal samples and HPRT1/TBP/GAPDH were the best trio for all pooled samples.
A limitation of the current study is the use of whole organs rather than specific cells to identify appropriate reference genes through the embryonic and post-natal periods. Indeed, there is a growing literature on the importance of bi-directional cross-talk between hematopoietic stem cells and their surrounding niche cells in both normal myeloid and lymphoid development, as well as in the development of malignancy (36-38). An extension of the current approach would be the use of cell sorting, through either magnetic bead or flow cytometrybased immuno-selection, to separate hematopoietic cells from the surrounding mesenchymal environment. These separate compartments could then be subjected to similar analyses to identify reference genes specific for those cell types of further interest. However, the reference genes we have identified here should prove useful in future studies directed at specific cellniche interactions in each of these disease processes.

Conclusions
This is the first study to detail the comprehensive selection of reference genes for normalization of expression in murine extramedullary hematopoietic tissues over multiple developmental stages. We conclude that there are optimal reference genes for individual tissues and developmental stages, but selection of a combination of reference genes affords flexibility in experimental design and analysis. Our study provides the basis for the selection of reference genes for future gene expression studies related to extramedullary hematopoiesis.