Identification of Thalidomide-Specific Transcriptomics and Proteomics Signatures during Differentiation of Human Embryonic Stem Cells

Embryonic development can be partially recapitulated in vitro by differentiating human embryonic stem cells (hESCs). Thalidomide is a developmental toxicant in vivo and acts in a species-dependent manner. Besides its therapeutic value, thalidomide also serves as a prototypical model to study teratogenecity. Although many in vivo and in vitro platforms have demonstrated its toxicity, only a few test systems accurately reflect human physiology. We used global gene expression and proteomics profiling (two dimensional electrophoresis (2DE) coupled with Tandem Mass spectrometry) to demonstrate hESC differentiation and thalidomide embryotoxicity/teratogenecity with clinically relevant dose(s). Proteome analysis showed loss of POU5F1 regulatory proteins PKM2 and RBM14 and an over expression of proteins involved in neuronal development (such as PAK2, PAFAH1B2 and PAFAH1B3) after 14 days of differentiation. The genomic and proteomic expression pattern demonstrated differential expression of limb, heart and embryonic development related transcription factors and biological processes. Moreover, this study uncovered novel possible mechanisms, such as the inhibition of RANBP1, that participate in the nucleocytoplasmic trafficking of proteins and inhibition of glutathione transferases (GSTA1, GSTA2), that protect the cell from secondary oxidative stress. As a proof of principle, we demonstrated that a combination of transcriptomics and proteomics, along with consistent differentiation of hESCs, enabled the detection of canonical and novel teratogenic intracellular mechanisms of thalidomide.


Introduction
Traditional approaches to toxicological testing typically involves exposure of chemicals to large numbers of animals during the crucial period of organ development and further investigations of foetuses for visceral and skeletal developments, these approaches are expensive and time consuming [1][2][3][4]. In order to provide costefficient and high throughput methods, a multitude of in vitro test systems have been proposed to assess the developmental toxicity of candidate drugs and environmental toxicants in the past 20 years. These platforms include primary in vitro cell cultures and ex vivo models using embryo cultures [5] Embryonic stem cells (ESCs) have the unique potential to differentiate into all somatic cell types. In this context, a mouse ESC test originally covering the three end points for predicting teratogenicity (ESC cytotoxicity, fibroblasts cytotoxicity, and the inhibition of ESC differentiation into cardiomyocytes) has been initiated [6,7] . Although these methods are able to predict toxicity of the drugs, hESC were introduced for toxicity testing in order to better reflect the human physiology and to avoid interspecies differences [8].
The embryotoxic drug thalidomide (Contergan) was launched in 1957 in Germany and was subsequently withdrawn in November 1961 after its teratogenic effects in humans were recognised. The clinical evidences showed that thalidomide causes various phenotypic malformations such as limb, ear, ocular, kidney, heart and gastrointestinal deformities (Reviewed in [9,10]). To determine the teratogenic potential of thalidomide, various in vivo assessments were applied to different animal species based on traditional clinical and histopathological measurements. Thalidomide induced distinct developmental adverse effects in different animals such as dogs, rats, mice and rabbits [11,12]. Moreover, congenital malformations induced by thalidomide were prominently found in rabbits where as in rats moderate effects were found and in mice no significant foetal changes were observed [11,13,14]. Transcriptional profiling of cynomolgus monkeys at 26-28 days of gestation, observed limb defects and downregulation of vasculature development-related transcripts [15].
To unravel the molecular mechanisms of thalidomide in embryonic development various in vitro models has been utilised, predominantly primary cells from non human origin and many hypotheses or mechanisms have been proposed for thalidomide teratogenecity [16][17][18][19]. It was demonstrated that generation of reactive oxidative species (ROS), oxidative DNA damage and perturbation of intra cellular signalling such as FGF, WNT and AKT are the causes of thalidomide induced limb deformities [9,20]. The generation of ROS leads to oxidation or alteration of glutathione content. This is crucial for detoxification of cells, especially oxidative stress conditions and regular embryonic development [21]. Although these studies described the required toxic dose and perturbations of organ development, these results cannot necessarily be extrapolated to other species or humans due to the known species-specificity of thalidomide [22]. Therefore, a consistent and predictive developmental toxicity model based on hESCs requires an in-depth insight into the molecular mechanisms that explain the adverse developmental potential of a drug in the concentration range applied under in vivo conditions. Omics approaches using ESCs as a model were proposed as an innovative way of drug safety testing (reviewed in [23]).
Multilineage differentiation of human embryonic stem cells (hESCs) can partially reproduce early human embryonic development [24]. Therefore, hESCs are a suitable tool to assess toxicant profiles to understand and predict the damage caused by potential therapeutic agents. To test the applicability of the hESCs as a proof of principle human in vitro embryotoxicity model, we applied a well-defined differentiation protocol. We studied thalidomide, a highly investigated in vivo developmental toxicant, at different concentrations with a combination of sensitive transcriptomics and proteomics approaches. Here, we demonstrated that low concentrations of thalidomide, normally observed in the plasma under in vivo conditions, induced the perturbation of genes and proteins participating in limb, heart development and oxidative stress-mediated nuclear cytoplasmic protein transport.

