Improving characterisation of human Multipotent Stromal Cells cultured in 2D and 3D: Design and evaluation of primer sets for accurate gene expression normalisation

Human Multipotent Stromal Cells (MSCs) are a valuable resource for regenerative medicine and are widely studied. They can be isolated from a variety of tissues and differentiate into multiple cell types (multi-potent). Many reports have been published using human MSCs and to be able to compare outcome, or be able to identify differences between MSCs, several cell surface markers have been proposed. Nevertheless, still many differences remain. Gene expression is known to be different between cell stage and origin. Furthermore, cells cultured on a culture dish (2D) show different gene expression profiles as compared to cells grown on scaffolds (3D). Even the RNA extraction method and the selection of genes used for normalisation have a role in gene expression profiling. To be able to compare gene expression data from samples cultured in different dimensions and RNA extracted using a variety of protocols we set out to define a set of reference genes suitable to normalise qPCR data from a very heterogeneous sample set. Hereto, Trizol was used to extract RNA from human MSCs cultured in 3D and 2D to validate newly designed and previously published primer sets. Subsequently, RNA from fresh human MSC samples and samples stored in RLT-buffer, Trizol or RNAlater was extracted using RNeasy and Trizol methods. All samples have been used to rank the candidate reference genes according to their stability after qPCR enabling identification of the most suitable reference gene(s) for normalisation of a heterogeneous sample set. The most stably expressed reference genes indicated superior normalisation of MSC marker gene expression over the least stable reference genes.


Introduction
Multipotent stromal cells, also known as mesenchymal stem cells, mesenchymal stromal cells or medicinal signalling cells [1] (MSCs), hold great promise for regenerative medicine and tissue engineering. MSCs can be isolated from virtually all post-natal organs and tissues in the human body [2], are plastic adherent whilst maintaining multipotency when expanded, can differentiate in vitro into several distinctive end-stage cell types and possess immunoregulatory and regenerative capacity through secretion of bioactive macromolecules [3]. Because of their relatively easy isolation from many tissues the identity of the true MSC is unknown. The definition of the MSC phenotype is focussed on the expression of cellular surface markers (CD105, CD90 and CD73) or the lack thereof (e.g. CD45, CD34, CD19 and CD11b) [4]. These minimal requirements resulted in numerous reports using MSCs with variable outcome regarding maintaining multipotency or differentiation potency and efficiency. Other reasons why differences between MSCs can be observed are the age of the donor [5] and the duration of culture of the MSCs [6].
To properly perform MSC batch analysis, gene expression profiling could be utilised and quantitative PCR (qPCR) is the gold standard. Low amounts of starting material are required for qPCR and it has a very wide dynamic range. This technique is very useful if the target gene is known. If it is unknown which gene(s) will be differentially expressed, microarray might be the preferred method, though its dynamic range is very low. An alternative to microarray is RNA sequencing (RNAseq). Although a bit more expensive than microarray, the much broader dynamic range, possible discovery of new genes or splice variants and the ability to examine DNA variations are only some of the benefits [7]. These techniques require the extraction of RNA from MSCs cultured in 2D (on culture dish in monolayer) or 3D (cultured in a scaffold of any particular material). Several methods have been developed to extract RNA from cells and are available through even more suppliers. Usually, at least 100ng good quality RNA is required per sample to perform RNAseq. For lower amounts or significantly degraded RNA some solutions are available like RNA amplification kits or kits developed for low quality RNA, respectively. Nevertheless, as these kits can introduce biases, their use should ideally be avoided [8]. Microarray or RNAseq results should always be confirmed by qPCR and therefore there needs to be additional RNA remaining from the samples that can be used for microarray or RNAseq. In general, 1μg RNA is used to obtain 20μl cDNA which is sufficient for the analysis of 4-10 genes by qPCR. To obtain high yields and the best quality of RNA for RNAseq and qPCR, protocols for the extraction of RNA need to be optimised depending on the culture condition and possibly on the transportation or storage method.
Quantitative PCR data need to be normalised using reference genes before accurate gene expression analysis can be performed. The candidate reference gene (RG) expression should not vary between samples under investigation or in response to experimental treatment [9]. Two of the most frequently reported RGs for gene expression data normalisation are betaactin (ACTB) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH). Recently, reports advised to forgo these genes identifying them as unsuitable for accurate gene expression normalisation [10,11], possibly due to the inability to design gene specific primers allowing to discriminate between the actual gene and their numerous (>60) pseudogenes [12]. Several other genes have been reported as normalisers in MSC gene expression analysis. Unfortunately, neither does a universal RG exist nor is there a fixed number of defined RGs to use and these should be determined experimentally [13]. In this report we aimed to identify suitable RGs for gene expression normalisation regardless of cell culture dimension, isolation, storage, or RNA extraction method. RNA has been extracted from two human MSC lines cultured in 2D and 3D by two different methods (RNeasy and Trizol) after storage in RLT-buffer, Trizol or RNAlater and evaluated for quality (gel electrophoresis) and yield (Nanodrop). Thirteen primer sets have been designed and evaluated as candidate RGs using four different ranking approaches. Most rankings indicated TBP amongst the most suitable of all candidate RGs for the normalisation of gene expression data from heterogeneous hMSC sample sets. Nevertheless, none of the rankings were the same nor could one single gene be nominated as suitable for normalisation of all comparisons. Therefore, as advocated already by others [9,14], we strongly recommend to use multiple candidate RGs and all or a selection of the validated primer sets described in this report could be suitable for future quantitative PCR studies involving human MSCs cultured in 2D and/or 3D.

