affymetrix.com/analysis/index.affx). Quality control, normaliza- tion and statistical modeling were performed by array group correlation, mixed model normalization and mixed model analysis respectively. The normalized expression values for each gene were standardized by linearly scaling the values across all samples of the time course to a mean of 0 with an SD of 1. Analysis of differential gene expression was based on loglinear mixed model of perfect matches [44]. A false discovery rate of a = 0.05 with Bonferroni- correction for multiple testing was used to set the level of significance. The raw and normalized data are deposited in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih. gov/geo/; accession no. GSE40942) and followed MIAME requirements. 4.7 Microarray data analysis The statistical analysis of the microarray data was based on the normalized mean expression values per probe at 6 time points with 2 replications at each time point (12 observations per probe). In order to identify subgroups of probes with similar expression profiles over time, a principal component analysis (PCA) of the co- variance matrix was carried out on the mean expression value for each probe at each time point. For each probe, factor scores for principal components 1, 2 and 3 were obtained by regression analysis of the 12 array results (6 time points in duplicate) for that specific probe to those components. The first principal component corresponded with the general expression level during the whole experiment, whereas the second and third component corre- sponded with changes over time. Since our interest was to identify genes associated with the changes that occur during differentiation from stromal cells towards chondrocytes, we focused our analysis on the second and third component. By construction, these factor scores had a mean of 0 with an SD of 1. Generally, the distribution over the factor scores showed a normal distribution with outliers. We used a cut-off of $3.29 to select outlying probes. This cut-off would select 0.1% of the probes, if the factors scores would follow a pure normal distribution that could be expected if the data were pure noise. The presence of replications allowed us to assess the statistical significance of the factor scores and to remove probes that were not significant at the a = 5% level. 4.8 Pathway and network analysis Using sets of probes emerging from PCA, a search for relevant KEGG pathways was performed using the DAVIDH Knowledge- base, a publicly available bioinformatics tool for functional annotation (http://david.abcc.ncifcrf.gov). To investigate known and predicted associations between the protein products of the genes that were more than 3.29-fold changed between week 0 and week 5, we used STRING9.0 (string-db.org). To visualize known and predicted associations between protein products of the genes that were differentially expressed between undifferentiated hfMSCs and pre-pubertal growth plate samples we used STRING9.0. Growth plates of 3 donor with Tanner stage B1 were included (one 12 year-old female with chondroid osteosarcoma, one 8 year-old female with osteosarcoma, one 10 year-old male with cerebral palsy). 4.9 Quantitative real-time polymerase chain reaction (qPCR) RNA was transcribed into cDNA using the First Strand cDNA Synthesis kit for qPCR (Roche Diagnostics) according to the hfMSCs Differentiate toward Growth Plate Cartilage manufacturer’s protocol. Specific primer sets (available at www. utwente.nl/tnw/DBE) were designed to amplify aggrecan (ACAN), pannexin 3 (PANX3), epiphycan (EPYC), collagen type II (COL2A1), and type X (COL10A1), SRY-box 9 (SOX9), WNT11, lymphoid enhancer-binding factor 1 (LEF1), Gremlin 1 (GREM1). b2-Microglobulin (B2M), and glyceraldehyde 3-phosphate dehy- drogenase (GAPDH) were used as housekeeping genes. The expression of these housekeeping genes was stable during differentiation of hfMSCs as shown by the microarray analysis. In order to test donor inter-variation, differentiated MSCs isolated from other fetal donors were used for qPCR analysis as well. All PCR reactions were performed in triplicate with 5 ng cDNA and according to the manufacturer’s protocol of the iQTM SYBRH Green Kit (Biorad) in a final volume of 20 ml. The cDNA was amplified using the following thermal cycling conditions: one cycle at 50uC for 2 min and 95uC for 10 min, followed by 40 cycles of 15 s at 95uC, 15 s at 60uC and 15 s at 72uC and finally a melting curve was generated using the mIQ Single-Color—Real-Time PCR System (BioRad Laboratories, Hercules, California, USA). Fluorescence spectra were recorded and the threshold cycle number (Ct) was calculated using the accompanying mIQ- software. For each time point mean Ct was calculated and from this value the fold difference in expression between undifferenti- ated hfMSCs and differentiating cells using the 22DDCt method was calculated using B2M as a reference as described by Schmittgen and Livak, 2008 [45]. For visualization, this value was log-transformed. Supporting Information Figure S1 There is a good correlation between qPCR and microarray expression data. For (A) ACAN, (B) COL2A1, (C) COL10A1, (D) SOX9, (E) PANX3, (F) EPYC, (G) WNT11, (H) LEF1 and (I) GREM1 during 5 weeks of chondrogenic differen- tiation of hfMSCs. qPCR data are expressed as delta delta CT values corrected for the housekeeping gene B2M. The primary y- axis (left) indicates the qPCR results as normalized mean fold expression on a log-scale. The secondary y-axis (right) indicates the microarray analysis results as least square means (lsm). (TIF) Table S1 List of genes selected with principal compo- nent analysis. (PDF) Table S2 List of top hits of up regulated/down regulat- ed genes at each time point with their fold change in expression compared to undifferentiated hfMSCs. (DOC) Table S3 List of top hits of differentially expressed genes of pre-pubertal GP with their fold change in expression compared to undifferentiated hfMSCs. (DOC) Author Contributions Conceived and designed the experiments: SvG JE JCHL CS JCvH JJG CK SS NG JMW GR JNP MK. Performed the experiments: SvG JE JCHL ED CS JCvH JJG CK SS NG JMW GR JNP MK. Analyzed the data: SvG JE JCHL ED CS JCvH JJG CK SS NG JMW GR JNP MK. Wrote the paper: SvG JE JCHL ED CS JCvH JJG CK SS NG JMW GR JNP MK. PLOS ONE | www.plosone.org 10 November 2012 | Volume 7 | Issue 11 | e44561