Ethics Statement
Experiments were approved by the governmental animal care and use office (Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen, Recklinghausen, Germany, Protocol No. 8.84-02.05. 30.11.030) and were in accordance with the German Animal Welfare Act.
hESCs culture conditions and differentiation H9 hESCs (WiCell, Madison, WI, USA) were cultured in DMEM-F12, 20% KO serum replacement, 1% non-essential amino acids, penicillin (100 units/ml), streptomycin (100 mg/ml) and 0.1 mM b-mercaptoethanol supplemented with 4 ng/ml basic fibroblast growth factor (bFGF). The cells were passaged with mechanical dissociation on irradiated mouse embryonic fibroblasts (MEF) which were prepared as described previously [25]. Prior to differentiation, the cells were maintained for five days in 60-mm tissue culture plates (Nunc, Langenselbold, Germany) coated with a hESC-qualified matrix (BD Biosciences, California, USA) in mTESR medium (Stem Cell Technologies). For the time kinetic multilineage differentiation, embryoid bodies (EBs) were prepared as described previously [24] with minor changes (60 to 70 clumps were added instead of 50 to 60), and the EBs were maintained for 21 days on a horizontal shaker ( Figure 1A). For the dose-response analysis, thalidomide (0.01 mM, 0.1 mM, 1 mM, 10 mM, or 70 mM) (Sigma, Steinheim, Germany) was added to the medium from day 0 (undifferentiated ESCs) to day 14. Equal volumes of the highest solvent concentration were added as controls. For 2DE analysis, a 10 mM thalidomide concentration was used. On alternating days, the medium was replaced with fresh medium containing the drug. All of the experiments for microarray and 2DE analysis were performed as three independent (n = 3) biological replicates.

Microarray labelling and hybridisation
To isolate total RNA for time kinetic analysis, samples were collected from undifferentiated hESCs and differentiated EBs for every 3 days from day 3 to day 21 ( Figure 1A). For the doseresponse analysis, samples were procured from hESCs, 14-day-old EBs treated with thalidomide (0.01 mM, 0.1 mM, 1 mM, 10 mM, or 70 mM) and the controls. The samples were homogenised in Trizol (Invitrogen, Darmstadt, Germany), and RNA was extracted using Trizol and CHCl 3 (Sigma, Steinheim, Germany). The total RNA was purified using an RNeasy mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. A Nanodrop (ND-1000, Thermo-Fisher, Langenselbold, Germany) was used to measure the RNA quantity. The integrity of the RNA was confirmed using denaturing agarose gel electrophoresis. For the transcriptional profiling of human cells, the Affymetrix Human Genome U133 plus 2.0 arrays were used. For RNA amplification, 100 ng total RNA was used with a Genechip 39 IVT Express Kit (Affymetrix, High Wycombe, United Kingdom) according to the manufacturer's instructions. For the aRNA purification, magnetic beads were used, and 15 mg of aRNA were fragmented. Fragmented aRNA (12.5 mg) was hybridised to microarrays in the Genechip Hybridisation Oven-645 (Affymetrix, High Wycombe, United Kingdom). The Affymetrix HWS kit was used for staining and washing with a Genechip Fluidics Station-450 according to the manufacturer's instructions. The stained arrays were scanned with an Affymetrix Gene-Chip Scanner-3000-7G, and the quality control analysis was performed with Affymetrix GCOS software.

Data analysis and statistical procedures
Background correction, summarisation and normalisation were performed with a Robust Multi-array Analysis. The raw dataset was normalised with the Quantile normalisation method (executable with a R-affy package) performed at the probe feature level. Probe sets having a signal of #6 (log2 scale), in any one of the conditions out of 21, were chosen for statistical analysis. The differentially expressed genes (DEG) were described by a linear model implementing R-LIMMA packages (Linear Models for Microarray Data). A one-way ANOVA calculation was performed, considering either the treatment or the days of differentiation as a factor. A pairwise comparison was performed, and the p-values of the scores of the moderated t-test calculation were used further for statistical correction. Benjamini-Hochberg method was used to adjust the raw p-values to control the false discovery rate. Differentially expressed transcripts were filtered with a FDRcontrolled p-value of #0.05 (95% confidence interval) along with fold change values at a threshold value of $62. To find the doseresponse or time kinetic effects across the dose or time, DEG were derived from F statistics. The distances of the samples were calculated using a Spearman Rank correlation. A principal component (PC) analysis was performed using the stats package in R with a prcomp function. The "x" attribute of the prcomp object was used to generate three-dimensional scatter plots. B, An unsupervised hierarchical clustering analysis demonstrates that a temporal gene expression pattern changes during differentiation and the major changes were observed between days 12 to 14. Three independent biological experiments were performed for statistical analysis (p#0.05). C, The dendrogams pattern of array cluster showed that a two major negative correlation for day 0 to 9 and 12 to 21 differentiation. D, The selected germ layer lineage specific BP (Enrichment p#0.01) for up-regulated transcripts demonstrated that the numbers of genes are optimal between days 12 and 15. The figure represents (i) endoderm-specific, ii) ectoderm-specific and iii) mesoderm-specific BP. The x-axis represents the days of differentiation, and the y-axis represents number of transcripts. doi:10.1371/journal.pone.0044228.g001