Materials and methods
All reagents were supplied by ThermoFisher Scientific (Hemel Hempstead, UK) unless otherwise specified.

Cell culture
All cell culture incubations were done in a humidified 5% CO2 incubator at 37˚C.
Human bone marrow derived MSCs with enforced TERT expression. MSC-hTERT cells [15] were seeded at 0.5 � 10 6 cells/well in DMEM + 10% FBS + 1% p/s (complete DMEM) in a 48-well culture dish or a 12-well culture dish. The fibrinogen-alginate and fibrinogen scaffolds were placed in a 48-well culture dish whereas the polycaprolactone-poly[N-isopropylacrylamide](PCL-PNIPAAM) beads were placed in a 12-well culture dish covering the entire bottom of the well before 0.5 � 10 6 cells were applied to the scaffolds and allowed to seed for 3h in the incubator. All wells were incubated for 48h in complete DMEM before performing Alamar Blue staining. After replacing the Alamar Blue staining solution, cells were cultured for another 16h in complete DMEM before adding Trizol and subsequent RNA extraction.
Human umbilical cord derived MSCs. This study was approved by the National Research Ethics Service South Central-Oxford C Research Ethics Committee (REC reference number 09/ H0606/5+5) and Oxford Radcliffe Biobank (reference number 16/A052). Written informed consent was obtained from all of the participants. Human umbilical cord was obtained from full-term (>36 weeks) caesarean section deliveries or vaginal birth with informed consent of the mother in the John Radcliffe Hospital. A portion of umbilical cord (length >20 cm) was placed into a sterile container (250 ml pots, Sterilin UK) with 100ml tissue collecting solution (PBS + 1% p/s). The collection container was kept at 4˚C for storage and brought to the laboratory for processing within 24 hours. The umbilical cord segments were dissected after washing twice in PBS. Wharton's Jelly tissue of umbilical cord was carefully separated from cord lining, vein and arteries first, and subsequently minced into 1 mm 3 fragments. The minced fragments were placed into a tissue culture dish and left for attaching without adding culture media for 30 minutes in the incubator. UC-MSC culture media composed of DMEM, 10% foetal bovine serum (FBS), 1% Glutamine and 2ng/ml basic fibroblast growth factor (bFGF; Peprotech, London, UK) was gently added to the culture dish. Media was changed once a week until cells migrating out from the fragments was detected. When outgrown cells reached 70-90% confluency, they were harvested by Trypsin dissociation for further culture or flow cytometry. After culture for 4 days in UC-MSC medium, cells were transferred to three 12-well culture plates (1ml per well) at a concentration of 0.1 � 10 6 cells/ml. After an additional culture for 3 days in fresh medium the cells were washed with PBS. Culture plates 2 and 3 were treated with 200μl Trypsin per well for 3 min at 37˚C before adding 800μl PBS. The total volume was transferred to tubes and centrifuged for 5 min at 13,000 rpm and subsequently 1ml of either Trizol, RLT-buffer (Qiagen, Manchester, UK) or RNAlater was added to the pallet. RNAlater samples were stored for 24h at 4˚C before transferring to either -20˚C or -80˚C. Cells from plate 1 were treated with 1ml of either Trizol or RLT-buffer (Qiagen). One sample each was used for immediate RNA extraction whereas the other samples were stored at either -20˚C or -80˚C until RNA extraction.