Gene expression and functional annotation analysis
To find the common DEG for thalidomide treatment an intersect analysis were performed. Briefly, the DEG from the temporal gene expression was used to calculate probable intersects among them. The counts are represented in the plot as number of genes. For temporal and dose response hierarchical clustering, the intensity values of the DEG from F statistics were used for uncentred correlation and average linkage clustering and visualised using Java Treeview. Transcripts were filtered with gene vector 0.2 and resulting 18567 transcripts for temporal and for dose response experiments 1076 genes from 70 mM were used for cluster analysis. The red colour represents up-regulated transcripts, and the blue represents down-regulated transcripts. To investigate the DEG, functional annotation and gene ontology (GO) clustering, the Database for Annotation, Visualisation and Integrated Discovery (DAVID) was used [26]. These analysis explore the biological processes (BP) related to the transcripts, which were differentially regulated in time kinetic and doseresponse analysis. To increase stringency, GO categories with less than five transcripts were eliminated and filtered further based on the EASE score enrichment p-value (P#0.01) and clustering stringency at the medium level as per DAVID.

Quantitative real time PCR
An independent dose-response experiment was performed with thalidomide (1 mM, 10 mM and 70 mM), and the total RNA was extracted as described above. cDNA synthesis was performed with the Super Script Vilo cDNA synthesis kit (Invitrogen, Darmstadt, Germany) using 2 mg of total RNA as starting material. The cDNA was diluted 20 times with nuclease-free water, and 100 ng was used as a template for PCR. The primer design was performed using Primer3, or the sequences were obtained online (http:// www.origene.com/). The assay was accomplished using the PlatinumH SYBRH Green qPCR SuperMix (Invitrogen, Darmstadt, Germany) with 0.5 mM of the primer using the Applied Biosystems 7500 FAST cycler. The primer sequences are listed in Table S1. All of the gene expression values were normalised to the reference b-actin. The mRNA expression levels were represented as a relative quantification.

Immunostaining
From an independent dose-response experiment on day 12 (1 mM, 10 mM, and 70 mM), EBs were dissociated with 0.05% trypsin for 10 minutes, and single cells were seeded on a fibronectin coated coverslip. On day 14, the cells were fixed with ice-cold methanol for 10 minutes at 220uC. The cells were then rehydrated three times (10 min each) using phosphate buffered saline with Ca ++ and Mg ++ (PBS). The samples were blocked with 5% bovine serum albumin (BSA) for 60 minutes and incubated with the primary antibody (b-Catenin (BD Biosciences, California, USA), GSTA1, GSTA2 (Abgent Europe, United kingdom), RANBP1 (Calbiochem, Darmstadt, Germany) overnight. The samples were then washed three times with PBS for 15 minutes and incubated with the respective secondary antibody conjugated to Alexa Fluor 488/594 (Millipore, Schwalbach, Germany) for 60 minutes. The samples were washed three times with PBS for 15 minutes and then mounted with Prolong Gold anti-fade with DAPI (Invitrogen, Darmstadt, Germany). The samples were observed with a Zeiss Axiovert 200 fluorescence microscope (Carl Zeiss Microscopy, Oberkochen, Germany).

Two-dimensional gel electrophoresis (2DE) and Proteomics analysis
The frozen cell pellets were dissolved in a sample buffer for isoelectric focusing (6 M urea, 3 M thiourea, 2% ampholytes, and 70 mM DTT). The proteins were separated by 2-DE according to a protocol previously used by Klose and Kobalz [27] with modifications by Zimny-Arndt [28]. The proteins (100 mg total protein amount) were anodically loaded and focused in the first dimension (isoelectric focusing) at 8,800 Vh in a gradient between pI 3 and pI 11. The isoelectric focusing gel was applied on top of a 15% SDS-PAGE acrylamide gel (gel size 20 cm630 cm; thickness of 1.0) to separate the proteins by size in the second dimension. The SDS-PAGE gels were fixed in 50% ethanol and 10% acetic acid for approximately 16 hours. The gels were silver stained with a mass spectrometry compatible silver staining kit (Proteome Factory PS-2001). The stained gel images were aligned using DECODON Delta2D software (http://www.decodon.com) within the biological replicates. To derive the differentiation markers, 14 day old EB gels were aligned with hESCs. The thalidomidespecific protein treatment gels were aligned with the untreated samples using warp transformations. For spot detection and quantification, the fusion gel was prepared from the respective replicates using the spot preservation and image fusion algorithm from Delta2D. The spot pattern was transferred from the fused image to all of the original images. A typical 2D image is shown in Figure S2. To derive the significantly regulated protein spots, the ttest was applied with a statistical value of p#0.01.

Tryptic digestion
The protocol for tryptic digestion was based on a protocol of Bertinetti et al with some modifications. 2DE-spots were picked and incubated with a shrinking solution consisting of 50 mM ammonium bicarbonate (Sigma Aldrich, Germany) in a 60% acetonitrile/water solution (Lichrosolv, hyper grade, Merck, Germany) for 5 minutes. The shrinking buffer was removed, and the gel was incubated with a swelling buffer containing 100 mM ammonium bicarbonate. Each gel band was placed into a 10-ml pipette tip that was then placed in a shortened 1000-ml tip and inserted into a borosilicate vial (CS-Chromatographie Service, Germany). This apparatus was placed in a 1.5 ml tube and centrifuged for 5 min at 14,000 rpm (Centrifuge 5415C, Eppendorf, Germany). The obtained gel slurry was dried by evaporation (RC10, Thermo Fisher Scientific, USA), and the proteins were digested overnight at 37uC after the addition of 20 ml trypsin (final concentration 10 ng/ml, dissolved in trypsin resuspension buffer, Promega, USA) to the gel. The reaction was stopped by the addition of a 5% formic acid and water solution (final concentration 0.1%), and the supernatant was taken from the gel slurry. The peptides were extracted twice by the addition of 10 ml of a solvent containing 65% acetonitrile with 35% formic acid (5% v/v in water) and 10 ml of 100% acetonitrile, respectively. The resulting peptide fractions were evaporated until they were completely dry. The fractions were then resolved in 10 ml of solvent (5% acetonitrile and 0.2% formic acid in water) prior to LC-MS analysis.

Protein identification by LC-MS-Analysis
Identification was performed on an Agilent 1100 LC/MSD-trap XCT series system. The electrospray ionization system was the Chip Cube system using a ProtID-Chip-43 (Agilent Technologies). Sample loading from the microtiter plate onto the enrichment column was performed at a flow rate set to 4 ml/min with two mobile phases at a ratio of 98:2 (mobile phase A: 0.2% FA in H2O; mobile phase B: 100% ACN). LC gradient was delivered with a flow rate of 500 nl/min. Tryptic peptides were eluted from the reversed-phase column into the mass spectrometer using a linear gradient elution of 2-30% B in 20 min. For MS experiments, the following mode and tuning parameters were used: Scan range: 300-2,000 m/z, polarity: positive, capillary voltage: 21,800 V, flow and temperature of the drying gas were 4 l/min and 325uC. The MS/MS experiments were carried out in auto MS/MS mode using a 4 Da window for precursor ion selection. After 3 MS/MS spectra, the precursor ions were actively excluded from fragmentation for at least 1 min. The generic files for database searching were generated by Data Analysis software version 4.0. For precursor ion selection, a threshold of 5 S/N was applied and the absolute number of compounds was restricted to 1,000 per MS/MS experiment (according to (F-23) with modifications. Protein identification was performed by Mascot software (Version 2.3) using the UniProtKB database (http://www.uniprot. org/)

Western blotting
For western blot analysis, the samples procured from the undifferentiated hESC, thalidomide treated (1 mM, 10 mM, and 70 mM) and untreated EBs at day 14 was used (An independent experiment was performed for validation). The samples were washed with PBS and lysed with ProteoJET (Fermentas, Germany) mammalian cell lysis reagent. The clarified samples were quantitated with the Bradford reagent (Sigma, Steinheim, Germany). The protein (30 mg) was separated using either 10% or 12% SDS-PAGE gels and then transferred to a PVDF membrane. The membranes were blocked with 5% BSA at room temperature and incubated with RANBP1 (Calbiochem, Darmstadt, Germany), KPNA2, DDAH2 (Abcam, Cambridge, UK), GSTA1, GSTA2, GATA6 and HAND1 (Abgent Europe, United kingdom), b-Catenin (BD Biosciences, California, USA) primary antibodies for overnight at 4uC on an orbital shaker. After washing with PBST (PBS with 0.1% Tween-20), the membranes were incubated with the respective secondary antibody (labelled with HRP) for 1 hour at room temperature. The antibodies were detected using the ECL Western Detection System (Fisher Scientific GmbH, Schwerte, Germany).

Time kinetic analysis of hESC differentiation
A global transcriptomic analysis for time kinetic differentiation was performed to determine the optimal differentiation time point that showed the highest expression levels of developmental genes ( Figure 1A). A total of 18567 DEG, derived from the 'F statistics' (Table S2A), were used for hierarchical cluster analysis. The gene cluster demonstrated that the major up-regulated transcripts occurring up to day 9 were repressed after day 12. The over expressed genes between days 12 and 15 were consistently induced up to day 21 ( Figure 1B). Furthermore, the gene array dendrograms results showed that days 0 to 9 and days 12 to 15 were negatively correlated (20.99) ( Figure 1C). To find the common genes between the time points, intersect gene analysis was performed (as mentioned in methods). A consistent expression of up-regulated genes were observed from days 12 to 21 and compared to days 6 to 12 (selected intersects are shown in Figure  S1, and the complete list given in Table S2B). To derive the biological rational for the up-regulated transcripts (Table S2C to S2I), DAVID GO analysis was performed. The BP for various organ developments (represented by the endoderm, ectoderm and mesoderm lineage) demonstrated that maximum numbers of genes were regulated between days 12 to 15 of differentiation ( Figure 1D). These results confirm that differentiation for 12 to 15 days can demonstrate the embryonic development with optimal gene expression of development related genes.

hESC differentiation for 14 days can demonstrates embryonic development
In order to uncover multiple embryonic development perturbations in presence of potential toxicants, it is critical to determine the optimal time point of hESC differentiation. From time kinetic experiment it was observed that between 12-15 days embryonic development related genes and biological processes were highly up-regulated, suggesting 14 days hESC differentiation as an optimal period for further study. The differentiation scheme is briefly described in Figure 2A. Comparison between ESCs after a 14-day differentiation in microarray analysis showed 5035 differentially expressed transcripts (Table S3A). The hierarchical cluster analysis of DEG showed two distinct clusters for 14 days of hESC differentiation representing up and downregulated transcripts ( Figure 2B). To further unravel the biological correlation of up-regulated genes, GO analysis were performed. BPs related to the neuron, skeletal system, vasculature and embryonic organ development were enriched among the upregulated genes after 14 days of hESC differentiation (Table  S4A). The neuron development BP gene expression pattern is represented in volcano plot ( Figure 2C). A reduction of pluripotency genes is crucial for embryonic development and the array results corroborates a down-regulation of POU5F1, NANOG, and LIN 28 ( Figure 2D). For proteomic analysis, Delta2D analysis (p#0.01) of 2DE gels (as mentioned in methods) for hESCs and 14-day old EBs identified 281 regulated spots. A total of 61 protein spots were selected for further mass spectrometry. The spots encompassed 33 up-regulated and 28 down-regulated proteins, at a fold change $61.5 (Table S5). The comparative analysis of genomics and proteomics for hESC differentiation unveiled 10 commonly regulated genes (Table 1). These included PAK2 and ENO2, which are highly expressed during neuronal development [29,30]. Among the 33 upregulated proteins, 11 were shown to be involved in neuronal differentiation and brain development, particularly the Reelin pathway regulators PAFAH1B2 and PAFAH1B3 ( Figure 2E). Within down-regulated proteins, POU5F1 regulators (such as TPIM28, PKM2 and RBM14) have been identified in 14-day-old EBs, suggesting a progressive differentiation of the hESCs ( Figure 2F).

Dose-response gene expression profiling in response to thalidomide
Five concentrations were selected for the dose-response study and the differentiation scheme in presence of treatment was briefly demonstrated in Figure 3A. The morphology of EBs for all concentrations (0.01 mM, 0.1 mM, 1 mM, 10 mM and 70 mM) showed no reductions or alterations in the EB size due to cytotoxicity ( Figure 3B). To determine the cytotoxic concentration, hESCs were treated with higher concentrations of thalidomide (up to 1 mM) and no cytotoxicity in Resazurin reduction assay was observed (data not shown).
To visualise the global gene expression pattern, principal component analysis (PCA) was performed. PCA projects the high content data with clustering of similar expression pattern of samples. Figure 3C shows merely any changes at the lowest concentration and at higher concentration significant difference from control was observed. The samples subjected to PCA demonstrate larger variances (PC#29.9%) at 1 mM and above concentrations. At a concentration of 0.01 mM, 5 DEG were observed and therefore was considered to be a low observed adverse effect level (LOAEL) concentration. The numbers of dysregulated transcripts increased progressively with thalidomide treatment and the highest thalidomide concentration (70 mM) yielded 803 down-and 223 up-regulated genes ( Figure 3D) (Table  S3A to S3F).

Transcriptomic signatures and GO assessments for thalidomide
To determine the effect of thalidomide across different concentrations, F-statistics (From thalidomide dose response study) were applied (FDR corrected p-value#0.05, 1-way ANOVA,) which resulted in 25433 DEG (Table S3G). The hierarchical clustering analysis of DEG showed dose dependent inhibition of gene signatures from 1 mM to 70 mM concentration ( Figure 3E). To identify the functional relevance of these genes, functional  annotation clustering analysis was performed separately for upand down-regulated transcripts for 1, 10 and 70 mM (Table S6A to S6F). The representative GO for down-regulated genes are involved in heart development, limb development, vasculature development, skeletal system development, and DNA-dependent regulation of transcription ( Table 2). The GO quantification of developmental biological processes showed that thalidomide induced significant developmental toxicity in a dose dependent manner ( Figure 3F). The complete lists of GO are provided in Tables S6A to S6F for 1 to 70 mM.

Transcription factor analysis
GO analysis of the down-regulated transcripts (at 70 mM) uncovered 96 transcripts for DNA-dependent, regulation of transcription ( Table 2). An independent functional annotation analysis shows that these factors are involved in embryonic organ development, skeletal system development, sensory organ development and pattern specification processes (Table S7). The list of genes contributing to each GO consists of Homeobox, Tbox, helix loop helix and members of the zinc finger transcriptional factor families.

Perturbation of heart, limb development and WNT related genes
The differential expression in heart and limb development BP coincides with clinical evidence of thalidomide toxicity [9,10]. The representative genes from the respective BP uncovered transcriptional factors (such as GATA6, HAND1 for heart development and, HOXA10, HOXA9, TBX3, and HOXB8 for limb development) which were up-regulated in untreated 14-day differentiated cells. In contrast, negative regulation of these genes was observed for thalidomide treatment ( Figure 4A). Additionally, we found that groups of genes involved in WNT signalling pathway, such as FZD10 (a WNT activator), AXIN (a b-catenin assembly stabiliser) and TCF/LEF (a WNT signalling target) were down-regulated by 2-fold ( Figure 4A). The Delta 2D analysis of 2DE gels for thalidomide (10 mM) revealed 6 down-regulated protein spots at p#0.01 and fold change $61.3 (Table S8). These spots were selected for mass spectrometry analysis ( Figure S3). The results showed a 1.9-fold down-regulation of dimethylarginine dimethylaminohydrolase 2 (DDAH2), which regulates the nitric oxide synthase via destruction of asymmetric dimethylarginines and is involved in limb development ( Figure 4B (II)). Regulation of selected markers for heart and limb development was confirmed with immunoblotting and RT-qPCR analysis ( Figure 4B, 4C). To find the interference of thalidomide in early cardiac development an independent temporal and dose response experiment was performed. The RT-qPCR analysis demonstrates cardiac specific transcription factors such as NKX2.5, GATA4, GATA6 and HAND1 were down-regulated ( Figure S4). The cellular localization of b-catenin is showed with immunostaining ( Figure 4D). The dysregulated genes related to anteroposterior (A-P) limb development in response to thalidomide are represented in Figure 4E.

Thalidomide-induced dysregulation of glutathione transferases (GST) and nucleocytoplasmic trafficking
Glutathion depletion or oxidation and ROS generation is another proven mechanism for thalidomide teratogenecity [20,31]. From microarray experiment we found GSTA1 was over expressed by 3 fold in 14 days old EBs and down-regulated by 2fold change for thalidomide. The GST family are involved in detoxifying xenobiotics and inactivating secondary metabolites (such as hydroperoxides and epoxides) during oxidative stress [32].
To further reconfirm glutathione family genes expression such as GSTA1, GSTA2 and GSTA3 real time PCR and immunoblotting analysis were performed which confirm dose dependent repression for thalidomide ( Figure 5A, 5B). To find the cellular localisation of GSTA1 and GSTA2, immunoflourescence analysis was performed for control which shows an enrichment of these respective proteins in nuclear membrane. In response to thalidomide positive signals were found in the nuclear membrane; however the enrichment of proteins collapsed ( Figure 5C).
The perturbation of generation and clearance of ROS in cellular machinery leads to the pronounced effect of oxidative stress. But the well known fact of oxidative stress induced nucleocytoplasmic traffic dysregulation was not studied previously in response to thalidomide. Microarray analysis showed repression of trafficking genes KPNA2, RANBP1, NUP98 and NUP160 by more than 21.5 fold (Table S3G). A mass spectrometry analysis showed the down-regulation of RANBP1 ( Figure 5D). A RAN binding protein present in the cytoplasm which is required for carriers nucleoporins and adaptor importin-a (KPNA2) [33]. The RT-qPCR and immunoblotting analysis confirmed the concentration dependent down-regulation of trafficking genes ( Figure 5D, 5E, 5F). Immunoflourescence analysis demonstrated that thalidomide induced alteration of cytoplasmic RANBP1 gradient ( Figure 5G).

Discussion
The present time kinetic experiment demonstrates that the optimum differentiation time point, which is relevant for developmental toxicity approaches, was reached after day 12. Additionally, the expression pattern was consistent up to day 21. The up-regulated transcripts of 14-day differentiated cells belong predominantly to various embryonic developmental biological processes, such as skeletal, neuronal, and respiratory system development. The significant number of transcripts expressed during differentiation emphasises a critical window to monitor molecular mechanisms and developmental pathways induced by developmental toxicants. Our comparative 2DE proteomics analysis resulted in 33 up-regulated proteins. Eleven of these represented proteins are involved in neuronal development ( Figure 2E). Among the 11 proteins, we found 2 proteins that belong to the Reelin pathway (PAFAH1B2 and PAFAH1B3) and that regulate brain development. PAK4 was also identified, which targets Rho GTPases, such as CDC42, that are required for neuronal development and axon outgrowth [29,34]. Up-regulation of RELN and CDC42 at the protein level also confirmed with mRNA levels uncovered by the microarray data ( Figure 2C).
To reveal the dose dependent gene expression pattern during hESC differentiation 5 concentrations which represent clinical plasma concentrations for therapeutic purposes were selected. As previously described, the mean steady state concentration for lower and higher doses of thalidomide ranges between 1 to 10 mM [35,36]. Highly expressed developmental genes at 14 days were down-regulated after treatment with increasing concentrations of thalidomide. Notably, BP and GO quantification results showed that thalidomide affected embryonic, heart and limb and development (Table 2, Figure 3F). Dose-response genomic data has an advantage because it demonstrates the specific dose(s) affecting the corresponding BP. The present study shows the alterations in BP associated with mammary glands at a dose of 0.88 mM (Table S9 BMD data), which was observed previously in the beagle dog [12]. Perturbation of heart and limb development was consistently shown for thalidomide both in vivo and in vitro [17,[37][38][39]. Therefore, these findings obtained under in vitro Table 2. Selected down-regulated significant (p#0.01) GO categories after thalidomide treatment. developmental conditions using hESCs suitably recapitulate the disturbances in heart and limb development observed in in vivo animal studies [11,15,39].
Canonical embryonic limb bud development in vertebrates proceeds along the proximal to distal (P-D), A-P and dorsal to ventral (D-V) axes [40]. P-D limb development is regulated by apical ectodermal ridge (AER) and FGF signalling [40]. The A-P Figure 4. Thalidomide perturbed heart, limb development and WNT signalling associated genes. A, The heat map shows the microarray expression pattern for the selected genes for heart, limb development (from Table S6D, S6F) and WNT signalling (from Table S9). Since multiple genes are involved in various developmental processes few redundant genes are present. The up-regulated genes in control were repressed in a dose dependent manner. B, The representative heart (I), limb development (II) associated genes were analysed using RT-qPCR analysis with an independent experiment. Thalidomide responsive genes were shown in I and II (*p-value#0.01, thalidomide-treated vs untreated 14-days old EBs) and 14 days old EBs gene expression were shown in (III) *p-value#0.01, 14 days old EBs vs ESCs . The error bars represents the SEM from 3 technical replicates. C, The mass spectrometry results showed down-regulation of DDAH2 gene (n = 3), (p#0.01). Figure   The error bars represent the SEM from 3 technical replicates from an independent experiment. B, Immunoblotting analysis of GSTA1/2 confirms thalidomide induced perturbation at protein concentration. C, The cellular localisation of GSTA1/2 showed pronounced concentration of GSTA on limb development depends on the expression of SHH and zone of polarising activity (ZPA). Using 14-day EBs, we observed high expression levels of the transcription factors HOXB8 and HAND2, which are essential for ZPA and SHH expression (a regulator of TBX3 and PTCH1) [41,42]. GREM1, a negative regulator of BMP, was downregulated by 10-fold. The similar expression pattern of 14-day EBs and A-P limb development genes in vivo suggest that this model is suitable to predict specific effects in A-P limb development. These genes were significantly dysregulated by thalidomide in a dose-dependent manner ( Figure 4E). Cereblon (CRBN) forms an E3 ubiquitin ligase complex with CUL4A and DNA binding protein1 and was identified as one of the molecular targets for thalidomide. It is also expressed during limb outgrowth in zebrafish and chicks, which depends on FGF signalling [43]. Because our model does not express genes involved in P-D limb development, such as FGF8 or CRBN, this mechanism of developmental toxicity cannot be confirmed. External addition of growth factors may enable our model to encompass P-D patterning.
Interestingly, expression of the homeobox transcription factors HOXA, HOXD and HOXB8 is also affected by thalidomide treatment in differentiated hESCs (Table S7, Figure 4C). It is well known that HOXA, HOXD and HOXB8 promote limb morphogenesis (reviewed in [44]) and early limb-bud formation [44,45]. WNT signalling controls the dynamic limb development processes such as initiation, patterning and differentiation [46]. The attenuation of markers involved in WNT signalling was observed in the present study ( Figure 4A, 4C, 4D) which coincides with the previous in vitro studies [17,18]. DDAH1/2 catalyses the hydrolysis of asymmetric dimethylarginine (an endogenous inhibitor of the nitric oxide synthase [47]) and is expressed in limb-bud patterning. We observed perturbations of DDAH2 at both the gene and protein level ( Figure 4B, 4C). Perturbations of A-P limb bud development genes in thalidomide-treated hESCs are shown in Figure 4E.
Transcription factor HAND1 expresses highly in developing left ventricle and shows profound expression in the right ventricle. The chick embryo studies showed knock down of HAND1 attenuates heart development [48]. GATA6 transcription factor expresses in the early mesoderm development during the heart specification and is reduced during the terminal differentiation [49]. The higher expression of these transcription factors in 14 days differentiation showed the specification of early embryonic development. A further differential expression with thalidomide demonstrates the causes of defects in early cardiac development which was consistent in an in vivo studies [15]. MSX1 plays an important role in cleft palate and craniofacial development, as demonstrated by MSX1-deficient mice [50]. The high gene expression levels of MSX1 in differentiated hESCs were repressed in the presence of thalidomide. These findings correlate well with cleft palate abnormalities observed in children's exposed to thalidomide [51,52].
Many molecular mechanisms have been proposed to explain the teratogenic effects of thalidomide including free radical oxygen species (ROS)-mediated oxidative DNA damage, oxidative stressmediated GSH depletion and NF-kB attenuation [16,20,21]. Most of the previous studies [21] mainly focus on the GSH depletion in response to thalidomide via ROS generation. Interestingly, in the present study we found reduced expression of GSTA ( Figure 5A, 5B, 5C) for thalidomide which is required for secondary cellular oxidative stress [32]. Oxidative stress affects basic cellular physiology, such as nuclear and cytoplasmic trafficking [33]. Any type of perturbation in the nuclear traffic machinery, a hallmark of oxidative stress, limits the efficiency of transport, which controls development, differentiation and transformation [53].
The down-regulation of RANBP1 (a RAN GTPase binding protein) ( Figure 5D) in 2DE analysis and gene expression pattern of importin-a (KPNA2) and nucleoporins (NUP98) was confirmed with real time and western blot analysis ( Figure 5E, 5F). Importinb and the adaptor importin-a are crucial components of the nuclear transport apparatus. Importins interact exclusively with the GTP bound form of the Ran GTPase. This process enables exportation of the importins from the nucleus in association with Ran-GTP, thereby regulating nucleocytoplasmic transport of cytosolic proteins through the nuclear pores [30]. A gradient of high Ran-GTP in the nucleus, versus high Ran-GDP in the cytoplasm, is established by cytoplasmic localisation of Ranbinding protein 1 (RanBP1) [30]. This results in a release of Ran-GTP from importins. It is well known that interference in the nucleocytoplasmic trafficking pathway affects embryonic development and differentiation [53], and it partially explains the teratogenic effects of thalidomide.
A recent publication [54] demonstrated that spontaneous differentiation of hESCs in the presence of toxicants dysregulated genes related to specific target tissues. Although this method can demonstrate the toxic effect of drugs for preliminary screening, a 7-day differentiation of the hESCs and a single-dose treatment assay limited the scope of the study. This study showed neither a critical dose for toxicity nor molecular mechanisms related to toxicants, including thalidomide. Differentiation of ESCs is a dynamic process, implicating the differential expression of thousand of genes in a very hierarchical manner. Therefore, missing the optimal time point or the appropriate concentration of the toxicant is critical for establishing a consistent and robust in vitro ESC-based toxicological model. Our independent time kinetic microarray study of hESC differentiation demonstrated that expression of embryonic development-related genes increased exponentially up to day 12. Following day 12, the global expression pattern was sustained for the remaining days (up to day 21) in hierarchical clustering ( Figure 1C). In previous study, we demonstrated that a 7-day differentiation of hESCs is not sufficient to monitor the developmental adverse effects of cytosine arabinoside [24]. In summary, this in vitro study demonstrated that the perturbation of limb and heart formation, the nucleocytoplasmic trafficking and GST expression (with respective markers at the genomic and proteomic level) may represent a critical embryotoxicity for thalidomide. One of the main advantages of this alternative approach is the application of human ESCs. Human ESCs better represent human cellular physiology, can improve the safety of drug testing and can partially explain the mechanism(s) of a specific toxicant. nuclear membrane in control is attenuated for 10 and 70 mM. The scale bar represents 20 mm. D, The genomic and mass spectrometry analysis showed down-regulation of RANBP1. E, Investigation of transporters genes for thalidomide treatment with RT-qPCR analysis (from independent experiment). The error bars represent the SEM from 3 technical replicates. F, Thalidomide mediated protein concentration changes were shown with immunoblotting analysis. G, The sub cellular localisation of RANBP1 is demonstrated with immunoflourescence analysis. The cytoplasmic enrichment of RANBP1 gradient is repressed with thalidomide treatment. The scale bar represents 20 mm. doi:10.1371/journal.pone.0044228.g005 Figure S1 The time kinetic experiment showed progressive expression of development related genes. To find the common genes among different days of differentiation intersect analysis was performed. We found regulated transcripts (n = 3), (p#0.05) constantly expresses after day 12. (TIF) Figure S2 A typical 2DE image of undifferentiated hESC (I) and 14-day-old EBs (II) are represented. The down (I) and up-regulated (II) spots are denoted in the gel pictures. IPisoelectric pH value. MW-molecular weight. The differentially regulated proteins (n = 3), (p#0.01) are from 3 biological replicates. (TIF) Figure S3 A typical 2DE image of control (I), thalidomide treatment (II) and percentage volume (III) for differentially regulated spots (n = 3), (p#0.01). The error bars represents SEM from 3 independent biological replicates. (TIF) Figure S4 Perturbation of thalidomide in early cardiac development. For time and dose response experiment, representative cardiac specific transcription factors were analysed with RT-qPCR Bar represents mean value from an independent experiment of 3 technical replicates (*p-value #0.01, thalidomidetreated vs untreated 14-days old EBs) and error bar shows SEM. Y-axis represents relative mRNA expression compared to control. X-axis shows temporal analysis for thalidomide treatment. (TIF)