Phenotype characterisation of UC-MSCs by flow cytometry
Fluorescent marker conjugated antibodies CD105, CD90, CD73, CD34, and CD31 (R&D Systems, Minneapolis, MN, USA) were used to confirm the phenotypic characteristics and to verify the lack of contaminating cells, according to the markers proposed by the International Society for Cellular Therapy [4]. Non-specific background signal was measured using an isotype control cocktail consisting of mouse anti-human IgG2A, IgG2B and IgG1 (R&D Systems).
Cells were harvested and suspended in flow cytometry buffer (PBS containing 0.5% (w:v) BSA) at a concentration of 4 × 10 6 cells/ml. Every 2x10 5 cells were incubated with 10μl of fluorescent marker conjugated antibody at 4˚C for 30 min in the dark. The cells were then washed and finally resuspended in 400μl of flow cytometry buffer for analysis. To determine the viability, 1μl Propidium iodide (PI; Sigma-Aldrich, Dorset, UK) was added to each tube just before the sample was analysed using BD FACS Canto-F60 cell sorter (BD Biosciences, Wokingham, UK). Data was analysed by Flowjo software (FlowJo LLC, Ashland, OR, USA).

Alamar blue staining
Fibrinogen-alginate and fibrinogen scaffolds were transferred to a new well. Cells cultured in 2D and on beads were washed with PBS once. Alamar Blue stain (10% v/v) was added to the wells and incubated for 3h at 37˚C before absorbance was measured at 570nm on a Spectra-Max i3x Multi-Mode Microplate Reader (Molecular Devices, San Jose, CA, USA). Blank wells without cells were used as negative controls.

RNA extraction
Samples stored in RNAlater were centrifuged and the supernatant removed before adding either Trizol or RLT-buffer (Qiagen). Cell lysates in Trizol were subjected to RNA extraction as described previously [16]. RNA from samples in RLT-buffer were extracted using the RNeasy kit per manufacturer's protocol. RNA quantity was determined using a Nanodrop One spectrophotometer. RNA quality was assessed by gel electrophoresis using approximately 0.5μg RNA per sample in a 1% agarose in 1% TAE buffer gel containing 3.5μl SybrSafe per 100 ml gel.

cDNA generation and standard curve composition
The QuantiTect Reverse Transcription kit (Qiagen) was used to convert 1-4 μg RNA into cDNA per manufacturer's instructions. Reaction volumes were supplemented with RNAse/ DNAse-free water to a final volume of 100μl. For the standard curve to evaluate PCR efficiency, S1 was generated with a pool of 5μg RNA converted into cDNA from all PCL-PNIPAAM beads and 2D 12-wells samples (n = 6) resulting in a final volume of 600μl. Subsequently a 4-fold dilution series was prepared starting with 150μl S1 in 450μl RNAse/DNAse-free water until S5. A no template control (NTC) consisted only of RNAse/DNAse-free water. All cDNA samples were stored at -20˚C until processed by qPCR. cDNA was used for qPCR without concentration measurements.

Quantitative PCR
Primers were designed (Primer-BLAST) if not available from previous papers or not within set ranges and compliant with in silico testing. All primers have been assessed for specificity (BLAST), dimers and hairpins (OligoAnalyzer 3.1; Tm of hairpin <60C, dimers ΔG5' < 5 kcal/mol; ΔG3' < 3 kcal/mol) and splice variant detection if applicable (Primer-BLAST). All amplicons have been screened for secondary structures [mFOLD; [17]]. Full details according to MIQE guidelines [13] can be found in Table 1.
Quantitative PCR was based on intercalating dye technology using SyGreen and a passive reference ROX (PCRBiosystems, London, UK) on an Applied Biosystems StepOnePlus Real-Time PCR system according to the manufacturer's instructions in 96-well plates and plate seals. Primers (Sigma-Aldrich) had a final concentration of 400nM each. Standards, NTCs and samples were all done in duplicate. For each reaction, 2μl cDNA template was used in a 20μl final reaction volume containing 3mM Mg 2+ and no Na + . Reactions started with 3 min at 95˚C, followed by 40 cycles of 15 s at 95˚C and 30 s at Tm. This reaction was followed by a melting curve, stepwise increasing temperature each 15 s by 0.5˚C, ranging from 65˚C to 95˚C. Optimal Tm was determined using a temperature gradient ranging from 56˚C to 66˚C on a 4-fold dilution series using cDNA derived from all PCL-PNIPAAM beads and 2D 12-wells samples. Baseline was set automatically and well specific during analysis. All samples for a particular gene were in the same plate and all plates contained the standard curve.

Statistical analysis
All statistical analyses were performed using GraphPad Prism 8. Raw Cq values obtained from StepOne software were used for validating the primer sets and ranking their stability using the web-based free analysis program RefFinder [18] integrating BestKeeper [19], the delta Ct method [20], GeNorm [9], and NormFinder [21]. Cited references describe the statistics and manipulations in detail. The ΔΔCq values for the target genes (i.e. NT5E, THY1, ENG) were normalised against different reference genes (or a combination of reference genes). Only samples with valid ΔΔCq levels across all studied sets were used for grouped correlation analysis in MATLAB. The reference gene sets that were used to normalise the target genes were either TBP, YWHAZ, and RPL13A (most stably expressed), PPIA, HPRT1, and B2M (least stably expressed), the two most stably expressed genes (TBP and YWHAZ), or a single reference gene (TBP or YWHAZ). The Pearson coefficient (R) of the target genes ΔΔCq values was calculated for each sample against the remaining sample set. The R for each sample against the samples from different groups (i.e., 2D UC-MSCs, 2D MSC-hTERT, or 3D MSC-hTERT) were averaged as the relative correlation against that particular group of samples and subsequently referenced as the grouped-correlation coefficient. The grouped-correlation coefficient for each sample was plotted on a 3-dimensional scatter plot.

Results and discussion
MSCs, immortalized by transduction with human telomerase reverse transcriptase (hTERT), referred to as MSC-hTERT [15] and MSCs derived from umbilical cord (UC) were cultured to obtain a variety of RNA samples. To generate these UC derived MSCs (UC-MSCs), fresh umbilical cord has been obtained and cells were able to adhere to plastic culture flasks. Flow cytometry indicates high levels of typical MSC markers CD73, CD90, CD105 and very low levels of CD31 and CD34 (Fig 1) indicating a genuine MSC phenotype [4]. RNA samples from MSC-hTERT and UC-MSC cultures were obtained using different extraction protocols (RNeasy and Trizol) after cell storage in different media (RLT, Trizol, and RNAlater) for several time periods (0 days, 7 days, and 30 days) to identify candidate RGs for accurate normalisation of gene expression in human MSCs cultured in 2D and 3D (Table 2 and Fig 2A). This variability in cell isolation, storage and RNA extraction method resulted in very similar RNA yields and quality ( Table 2). One sample ( Table 2, III.B2) indicated no or very little RNA and was excluded for further analysis. Although the A260/A280 ratios for the MSC-hTERT samples were somewhat lower, RNA yield was proportional to the starting cell amount when compared with the UC-MSCs.
Viability in cell culture can be assessed using Alamar Blue. The viability of MSC-hTERT samples cultured in 2D in 12-well dishes was comparable to cells cultured in 48-well dishes (Fig 2A). Whereas both types of fibrin scaffolds were transferred to a new plate before Alamar Blue assay, the PCL-PNIPAAM beads were not. Therefore the Alamar Blue level and the RNA yield (Fig 2A) might also have come from cells not attached to the beads but to the surface of the culture dish. The fibrin and fibrin-alginate scaffolds are porous and Trizol might not be able to lyse all cells or be aspirated as much as from the beads or the 2D culture dishes which might explain a reduced RNA yield from the fibrin-alginate scaffold when compared to the beads and 2D cultured cells. Although the fibrin scaffold showed low levels of Alamar Blue staining the RNA yield seemed relatively higher than from the others. Unfortunately the electrophoresis indicated lower amounts of RNA compared to the other samples ( Fig 2B) and indeed the Nanodrop readings might have been skewed possibly because of contaminants such as phenol.
The remaining samples (N = 42) were all very different in cell type, culture dimension, storage and extraction method. Because of this variable sample composition 13 genes (ACTB, B2M, GAPDH, GUSB, HMBS, HPRT1, PPIA, PUM1, RPL13A, RPLP0, TBP, TFRC, and YWHAZ) were selected to develop primers and evaluate their suitability as reference genes for the accurate normalisation of gene expression data regardless of prior experimental differences. These genes were selected because of their use in many other reports. However, to make sure that the primers met the in silico design requirements (e.g. amplicon size 70 -250bp, secondary structure energy >-3.0 for primer annealing sites on amplicon, primers separated by at least one intron, etc.) we designed those for which we could not identify previously reported appropriate primer sets (Table 1). Primers were validated by qPCR for their efficiency and optimal annealing temperature (Tm) using a standard curve consisting of a 4-fold dilution series composed of cDNA converted from RNA extracted from cells grown on 12-well culture dishes and beads. From the standard curve plot the slope can be determined and efficiency (Eff.) was calculated according to the equation; Some genes were more abundantly expressed than others (Fig 3A) though there was no correlation between amplicon size and efficiency (Fig 3C; r 2 = 0.007705, p = 0.7755) or quantification cycle (Cq, a.k.a. Ct or Cp) (Fig 3D; r 2 = 0.03130, p = 0.5631). Also the primer sets with potential annealing difficulties because of possible folding of the amplicon at the annealing site performed equally well as the others (Table 1 and Fig 3B, 3C and 3D). Despite careful design, some primer sets were able to form primer-dimers causing an amplification signal in the  NTCs. No primer-dimers were detected in the samples. Nevertheless, the presence of primerdimers always needs to be evaluated when performing qPCR since they could have had an effect on the overall amplification signal. All primer sets but HMBS showed efficiencies between 91-107% (Table 3 and Fig 3) and the remaining primer sets were used to evaluate  expression levels in the available samples to determine their stability. Several samples (n = 13) showed undetectable levels of gene expression when analyzing B2M, HPRT1, PUM1, RPL13A, and TFRC and were excluded for further analysis (Fig 4A). The remaining samples (n = 29) were used to rank these genes (n = 12) according to their stability using the currently available major computational programs GeNorm [ [9]; Fig 4B] [18] resulting in an overall comprehensive ranking ( Fig 4F).
All but BestKeeper ranked TBP as the best reference gene and YWHAZ in the top 3. RPL13A was found in the top 5 of all rankings whereas PUM1 was ranked second using Normfinder and the Delta Ct method though amongst the lesser stable genes according to Genorm and BestKeeper. The most commonly used genes for normalisation are ACTB and GAPDH which are found to be moderately stable except for GAPDH ranking third when using Genorm. BestKepper ranked GUSB first whereas all the others ranked this gene as the  (Table 1) are depicted as open triangles. https://doi.org/10.1371/journal.pone.0209772.g003 Reference genes for hMSC qPCR least stable. Also, HPRT1 was ranked amongst the least stable genes in all but BestKeeper ranking.
As pseudogenes can skew the gene expression of the genuine gene their number in the genome should ideally be absent. The NCBI Gene databank (www.ncbi.nlm.nih.gov/gene) indicates that only B2M, PUM1, TBP, and TFRC are devoid of any known pseudogenes whereas the others, and GAPDH and PPIA in particular, have some or numerous known pseudogenes present in the genome (Table 1). According to Liao and colleagues the actual number of pseudogenes for human ACTB is even higher at 64 pseudogenes [12].
As some of the samples were from cultured MSC-hTERT and other samples from hUC-MSC we established a comprehensive ranking for these samples separately. Comparing these rankings (2D MSC-hTERT vs. 2D UC-MSC), striking differences could be observed ( Fig  4G and 4H). Not only was YWHAZ ranked second in UC-MSC whilst ranked second last for MSC-hTERT, also B2M was inversely ranked between the two cell lines in 2D cultures. Furthermore, comparing the 2D cultured MSC-hTERT (Fig 4G) with the 3D cultured MSC-hTERT (Fig 4I), revealed the inverted ranking of TFRC; ranked last in 2D MSC-hTERT whilst simultaneously the best reference gene in 3D culture.
When comparing the comprehensive ranking of all samples (Fig 4F) with only 2D (Fig 4J) or only 3D samples (Fig 4I), TFRC is ranked in the middle except for the 3D samples where it is considered the most stably expressed gene. Cancer stem cells (CSCs), e.g. from sarcoma, a cancer arising from mesenchymal cells, resemble MSCs in respect to self-renewal potential and differentiation capacity and can be used as a 3D model. PPIA, GAPDH, and YWHAZ have been shown to be amongst the best suitable reference genes for gene expression normalisation in CSC samples [22]. For our 3D cultured MSCs both PPIA and GAPDH ranked in the top 4 whereas YWHAZ ranked amongst the least stable genes. These differences indicate that careful consideration is required when selecting reference genes even when cell types are fairly similar or when culture dimension (e.g. 2D or 3D) is different. Reference genes for hMSC qPCR Although in the top 3 for all samples together, YWHAZ is an excellent candidate reference gene in hUC-MSC samples, in MSC-hTERT it is amongst the lesser stable genes and for PPIA and GAPDH the opposite is true. Furthermore, YWHAZ ranks in the top 3 for UC-MSCs ( Fig  4H) and overall (Fig 4F) whereas it is in the bottom three in both 3D cultured MSC-hTERT ( Fig  4I) and 2D-MSC-hTERT (Fig 4G) but back in the top 5 for all MSC-hTERT (Fig 4K) samples. These differences in rankings illustrate the necessity to use multiple candidate reference genes to determine their stability and their usefulness for accurate gene expression normalisation for each and every experiment. As indicated some samples showed no detectable expression levels for some genes (Fig 4A). For the 2D cultured UC-MSCs 1-8 samples (out of 27) no expression was detected in five of the 12 genes and only one sample (out of 9) of the 3D cultured MSC-hTERT showed no detectable expression for two candidate reference genes. The MSC-hTERT cells cultured in 2D (n = 6) showed detectable expression levels for all samples in all genes. Two of these genes (PUM1 and RPL13A) ranked in the overall top five (Fig 4F) and the 2D-UC-MSC (Fig 4H) ranking whereas they dropped to some extent in the other rankings. This further indicates that careful selection and assessment of candidate reference genes is needed depending on cell type.
The literature advises the use of at least three control genes for normalisation [9]. To illustrate the effect of using inappropriate reference genes we performed additional qPCR for MSC marker genes Endoglin (ENG or CD105), THY1 (CD90), and 5'-Nucleotidase Ecto (NT5E or CD73) (primer details in Table 1) and implemented normalisation using the three most stably (TBP, YWHAZ, and RPL13A) or least stably (PPIA, HPRT1, and B2M) expressed reference genes. The respective standard deviations increased when the least stably expressed reference genes were used (S1 Fig). The effect was less pronounced or even adverted for the MSC-hTERT samples since the three most stably expressed reference genes for these samples were TBP, PPIA, and TFRC. The top three ranked genes for the UC-MSCs were the same as for the overall comprehensive top three ranked genes.
The three sample groups (2D UC-MSC, 2D MSC-hTERT, and 3D MSC-hTERT) investigated showed much better correlation between the delta-delta Cq (ΔΔCq) values when using the three most stably expressed reference genes compared to the least stably expressed reference genes (Fig 5 and S2 Fig). The coefficient of variance (CoV) was also much lower when using the best reference genes (Table 4). Nevertheless, the CoVs for the three groups further decreased when only TBP and YWHAZ were used for normalisation and even more so when using only YWHAZ. In contrast, using TBP alone for normalisation resulted in an increase of the CoV despite this gene being the most stably expressed gene for most of the rankings. As noted earlier, the use of multiple reference genes is advised. Although gene expression normalisation with only YWHAZ resulted in the lowest CoV between samples in all groups, we would still recommend the use of multiple reference genes since the results further indicate that the best ranked gene might not be the best single reference gene for normalisation.
Potentially, all candidate reference genes are still suitable for gene expression normalisation. Less stable in this experiment could be the favorite gene for another experiment. All or a selection of the validated primer sets can be used for future experiments to perform accurate gene expression normalisation, although it is advised to always include TBP in these assessments. Gene expression detection range (Cq) was analyzed (A) in 2D UC-MSC (blue; n = 27), 2D cultured MSC-hTERT (green; n = 6) and 3D cultured MSC-hTERT (red; n = 9) samples. Samples with no Cq-values (number states amount) are shown as Cq > 40 (diamonds). Black, dotted line is maximum cycles in qPCR run. Boxes represent interquartile range (IQR = Q1 -Q3). Single, triangles represent outliers (Q3+1.5xIQR or Q1-1.5xIQR) if present. Whiskers represent maximum and minimum value after omitting outliers. Means are depicted by dots and medians by horizontal lines in the boxes. Gene expression stability according to Genorm (B), BestKeeper (C), Normfinder (D), and the Delta Ct method (E) were determined as well as the overall comprehensive ranking (F). Comprehensive rankings for all UC-MSC samples (G), only the 2D cultured MSC-hTERT (H), only 3D cultured hMSC-hTERT samples (I), all 2D cultured samples (J), and all hMSC-hTERT samples (K) were calculated using RefFinder. https://doi.org/10.1371/journal.pone.0209772.g004 Reference genes for hMSC qPCR

Conclusion
To accurately normalise gene expression data a panel of 12 reference genes has successfully been developed and assessed for their stability in hMSC cultured in 2D and 3D. Differences between 2D and 3D cultured samples and the two cell lines used can be observed in the ranking of the candidate reference genes indicating the influence of experimental properties and careful selection of reference genes for normalisation. The inconsistency between rankings further illustrates the necessity of using multiple reference genes for normalisation and moreover, no universally suitable gene for normalisation exists. This panel of validated candidate reference genes is available for gene expression analysis of MSCs although determining which genes to use needs to be established per experiment. Devoid of any known pseudogenes and amongst the most stable genes in all rankings, TBP is recommended for inclusion in any reference gene assessment in hMSC culture experiments.