Implications of TGFβ on Transcriptome and Cellular Biofunctions of Palatal Mesenchyme.

Development of the palate comprises sequential stages of growth, elevation, and fusion of the palatal shelves. The mesenchymal component of palates plays a major role in early phases of palatogenesis, such as growth and elevation. Failure in these steps may result in cleft palate, the second most common birth defect in the world. These early stages of palatogenesis require precise and chronological orchestration of key physiological processes, such as growth, proliferation, differentiation, migration, and apoptosis. There is compelling evidence for the vital role of TGFβ-mediated regulation of palate development. We hypothesized that the isoforms of TGFβ regulate different cellular biofunctions of the palatal mesenchyme to various extents. Human embryonic palatal mesenchyme (HEPM) cells were treated with TGFβ1, β2, and β3 for microarray-based gene expression studies in order to identify the roles of TGFβ in the transcriptome of the palatal mesenchyme. Following normalization and modeling of 28,869 human genes, 566 transcripts were detected as differentially expressed in TGFβ-treated HEPM cells. Out of these altered transcripts, 234 of them were clustered in cellular biofunctions, including growth and proliferation, development, morphology, movement, cell cycle, and apoptosis. Biological interpretation and network analysis of the genes active in cellular biofunctions were performed using IPA. Among the differentially expressed genes, 11 of them are known to be crucial for palatogenesis (EDN1, INHBA, LHX8, PDGFC, PIGA, RUNX1, SNAI1, SMAD3, TGFβ1, TGFβ2, and TGFβR1). These genes were used for a merged interaction network with cellular behaviors. Overall, we have determined that more than 2% of human transcripts were differentially expressed in response to TGFβ treatment in HEPM cells. Our results suggest that both TGFβ1 and TGFβ2 orchestrate major cellular biofunctions within the palatal mesenchyme in vitro by regulating expression of 234 genes.

Introduction

Cleft lip and/or palate is one of the most prevalent birth defects worldwide (1 in 800 live births; Schutte and Murray, 1999; Spritz, 2001), and is caused by failures in palate development. The formation of a continuous palate is a complex process composed of multiple steps, including palatal shelf growth, elevation, attachment, and fusion. Palatogenesis in the human spans from approximately gestational day 48 to 59 and the outgrowth of the secondary palate can generally be detected around day 49. During day 54–55, the palatine processes rapidly elevate, assuming a horizontal position which allows them to grow toward each other, attach, and fuse (Wyszynski, 2002). In general, with slight variation among strains, the stages of palatogenesis in mice [12.5–16.5 days post coitum (dpc)] are extremely similar and comparable to that of humans; therefore, mice have been used as a model to study human palate development (Ferguson, 1988). The failure of palatal shelves to grow and adhere after elevation is the most common type of cleft palate defect documented in murine models (Chai and Maxson, 2006). The palatal cellular components originate from the cranial neural crest (CNC)-derived palatal mesenchyme, concealed with a veneer of pharyngeal ectoderm-derived epithelium (Ito et al., 2003; Nakajima et al., 2010). A precise and time-sensitive regulation of various mesenchymal biofunctions, such as cellular movement, cell death (apoptosis), cell morphology, cell cycle progression, development, and growth and proliferation, is fundamental for the proper development of the palate. These cellular functions are coordinated by numerous genes encoding a range of growth factors, signaling mediators, transcriptional factors, cytokines, and extracellular matrix proteins (Richman and Tickle, 1989; Greene and Pisano, 2004, 2005). Therefore it is immensely important to explore the genes and the molecules that regulate the plethora of these biofunctions to understand cellular behavior during palatogenesis.

The TGFβ family consists of more than 30 ligand proteins, including activins, BMP, and TGFβ cytokines, regulating a wide variety of biological processes such as cellular development, morphology, movement, growth and proliferation, survival, mitotic regulation, apoptosis, and epithelial–mesenchymal transition (EMT). Although the three isoforms of TGFβ (β1, β2, and β3) are highly conserved between species (Rotzer et al., 2001) and share 71–76% sequence identity, these ligands have isoform-specific activities that cannot be compensated by other family members (Iwata et al., 2011). Based on knockout mouse models, TGFβ isoforms have been found to be essential for normal morphogenesis of the palate (Iordanskaia and Nawshad, 2011). TGFβ1 null mice are embryonically lethal and die before 11 dpc (Brunet et al., 1995), so its role in palate development cannot be evaluated. TGFβ2 knockout mice, which are also embryonic lethal at 18 dpc, have defects in their mandible and maxilla, with 23% of cases resulting in cleft palate (Sanford et al., 1997); whereas all TGFβ3 null mice develop cleft palate (Kaartinen et al., 1995; Proetzel et al., 1995) and die within 24 h after birth. Furthermore, TGFβ mutations and expression patterns have been shown to be associated with the occurrence of cleft lip and/or palate (Stoll et al., 2004). Mutations in TGFβ1 and TGFβ3, as well as their variants, are associated with cleft palate in humans (Lidral et al., 1998; Mitchell et al., 2001; Kim et al., 2003; Vieira et al., 2003; Rullo et al., 2006; Guo et al., 2010; Salahshourifar et al., 2011). Moreover, in vitro studies of human tissues showed that TGFβ1 and β3 are differently expressed and correlated with the cleft lip and/or palate phenotype (Bodo et al., 1999). Overall, these findings underscore the crucial function of TGFβ isoforms in the optimal regulation and completion of palate development.

The TGFβ isoforms are expressed in the early stages of mouse palate development (Fitzpatrick et al., 1990; Pelton et al., 1990; Gehris et al., 1991; Gehris and Greene, 1992). During the sequential steps of palatogenesis, TGFβ1 is expressed both in the epithelial and mesenchymal components of the palatal shelves (Fitzpatrick et al., 1990; Pelton et al., 1990). TGFβ2 is predominantly expressed in the mesenchymal cells, particularly immediately adjacent to the epithelium, with few epithelial cells also expressing TGFβ2 transcripts (Fitzpatrick et al., 1990; Pelton et al., 1990; Gehris and Greene, 1992). Intense and distinct localization of TGFβ3 has been detected in the medial edge epithelium (MEE) of the palate (Fitzpatrick et al., 1990; Pelton et al., 1990). In our lab, we have demonstrated that TGFβ3 is also expressed in the palatal mesenchyme, albeit at a lower level compared to the palatal epithelium (Unpublished data). It has also been shown that each TGFβ ligand can signal via different receptor complexes and downstream signaling molecules resulting in divergent cellular functions and behavior (Abbott and Pratt, 1988; Iwata et al., 2011). Therefore the isoforms of TGFβ may behave uniquely in the palatal mesenchyme compared to the palatal epithelium.

In this study, we investigated the crucial roles of TGFβ1, β2, and β3 on the regulation of palatal mesenchyme transcriptome and various cellular biofunctions, such as growth and proliferation, development, cell morphology, movement, cell cycle, and cell death. Using human embryonic palatal mesenchymal (HEPM) cells and bioinformatics tools, we analyzed how these isoforms regulate differential expression of transcripts and gene interaction networks within the palatal mesenchyme. Using microarray genechips, we found that expression of only 566 genes, which corresponds to >2% of the overall human transcriptome, were differentially expressed in TGFβ-treated HEPM cells with statistical significance; including candidate genes recognized as inducers of cleft palate either in human or mouse (EDN1, INHBA, LHX8, PDGFC, PIGA, RUNX1, SNAI1, SMAD3, TGFβ1, TGFβ2, and TGFβR1). We also found that exogenous TGFβ can regulate different molecules of the canonical TGFβ signaling pathway. Our results suggest that both TGFβ1 and TGFβ2, in concordance, regulate expression of 234 genes that govern cellular biofunctions which are crucial for the complete formation of the palate.

Materials and Methods

Cell culture and treatment

The HEPM cell line (CRL 1486) was obtained from the American Type Culture Collection (ATCC; Rockville, MD, USA) from a single human abortus at the time of palatal shelf elevation, but prior to epithelial contact. HEPM cells represent the undifferentiated fibroblast-like cells from embryonic palatal shelves. The cells were cultured at 37°C with 5% CO2. The culture medium consisted of DMEM (Mediatech, Manassas, VA, USA) supplemented with 10% (v:v) fetal bovine serum (FBS) and 1% (v:v) penicillin streptomycin. For the microarray and qPCR studies, following 24 h serum starvation, HEPM cells were treated either with 0.5% FBS (Control) or with recombinant TGFβ (R&D Systems, Minneapolis, MN, USA) at 10 ng/ml TGFβ1, TGFβ2, or TGFβ3 for 24 h in triplicate and harvested for RNA extraction.

RNA extraction

Total RNA from the TGFβ-treated HEPM cells was extracted using the RNeasy spin column RNA purification kit (Qiagen, CA, USA) according to the manufacturer’s instructions. RNA yield and integrity were initially evaluated using a NanoDrop 2000c spectrophotometer (NanoDrop, Wilmington, DE, USA). Samples with a concentration above 40 ng/ml, OD260:280 of 1.8–2.0, and OD260:230 above 1.8 were retained. Samples intended for microarray analysis were further evaluated using the Agilent 2100 Bioanalyzer and RNA 6000 Pico kit (Agilent, Santa Clara, CA, USA) according to the manufacturer’s instructions. Samples with electropherograms exhibiting sharp 18S and 28S rRNA peaks and no evidence of degradation were retained.

Microarray processing

RNA samples from each treatment group of triplicates (Total:12) were profiled using the Affymetrix GeneChip® Human Gene 1.0 ST Array (Affymetrix Inc., Santa Clara, CA, USA) representing 28,869 genes. In this chip, each of the genes was represented on the array by approximately 26 probes spread across the full length of the gene, providing a more complete and accurate interpretation of gene expression than 3′ based expression array designs (Pradervand et al., 2008). The microarray procedure was performed according to the standard Affymetrix gene chip analysis protocol at the UNMC Microarray Core. Three hundred nanograms of total RNA was reverse transcribed and converted to double stranded cDNA using the Ambion WT Expression kit (Ambion, Austin, TX, USA). The samples were subsequently subjected to fragmentation and biotin end-labeling using the Affymetrix WT Labeling Kit and hybridized per the manufacturers suggestions to the Affymetrix Human Gene ST 1.0 arrays. The wash, stain, and scan procedures were carried out on the Affymetrix Hybridization Oven, FS 450 wash station, and the 7000G scanner per the manufacturer’s recommendations (Affymetrix, Santa Clara, CA, USA).

qPCR

Synthesis of cDNA prior to qPCR was performed using the SuperScript® III First Strand SuperMix (Invitrogen, Carlsbad, CA, USA) with 5 μg input total RNA per 25 μl reaction according to the manufacturer’s instructions. Negative control samples were performed in parallel by omitting reverse transcriptase. The sequences of primers for genes of interest (Table A1 in Appendix) were designed using the PrimeTime qPCR Assay [Integrated DNA Technologies (IDT), Iowa City, IA, USA) and synthesized by Eurofins MWG Operon (Huntsville, AL, USA). All qPCR experiments were performed using cDNA synthesized from the same batch and starting amount of total RNA. Negative controls containing no cDNA template were included for each gene within each qPCR run. To avoid variation in amplification conditions across runs, reactions for all experimental conditions (i.e., all untreated and TGFβ-treated HEPM cells) for each gene of interest were performed as triplicates in the same 96-well plate. Cycling conditions were designed according to the RealMasterMix SYBR ROX handbook (5 PRIME, Gaithersburg, MD, USA) as initial activation at 95°C for 10 min; 40 cycles of 95°C for 15 s, 60°C for 1 min. Amplification specificity for each gene was confirmed by a single distinct melting curve. Relative quantification of gene expression was determined by normalization of the samples with the internal reference gene, GAPDH. Data from the triplicates was evaluated and analyzed for significance by the pair-wise fixed reallocation randomization test using the REST 2009 software (Pfaffl et al., 2002). A p-value of ≤0.05 was considered significant.

Table A1

Primers for q-RTPCR.

Gene name	Entrez ID	Fw primer	Rv primer
AMIGO2	NM_001143668.1	CAACATCACCAGCATTTCCAC	TCACCGTCTTCAGCTTATTGG
NOX4	NM_016931	TCACAGAAGGTTCCAAGCAG	ACTGAGAAGTTGAGGGCATTC
TNFAIP6	NM_007115.3	AATACAAGCTCACCTACGCAG	GGTATCCAACTCTGCCCTTAG
DACT1	NM_016651.5	TTGGAGGAGAACATCTTGCTG	GTCTTTTCTACATCCAGTCTCAGG
ANKRD1	NM_014391.2	GGTGAGACTGAACCGCTATAAG	GGCTGTCGAATATTGCTTTGG
PMEPA1	NM_020182	GCAAACGCTCTTTGTTCCAG	ACCATCACCATCATCACCAC
MSC	NM_005098.3	CGCTATGAGAACGGCTACG	CCCATCAAGTGAGTTCCAGTC
EDN1	NM_001955	CTTCGTTTTCCTTTGGGTTCAG	GCTCAGCGCCTAAGACTG
VDR	NM_000376.2	CACTATTCACCTGCCCCTTC	CTTCCTCTGCACTTCCTCATC
ITGB3	NM_000212.2	CAAGTGTGAATGTGGCAGC	TTTTCGTCATGTAGGGCTCC
KCNH1	NM_172362.2	GCCTTCTCCCATTCCTTCTC	CCTCATTCTTTCGTTTCATGCG
HAPLN1	NM_001884.3	AGTGTGAGGTGATTGAAGGATTAG	CTGCGCCTCGTGAAAATTG
NEDD9	NM_006403.3	TGTAGGAAAACGGCTCAACC	CCCTGTGTTCTGCTCTATGAC
GDF5	NM_000557.2	ACAGAAAGGGAGGCAACAG	CTTCTCCAGGGCACTAATGTC
SECTM1	NM_003004.2	GACACCAGAGAAATAACAGACAAG	GTACCAGGCGAACATGACC
NPTX1	NM_002522.3	GCCCTATCACCCCATCAAG	TTGCGGTCCCAGATGTTG
SMAD3	NM_005902	TCCATCCCCGAAAACACTAAC	CATCTTCACTCAGGTAGCCAG
OAS1	NM_016816	CATCTGTGGGTTCCTGAAGG	GAGAGGACTGAGGAAGACAAC
TBC1D8	NM_001102426	CCCAGGTTTATTCTCCCATAGC	AGCCTTTACCTTCCCTTTGAC
TRERF1	NM_033502.2	GTGATGGTTGCTCTGGAAATG	TCTTTGCTGTAAGTGGCTAGTG
VSIG1	NM_182607	GCTCTCAACTAACCTCCACAC	AGATCCAACAGTCACGTTCAC
TGFBR3	NM_003243.4	CGGAAACATCACCTTCAACATG	CCCAGTTCTTGTTCAGCCTTAG
CA12	NM_001218	CAGGTCCAGAAGTTCGATGAG	CAATACAGATGCCAAGAATGCC
ARHGAP28	NM_001010000	AGTAAAAGGACGAGACAATGGG	ATTCCTTCAGATTCCAGACCTG
HGF	NM_000601.4	GCTATACTCTTGACCCTCACAC	GTAGCCTTCTCCTTGACCTTG
NR4A2	NM_006186	CTGGCTGTTGGGATGGTC	TGTGGGCTCTTCGGTTTC
IFIT2	NM_001547.4	GAAGATTTCTGAAGAGTGCAGC	CCTCCATCAAGTTCCAGGTG

Data analysis: Normalization and differential gene expression

Scanned array images were analyzed by dChip applying a smoothing spline normalization method prior to obtaining model-based gene expression indices, a.k.a. signal values. There were no outliers identified by dChip, so all samples were carried on for subsequent analysis. When comparing two groups of samples to identify enriched genes in a given group, we used the lower confidence bound (LCB) of the fold change (FC) between the two groups as the cut-off criteria. If 90% LCB of FC between the two groups was above 1.2, the corresponding gene was considered to be differentially expressed. LCB is a stringent estimate of the FC and has been shown to be the better ranking statistic (Li and Wong, 2001). By use of LCB, we can be 90% confident that the actual FC is some value above the reported LCB. Researchers demonstrated that (Ramalho-Santos et al., 2002) selecting genes that have a LCB above 1.2 most likely corresponds to genes with an “actual” FC of at least three in gene expression. For all genes scored, the FC was calculated by dividing the experiment value by the control value. If this number was less than one, the (negative) reciprocal is listed (e.g., 0.75 or a drop of 25% from control is reported as −1.25 FC). The reported FCs are the average of the independent experiments. All of the significantly changed genes (FC >1.2-fold down or up) were selected for the biofunctional analysis. p-Values less than 0.05 indicated a statistically significant, non-random association. The complete list of genes passing dChip’s statistical filter for each treatment with their respective FCs and p-values was deposited to the gene expression omnibus (GEO) of NCBI (Accession number: TBD). To compare the results of our multiple microarray experiments, we drew Venn diagrams which depicted the relations among the TGFβ treatment groups as intersections or uniqueness using the GeneVenn (Pirooznia et al., 2007).

Clustering and functional analysis

Clustering of samples and genes was performed using the unweighted pair group method with arithmetic-mean (UPGMA), a hierarchical clustering technique used to construct a similarity tree, and principal components analysis (PCA; Sneath, 1973; Handl et al., 2005). Heatmaps with dendrograms are constructed using Bioconductor of R statistical software (Gentleman et al., 2004; Reimers and Carey, 2006). Samples and genes were clustered using the normalized and modeled expression values obtained by dChip analysis. The expression data matrix was row-normalized for each gene prior to the application of average linkage clustering and the Pearson’s correlation was used as the distance measure. In the heatmaps, each column corresponds to the averaged expression profile of a treatment in triplicate, and each row corresponds to an mRNA. The upper portion (green-to-red) represents downregulated transcripts and the lower portion (red-to-green) represents upregulated transcripts. The increasing intensities of red signify a higher expression in the given sample of a specific mRNA, whereas the increasing intensities of green indicate a lower expression of mRNA and black indicates mean level expression. Functional analysis was comprised of finding gene ontology (GO) categories in the gene clusters of interest that warrant further investigation. GO annotations are used for the biological interpretation of detected clusters (Ashburner et al., 2000).

Biological interpretation and pathway analysis

The effects of TGFβ on cellular functions of HEPM cells were examined in the context of differentially expressed human genes using the ingenuity pathway analysis (IPA; Ingenuity Systems, CA, USA) a web-delivered application used to discover, visualize, and explore relevant networks. Following statistical analysis and filtering of the microarray data using dChip, Affymetrix probe identifiers, and fold-values of each treatment group were uploaded to IPA as a dataset. Each identifier was mapped to its corresponding gene object in the ingenuity knowledge base (IKB). Interactions were then queried between these datasets and all other gene objects stored within IPA to generate a set of direct interaction networks that were overlapped. The significant genes were categorized, compared to functional categories in the IPA database, and ranked according to their p-values. p-Values less than 0.05 indicate a statistically significant, non-random association between a set of significant genes and a set of all genes related to a given function in the IKB (Mori et al., 2009). Through the assessment of differentially expressed genes, cellular processes that are most significantly altered in our dataset were: cellular growth and proliferation, development, morphology, movement, cell cycle, and cell death. Using these IPA analyses, molecules within the cellular biofunction categories were filtered and imported to Excel (Microsoft, Seattle, WA, USA) for further graphical analysis of the networks.

Results

TGFβ regulates differential expression of more than 2% of the palatal mesenchyme transcriptome

We performed gene expression profiling on the Affymetrix GeneChip® Human Gene 1.0 ST Array representing 28,869 genes. We explored the alterations in gene expression in HEPM cells extracted from human embryonic palatal shelves in response to 10 ng/ml TGFβ1, TGFβ2, and TGFβ3 for 24 h. Statistically significant altered genes were filtered using dChip through normalization and modeling analysis. Out of 28,869 genes represented within the GeneChip®, 301 genes demonstrated upregulated and 200 genes demonstrated downregulated expression in response to TGFβ1 treatment. Similarly, 190 genes were upregulated and 77 genes were downregulated in response to TGFβ2 treatment. Interestingly, only four genes demonstrated increased and four genes demonstrated decreased expression patterns in response to TGFβ3 treatment. The numbers of differentially expressed genes, indicating individual and overlapping aggregations of TGFβ isoforms, are depicted in a Venn diagram (Figure 1).

Figure 1

Venn diagrams. Multiple microarray experiments were designed to detect differential expression of transcripts with TGFβ treatment and were compared using Venn diagrams. List of differentially expressed genes with their accession IDs were uploaded to GeneVenn, query was run, and the detected numbers and gene lists were exported to MS Excel. Graphical figure is illustrated using the Adobe Photoshop (Adobe Inc, CA, USA). Relations among the list of altered genes within TGFβ treatment groups are depicted as intersections or uniqueness. A significant overlap of differentially expressed genes was observed between the TGFβ1 and TGFβ2 treatment groups both in the upregulated (A) or downregulated (B) gene lists.

The genes (n = 27) displaying the most significant changes in expression patterns in response to TGFβ1 and TGFβ2 stimulation with their FCs and p-values are shown in Table 2. According to dChip normalization and modeling, the level of differential gene expression by TGFβ3 was detected to be statistically non-significant; therefore TGFβ3 has been excluded from the FC analysis. Transcripts, which were differentially expressed both by TGFβ1 and TGFβ2, are highlighted in the table (Table 2). Graphical representation of FC values detected by microarray analysis is shown in Figure 2A. In order to validate the gene expression profiling results of these most significantly up- and downregulated transcripts, expression levels of 27 genes were quantified by qPCR and compared to the levels determined by microarray (Figure 2B). Expression profiles of all 27 genes tested were found to be in agreement between the two analytical assay results. We identified a strong correlation of gene expression pattern between microarray and qPCR data in terms of up- and downregulation, but not relative level of expression.

Table 2

The most significantly altered genes within the palatal mesenchyme with TGF-β treatment.

TGF-β1			TGF-β2			TGF-β3
TOP-10 UPREGULATED GENES
Gene name	FC	p-Value	Gene name	FC	p-Value	Gene name	FC	p-Value
AMIGO2	22.89	4.88E-03	AMIGO2	19.66	1.30E-05	KIAA1654	1.56	2.18E-02
NOX4	13.87	2.28E-02	NOX4	11.12	2.72E-02	NTM	1.47	3.91E-02
TNFAIP6	9.36	1.02E-02	ANKRD1	10.33	2.38E-03	LAMA1	1.46	3.58E-02
DACT1	8.02	7.88E-03	KCNH1	7.76	2.14E-02	LRP11	1.41	3.78E-02
ANKRD1	6.52	4.48E-02	EDN1	6.48	1.65E-02
PMEPA1	5.69	1.56E-03	DACT1	6.47	2.54E-03	Common upregulated
MSC	5.55	2.01E-02	PMEPA1	5.55	2.23E-02
EDN1	5.23	1.44E-02	HAPLN1	4.89	1.37E-02
VDR	5.19	1.36E-03	VDR	4.75	6.66E-04
ITGB3	5.19	1.04E-02	NEDD9	4.54	4.67E-03
TOP-10 DOWNREGULATED GENES
GDF5	-4.08	2.07E-02	TBC1D8	-2.570	6.70E-03	GEN1	-1.40	2.75E-02
>SECTM1	-3.33	2.44E-02	>SECTM1	-2.540	4.37E-02	POLK	-1.49	1.40E-02
>NPTX1	-3.25	1.61E-02	>SMAD3	-2.510	1.63E-02	FLVCR2	-1.52	4.42E-02
>SMAD3	-3.24	1.09E-02	>OAS1	-2.290	3.95E-02	CPA4	-1.70	1.40E-02
>OAS1	-3.11	2.27E-02	ARHGAP28	-2.230	4.42E-02
>TBC1D8	-3.06	4.50E-03	>NPTX1	-2.220	2.59E-02	Common downregulated
TRERF1	-2.75	3.35E-03	>VSIG1	-2.190	3.11E-02
>VSIG1	-2.75	3.00E-02	HGF	-2.150	9.75E-03
TGFBR3	-2.62	2.82E-03	NR4A2	-2.020	4.56E-02
CA12	-2.59	7.71E-03	IFIT2	-1.970	3.39E-02
GENES SUSCEPTIBLE FOR CRANIOFACIAL DEFORMITIES
EDN1	5.23	1.44E-02	EDN1	6.48	1.65E-02
INHBA	2.54	7.20E-03	INHBA	2.21	1.62E-02
LHX8	1.62	1.84E-02	LHX8	NS	NS
PDGFC	2.21	6.21E-03	PDGFC	2.06	5.95E-03
PIGA	1.51	1.71E-02	PIGA	NS	NS
RUNX1	1.78	1.06E-03	RUNX1	1.68	1.39E-03
SMAD3	−3.24	1.09E-02	SMAD3	−2.51	1.63E-02
SNAI1	1.75	1.39E-02	SNAI1	1.63	4.50E-02
TGFβ1	1.91	3.26E-03	TGFβ1	1.59	2.29E-03
TGFβ2	NS	NS	TGFβ2	1.78	4.22E-02
TGFβR1	1.49	2.05E-02	TGFβR1	−1.81	7.33E-03

NS, Non-significant expression detected by dChip. Red (A): Genes upregulated in response to both TGFβ1 and TGFβ2. Green (B): Genes downregulated in response to both TGFβ1 and TGFβ2.

Figure 2

Expression levels of the most significantly up- and downregulated genes based on microarray and qPCR analyses. The graphical representation of genes (n = 27) displaying the most significant changes in expression patterns in response to TGFβ stimulation with their FC values based on microarray analyses were shown (A). Only TGFβ1 and TGFβ2 treatments were represented in (A) due to their statistically significant p-values detected by dChip. (B,C) Gene expression profiles of all 27 genes were independently validated using qPCR and data was analyzed using REST 2009 (Pfaffl et al., 2002). Expression variation for each gene was visualized in a whisker-box plot. Transcripts were sorted alphabetically.

Furthermore, direct or indirect relations of these 27 genes with TGFβ signaling were detected using IPA. Among the significantly upregulated transcripts, we detected that ANKRD1, EDN1, IGFBP3, ITGB3, NOX4, NEDD9, PMEPA1, TNFAIP6, and VDR were directly related with TGFβ; while AMIGO2 and MSC were indirectly associated with TGFβ via other molecules (Figure 3A). Within the significantly downregulated gene group HGF, NR4A2, SMAD3, and TGFBR3 presented a direct relation with TGFβ; while transcripts of TBC1D8, GDF5, IFIT2, and OAS1 presented an indirect association with TGFβ via other molecules (Figure 3B).

Figure 3

Direct and indirect relations of most differentially expressed genes with TGFβ. Affymetrix probe identifiers and FC values of significantly altered transcripts (n = 27) were uploaded to IPA and each identifier was mapped to its corresponding gene object in the IKB. Interactions were then queried between these gene objects and all other gene objects stored within IPA to generate a set of direct interaction networks that were merged. Networks with the highest number of molecules were selected, overlaid with TGFβ molecules, and organized using Path Designer of IPA. Upregulated transcripts were highlighted with orange (A); and downregulated transcripts were highlighted with green (B). Types of molecules were annotated in the legend in the box.

We focused on the differential expression of genes involved in cellular biofunctions with a criterion of FC >1.2 using the IPA software. While several genes overlapped based on cellular functions, we were able to cluster them according to their GO annotations as: cell cycle, cell death, cell morphology, growth and proliferation, development, and movement (Table 1). The individual heatmaps were generated according to hierarchal clustering of genes which are involved in the cellular processes (Figure 4). According to IPA analysis, there were 66 differentially expressed genes detected to have a role in cellular movement; 48 were upregulated, and 18 were downregulated in response to TGFβ treatment. Eighty genes were found to play a role in apoptosis (cell death); 55 demonstrated increased expression, while 25 displayed lower expression levels). Fifty genes function in cell morphology; 12 were downregulated and 38 were upregulated. Cellular development included 91 genes, among which 53 were upregulated and 38 were downregulated. The Cell cycle GO was represented by 37 genes; 20 were upregulated and 17 were downregulated. The GO cluster containing the highest number of differentially expressed genes (n = 93) compared to other cellular biofunction clusters was the cellular growth and proliferation GO, including 57 upregulated and 36 downregulated transcripts.

Table 1

Differentially expressed molecules clustered into GO of cellular biofunctions.

Function	Molecules
Cell cycle	TGFBR1, PCGF2, SGK1, SMAD3, S100A4, MKI67, DDB2, LIMK1, CAMK2D, BOP1/LOC727967, TGFB1, HIPK2, POLK, NEDD4L, ALDH3A1, PDGFRB, TP53, STAT6, DLG1, AMACR, NOX4, mir-199, SNAI1, STK38L, ANLN, JUNB, INHBA, KAT2B, IGF2, IRS1, BHLHE40, MBD4, SIK1, KIF11, RUNX1, SULF1, TPM1, GADD45B, USP2, SEPT9, FAS, EIF4EBP1, KSR1, CDC25B, EDN1, ASPM, CCDC99, MPHOSPH6, HGF, FOXO3, TOP2A, NEK7, RPS6KA2, NUSAP1, TNC, PPARD, MGMT, GNAI1, SMAD7, XBP1, GRB10, IVNS1ABP, CENPI, CSF1, KIF20B, CAT, KLF5, ZEB2, SKIL, NCOR2, mir-27, NEDD9, IRAK4
Cell death	CTGF, TGFBR1, PCGF2, APH1B, SGK1, TGFBR3, DDB2, NPTX1, BNIP3, SGPL1, GDF5, PRDM1, SERPINE1, HIPK2, ALDH3A1, MVP, ITGB5, PDGFRB, TP53, STAT6, mir-199, NOVA1, NDRG1, DDIT4, FOXP1, SLC2A1, SNAI1, JUNB, CDH2, CD9, IRS1, BHLHE40, MBD4, SFRP1, FLNB, TCF4, GADD45B, ZFP36, IL17RD, PTPN13, HYOU1, TAOK3, EIF4EBP1, KSR1, EFNA5, RYBP, TOP2A, PARP4, RPS6KA2, TLR3, VDR, FAP, OAS1, UCP2, TNC, TGFBI, MGMT, XBP1, FAIM2, OBFC2A, MMP10, TRIB3, P2RX7, PTGES, SEMA3A, KLF7, CSF1, BMF, ZEB2, CTH, NCOR2, NEDD9, SGCD, PTGR1, AMIGO2, IRAK4, CHRNA1, EPHB2, SMAD3, DPYSL4, UNC5B, S100A4, RBP1, CAMK2D, TGFB1, POLK, SMAD1, RAI14, ITGB3BP, PTPRE, NOX4, PDE4B, CDCA2, INHBA, TMEM57, ITPK1, KAT2B, ANKRD1, IGF2, NUAK1, ALDH3B1, JAG1, SIRPA, PMEPA1, TNFAIP8, RUNX1, KDM3A, SULF1, TPM1, DPYSL3, KLF10, LIMS1, USP2, FAS, PRDX6, CDC25B, SYVN1, EDN1, HGF, FOXO3, NEK7, EYA4, ADRB2, SCN2A, PPARD, BGN, SMAD7, IKBKE, DLX2, ITGB3, GRB10, TSHZ3, IVNS1ABP, EBF1, B4GALT1, LAMA4, KLF5, CAT, SH3RF1, SKIL, ATP2B4
Cellular growth and proliferation	TGFBR1, CTGF, PCGF2, TGFBR3, DDB2, MKI67, BNIP3, SGPL1, MFAP2, GDF5, BOP1/LOC727967, AFAP1, OSMR, PRDM1, HIPK2, SERPINE1, SPOCK1, ALDH3A1, ITGB5, MVP, PDGFRB, TP53, STAT6, AMACR, FOXP1, SLC2A1, FBLN5, NDRG1, ANLN, GRK5, JUNB, VLDLR, CDH2, CD9, IRS1, GLUL, SFRP1, TRIM22, TNFAIP6, TCF4, SLC4A2, GADD45B, ZFP36, PTPN13, BMPR2, SEPT9, EIF4EBP1, KSR1, PTPRJ, TLR3, VDR, IRX3, GLMN, TNC, TGFBI, MGMT, XBP1, TOB1, PDLIM2, P2RX7, SLC29A1, SEMA3A, PTGES, NOV, CSF1, KIF20B, BMF, TBC1D8, ZEB2, CTH, NCOR2, NEDD9, IRAK4, CHRNA1, EPHB2, SMAD3, UNC5B, S100A4, mir-221, PBX1, RBP1, TSC22D3, LIMK1, STARD10, CAMK2D, TGFB1, RORA, ITGA11, SLC7A5, SMAD1, PTX3, DLG1, PFKFB3, PTPRE, NOX4, STK38L, INHBA, KAT2B, ANKRD1, IGF2, MAPRE2, JAG1, SIRPA, KIF11, PMEPA1, TNFAIP8, RUNX1, TPM1, KLF10, PDGFC, PLEKHO1, FAS, CDC25B, AEBP1, EDN1, HGF, FOXO3, ENPP2, ADRB2, PPARD, BGN, SMAD7, IKBKE, MSI2, EDF1, PIGA, AFF1, GRB10, ITGB3, IVNS1ABP, B4GALT1, LAMA4, KLF5, CAT, SKIL, C5orf13, SHMT2
Cell morphology	EPHB2, ANK1, BNIP3, ITGA11, TGFB1, PRDM1, CDC42EP1, SERPINE1, NEDD4L, ITGB5, TP53, DLG1, NOX4, INHBA, IGF2, CD9, IRS1, CDC42SE1, JAG1, SIRPA, TPM1, FLNB, SULF1, SLC4A2, LTBP2, FAS, KSR1, EDN1, EFNA5, HGF, FOXO3, VDR, ADRB2, TNC, TGFBI, PPARD, ADCY3, SMAD7, TOB1, P2RX7, PVR, PIGA, EDF1, ITGB3, SEMA3A, CAP2, CSF1, LAMA4, CAT, NEDD9
Cellular development	GFBR1, CTGF, PCGF2, TGFBR3, DDB2, BNIP3, SGPL1, GDF5, OSMR, PRDM1, HIPK2, SERPINE1, ITGB5, PDGFRB, TP53, AMACR, STAT6, FOXP1, FBLN5, NDRG1, SNAI1, JUNB, VLDLR, CDH2, CD9, IRS1, BHLHE40, GLUL, SFRP1, FLNB, TCF4, MSC, SLC4A2, GADD45B, ZFP36, PTPN13, NID2, IL17RD, BMPR2, SEPT9, SMURF1, KSR1, EFNA5, TOP2A, DAGLB, TLR3, VDR, IRX3, GLMN, TNC, TGFBI, NRP2, MGMT, XBP1, TOB1, TRIB3, P2RX7, PVR, SEMA3A, NOV, CAP2, CSF1, KIF20B, BMF, ZEB2, LBH, NCOR2, IRAK4, EPHB2, SMAD3, UNC5B, mir-221, S100A4, PBX1, SYNE1, RBP1, BPGM, TSC22D3, ANK1, LIMK1, RORA, TGFB1, ITGA11, LHX8, SMAD1, DLG1, NOX4, LMO4, INHBA, IFRD1, KAT2B, IGF2, BHLHE41, CDC42SE1, JAG1, SIRPA, PMEPA1, TNFAIP8, RUNX1, KDM3A, TPM1, GLIS3, KLF10, FAS, PLEKHO1, CDC25B, EDN1, ASPM, HGF, FOXO3, MKL2, PHGDH, LONP1, ADRB2, VEZF1, PPARD, BGN, SMAD7, IKBKE, DLX2, MSI2, EDF1, PIGA, GRB10, ITGB3, AFF1, IVNS1ABP, EBF1, B4GALT1, LAMA4, CAT, KLF5, SKIL, C5orf13
Cellular movement	CTGF, TGFBR1, TGFBR3, NPTX1, SGPL1, GDF5, SCHIP1, AFAP1, SERPINE1, HIPK2, ITGB5, PDGFRB, TP53, STAT6, mir-199, FBLN5, SLC2A1, NDRG1, SEMA5A, SNAI1, ANLN, JUNB, VLDLR, CDH2, CD9, IRS1, SFRP1, FLNB, LTBP2, BMPR2, SEPT9, PARP9, PTPRJ, PCDH10, EFNA5, PPAP2B, TOP2A, CHST10, VDR, TLR3, FAP, UCP2, TNC, NRP2, MMP10, LRRC15, PDLIM2, PVR, P2RX7, SEMA3A, PTGES, NOV, CSF1, KIF20B, ZEB2, ANXA3, NEDD9, IRAK4, EPHB2, SMAD3, UNC5B, S100A4, LIMK1, TGFB1, NEDD4L, NOX4, STK38L, PDE4B, INHBA, IGF2, NUAK1, JAG1, SIRPA, COL7A1, TNFAIP8, SULF1, TPM1, LIMS1, FAS, SP100, CDC25B, EDN1, HGF, FOXO3, ENPP2, ADRB2, NUSAP1, PPARD, BGN, SMAD7, GNAI1, DLX2, ITGB3, IVNS1ABP, EBF1, B4GALT1, CAT, KLF5, ASAP2, C5orf13, ATP2B4

Figure 4

Hierarchical clustering (heatmap) analysis of differentially expressed transcripts. Altered genes were clustered based on cellular biofunctions and heatmaps were constructed using Bioconductor. Coexpressed groups of genes were illustrated with dendrograms and progressively and distinctly up- or downregulated genes were ranked accordingly. Each column corresponds to the expression profile of a treatment (either TGFβ or Control) in triplicate, and each row corresponds to a transcript. The increasing intensities of red signify a higher expression in the given sample of a specific mRNA, whereas the increasing intensities of green indicate a lower expression of mRNA and black indicates mean level expression. The list of genes for each particular function (A–F) was shown in Table 1.

A representation of the network-based interactions of the differentially expressed genes, according to their cellular biofunctions, and their molecule type, is shown in Figure 5. Due to the high number of genes in each biofunction, we preferred the genes with an FC >2.0 for network interactions. In the cell cycle network, only three genes (SMAD3, let-7, and EPGN) were downregulated, and seven genes (mir-199, IGFBP3, SKIL, NOX4, IVNS1ABP, EDN1, and NEDD9) were upregulated. Within the cellular movement network there were 4 significantly downregulated genes (SMAD3, SERPINB2, GDF5, and NPTX1) and 15 upregulated genes (ADAM19, VDR, VLDLR, NUAK1, NEDD9, LRRC15, mir-199, IGFBP3, EDN1, MAP2, ADAM12, LTBP2, IVNS1ABP, NOX4, and ITGB3) in response to TGFβ treatment. In the network of cellular growth and proliferation, we detected 8 genes (TRA2A, EPGN, GDF5, SERPINB2, SMAD3, let-7, GPAM, and mir-218) which were downregulated, while 17 genes (KCNH1, ITGB3, PMEPA1, NOX4, IVNS1ABP, VLDLR, IGFBP3, ADAM12, EDN1, AK4, VDR, ANKRD1, TGFβI, NEDD9, SKIL, TNFAIP6, and SPOCK1) were upregulated. Downregulated genes belonging to the cell morphology GO network include NPTX1 and SMAD3; and upregulated genes (n = 11) include LTBP2, VDR, ASPN, TGFβI, IGFBP3, EDN1, MAP2, VLDLR, NEDD9, NOX4, and ITGB3. In the cellular development network, 6 genes (mir-218, TRA2A, GDF5, SERPINB2, SMAD3, and let-7) were downregulated and 12 genes (EDN1, PMEPA1, SKIL, VDR, ASPN, IVNS1ABP, TGFβI, VLDLR, NOX4, IGFBP3, ITGB3, and KCHN1) were upregulated. Within the cell death (apoptosis) GO network, 5 genes (NPTX1, GDF5, SERPINB2, SMAD3, and GPAM) were downregulated and 14 genes (NUAK1, SKIL, VDR, AMIGO2, ANKRD1, TGFβI, ADAM12, mir-199, IVNS1ABP, IGFBP3, NEDD9, NOX4, EDN1, PMEPA1, and ITGB3) were upregulated.

Figure 5

Network of cellular biofunctions and corresponding transcripts. Networks were constructed with the IPA software using the selected genes matching the criterion of FC >2.0 and involved in cellular biofunctions. Several differentially regulated genes from the study were used to construct a gene association map for predicting various cellular and molecular events operating within the palatal mesenchyme (A–F). Red symbols specify upregulated expression of genes, whereas green symbols indicate downregulated genes. The color darkness represents the FC intensity.

Expression levels of orofacial development genes are also affected by TGFβ treatment

A table of genes vital to the morphogenesis of orofacial development was constructed based on the most recent literature review (Table A2 in Appendix; Wilkie and Morriss-Kay, 2001; Lidral and Moreno, 2005; Gritli-Linde, 2007; Dixon et al., 2011). Genes from the table of cellular biofunctions (Table 1) were merged with the Table A2 in Appendix to filter the differentially expressed genes according to putative candidates for palatal deformities. Some of the genes listed were identified as being important for palatogenesis in mouse, but not directly associated with human cases of cleft palate and thus remain “putative” cleft palate genes in the latter, as labeled with asterisk (*) in the Table A2 in Appendix. Overlapping genes (n = 11) were presented in Table 2 with their FC and p-values. Among these genes, EDN1, INHBA, TGFβ1, TGFβ2, RUNX1, LHX8, SNAI1, PDGFC, and PIGA were upregulated; while only SMAD3 was downregulated in response to TGFβtreatment. Intriguingly, TGFβR1 was upregulated with TGFβ1 treatment and downregulated with TGFβ2 treatment; which correlates with isoform-specific receptor requirement in TGFβ signaling (Rojas et al., 2009).

Table A2

Genes susceptible for cleft palate.

Abca4	Fst	Mthfr	Tbx1
Acvr2*	Gabrb3	Myf5	Tbx22
Acvr2a	Gad1*	Myod1*	Tbx3*
Apaf1*	Gad2*	Nat2	Tcof1
Bmp	Gas1	Ofd1	Tgfα
Bnc2*	Gli	Pax9	Tgfβ1*
Cask*	Gli3	Pdgfc*	Tgfβ2*
Cdkn1c*	Gpc3	Pdgfr-α*	Tgfβ3
Col11a1	Hic1*	Pds5*	Tgfβr1*
Col11a2	Hoxa2	Pds5b*	Tgif1
Col2a1	Hspg2	Piga*	Tgif2
Crk*	Hyal1	Pitx1	Tp63
Cspg	Inhba*	Pitx2	Tshz1*
Dhcr7*	Insig1*	Prrx1*	Vax1*
Dlx1*	Insig2*	Ptprf*	Vcan*
Dlx5*	Integrin*	Ptprs*	Viaat*
Dtdst	Irf-6	Pvrl1	Vlk*
Edn1*	Itgav*	Ryk	Wnt
Edna*	Itgb6*	Satb2	Wnt9b*
Egf	Itgb8*	Shh	Zeb1
Egfr	Jag2	Shox2*
Esr1	Lgr5*	Sim2*
Fgf	Lhx8*	Six3
Fgf-10*	Mafb	SMAD2
Fgfr	Meox1*	SMAD3
Fgfr-2*	Mid1	SMAD4
Folr1*	Mmp	Snail 1
Foxc2	Mn1*	Snail 2*
Foxe1	Mnt*	Snail 3*
Foxf2	Msx1	Sox9

*Putative genes, which may cause cleft palate in human.

Differentially expressed putative cleft palate genes in response to TGFβ1 and TGFβ2.

Ingenuity pathway analysis, a knowledge based program to generate relevant and interacting biological networks, was used to determine how the differentially expressed overlapping genes might interact to facilitate cellular biofunctions. A direct interaction network of the orofacial deformity susceptibility genes, overlaid with clustered cellular biofunctions crucial for development and maturation of the palatal mesenchyme was generated and is presented in Figure 6. Finally, the localization of differentially expressed genes within the canonical TGFβ pathway was illustrated by using the IPA Path Designer in Figure 7. The overlapping genes were: TGFβ1, TGFβ2, INHBA (activins/inhibins), BMPR2, SMAD7, SMURF1, SMAD 1/5/8, SMAD2/3, VDR, and PAI-1; suggesting that not all of the TGFβ pathway molecules are significantly regulated by exogenous TGFβ treatment. The list of genes regulated by differentially expressed TGFβ signaling molecules and their level of altered expression are represented at Table 3. The full list of these genes, including non-significantly altered, is provided in Table A4 in Appendix.

Figure 6

Merged interaction network of cleft palate susceptibility genes and biofunctions. Affymetrix probe identifiers and fold-values were uploaded to IPA and each identifier was mapped to its corresponding gene object in the IKB. Interactions were then queried between these gene objects and all other gene objects stored within IPA to generate a set of direct interaction networks that were merged. Cellular functions that directly interact with differentially expressed putative candidate genes for cleft palate in response to TGFβ treatment of HEPM cells were included in the network. The 11 transcripts meeting the 0.05 criteria of FDR in response to TGFβ treatment were analyzed for the direct interaction network using IPA. Upregulated genes (EDN1, INHBA, TGFβ1, TGFβ2, RUNX1, LHX8, SNAI1, PDGFC, and PIGA) were highlighted with orange and downregulated gene (SMAD3) was highlighted with green. TGFβR1 was upregulated with TGFβ1 treatment and downregulated with TGFβ2 treatment, thus highlighted with cyan. Fx, cellular function; CP, canonical pathway.

Figure 7

TGFβ signaling pathway overlaid with differentially expressed genes. Cellular localization of differentially expressed genes active in canonical TGFβ pathway was illustrated using the IPA Path Designer. Genes in Table 1 were uploaded to IPA as a dataset and analyzed using default settings. Acquired networks were overlaid with TGFβ canonical pathway. The overlapping genes are: TGFβ1, TGFβ2, INHBA (activins/inhibins), BMPR2, SMAD7, SMURF1, SMAD 1/5/8, SMAD2/3, VDR, and PAI-1. Green color represents downregulated expression, whereas red color represents upregulated expression of genes. Color intensity reflects the level of expression. Types of molecules are annotated in the legend in the box.

Table 3

List of differentially expressed genes targeting TGFβ signaling molecules.

Gene name	TGFβ1 FC	TGFβ2 FC	Target TGFβ molecule
BMPR2	1.42	1.38	SMURF1
CTGF	2.00	1.95	SMAD7
			SMAD3
			TGFβ1
			TGFβ2
ITGB3	5.19	4.42	PAI1
SERPINE1 (PAI1)	1.83	1.96	SMAD7
			PAI1
			SMAD3
			TGFBR1
			TGFβ1
			TGFβ2
SMAD1	1.48	NS	SMURF1
			BMPR2
			TGFBR1
SMAD3	−3.24	−2.51	SMAD7
			SMAD3
			BMPR2
			TGFBR1
			TGFβ1
			TGFβ2
SMAD7	4.17	4.04	SMAD7
			SMURF1
			INHBA
			TGFβ1
			TGFβ2
TGFB1	1.91	1.59	PAI1
			SMAD3
TGFBR1	1.49	NS	SMAD7
			SMURF1
			TGFBR1

NS, non-significant expression detected by dChip; FC, fold-change from microarray analysis.

Table A4

The full list of genes regulated by differentially expressed TGFβ signaling molecules and their level of altered expression.

Molecule	Gene	TGFβ1	TGFβ2
SMAD7	Alpha1	NS	NS
	CDKN1A	NS	NS
	CDKN2B	NS	NS
	COL1A1	NS	NS
	COL3A1	NS	NS
	Collagen	NS	NS
	CTGF	2.00	1.95
	MAPK	NS	NS
	P38	NS	NS
	SERPINE1	1.83	1.96
	SMAD2	NS	NS
	Smad2/3	NS	NS
	SMAD3	−3.24	−2.51
	SMAD7	4.17	4.04
	TGFBR1	1.49	NS
	TGFBR2	NS	NS
SMURF1	BMPR2	1.42	1.38
	MAP3K2	NS	NS
	RHOA	NS	NS
	RUNX2	NS	NS
	SMAD1	1.48	NS
	SMAD1/5/8	NS	NS
	SMAD4	NS	NS
	SMAD5	NS	NS
	SMAD7	4.17	4.04
	TGFBR1	1.49	NS
	TRAF1	NS	NS
	TRAF4	NS	NS
BMPR2	ID1	NS	NS
	ID2	NS	NS
	MAP3K7	NS	NS
	MAPK	NS	NS
	P38	NS	NS
	PKA	NS	NS
	Ras	NS	NS
	SHC1	NS	NS
	SMAD1	1.48	NS
	Smad1/5/8	NS	NS
	SMAD3	−3.24	−2.51
	SMAD5	NS	NS
TGFBR1	CDKN1A	NS	NS
	Erk1/2	NS	NS
	MAP3K7	NS	NS
	Ras	NS	NS
	SERPINE1	1.83	1.96
	SMAD1	1.48	NS
	SMAD2	NS	NS
	SMAD3	−3.24	−2.51
	SMAD4	NS	NS
	Tab1–Tak1	NS	NS
	TGFBR1	1.49	NS
	Gene	TGFβ1	TGFβ2
SERPINE1 (PAI1)	F2	NS	NS
	FLT1	NS	NS
	FN1	NS	NS
	ITGAV	NS	NS
	ITGB3	5.19	4.42
	Laminin	NS	NS
	LRP1	NS	NS
	PLAT	NS	NS
	PLAU	NS	NS
	PLAUR	NS	NS
	PLG	NS	NS
	SERPINE1	1.83	1.96
	TGFB1	1.91	1.59
	VTN	NS	NS
VDR	BGLAP	NS	NS
	Ca2 +	NS	NS
	CYP24A1	NS	NS
	CYP27B1	NS	NS
	FGF23	NS	NS
	IL2	NS	NS
	PTH	NS	NS
	REN	NS	NS
	Renin	NS	NS
	RXRA	NS	NS
	S100G	NS	NS
	SLC34A1	NS	NS
TGFβ1	ACTA2	NS	NS
	CDH1	NS	NS
	CDKN1A	NS	NS
	CDKN2B	NS	NS
	COL1A1	NS	NS
	COL1A2	NS	NS
	CTGF	2.00	1.95
	FN1	NS	NS
	FOXP3	NS	NS
	SERPINE1	1.83	1.96
	SMAD2	NS	NS
	SMAD3	−3.24	−2.51
	SMAD4	NS	NS
	SMAD7	4.17	4.04
	Gene	TGFβ1	TGFβ2
SMAD3	Alp	NS	NS
	CDKN1A	NS	NS
	CDKN2B	NS	NS
	COL1A2	NS	NS
	CTGF	2.00	1.95
	FSHB	NS	NS
	IFNG	NS	NS
	IL2	NS	NS
	NOS2	NS	NS
	SERPINE1	1.83	1.96
	SMAD3	−3.24	−2.51
	TBX21	NS	NS
	TGFB1	1.91	1.59
INHBA	ACVR1B	NS	NS
	ACVR2B	NS	NS
	BAX	NS	NS
	CCND2	NS	NS
	CDKN1A	NS	NS
	CDKN2B	NS	NS
	CGA	NS	NS
	CYP11A1	NS	NS
	FSH	NS	NS
	FSHB	NS	NS
	FSHR	NS	NS
	FST	NS	NS
	GNRHR	NS	NS
	INSULIN	NS	NS
	SMAD7	NS	NS
TGFβ2	CDKN1A	NS	NS
	CDKN2B	NS	NS
	COL1A1	NS	NS
	COL1A2	NS	NS
	CTGF	2.00	1.95
	FN1	NS	NS
	FOS	NS	NS
	FOXP3	NS	NS
	SERPINE1	1.83	1.96
	SMAD2	NS	NS
	SMAD3	-3.24	-2.51
	SMAD4	NS	NS
	SMAD7	4.17	4.04
	TNF	NS	NS

NS, non-significant expression detected by dChip.

Discussion

Following vertical positioning of the palatal shelves along both sides of the tongue, each shelf elevates to a horizontal position above the tongue. In order to explain this reorientation and extension several hypotheses have been proposed, including: mesenchymal proliferation and organization (Jin et al., 2010), mesenchymal and muscular contraction (Innes, 1978), generation of hydrostatic forces by hyaluronan (Goudy et al., 2010), and epithelial reorganization (Babiarz et al., 1979). Taken as a whole, cellular changes within the palatal mesenchyme, which composes the majority of the palatal shelves, governs the optimum growth and elevation of the palatal shelves prior to adherence and fusion (Ito et al., 2003; Iwata et al., 2010, 2011). HEPM cells, which were extracted from a single human abortus at the time of palatal shelf elevation (Yoneda and Pratt, 1981), have previously been used in several studies, including: determining the mechanisms of cleft palate development induced by various teratogens (Dhulipala et al., 2004, 2011), the role of various growth factors active in the regulation of palate development (Takechi et al., 2008), and the teratogenic potential of environmental agents (Watanabe et al., 1990). Therefore, HEPM cells serve as a physiologically relevant model to study differential expression profiling of genes which function in palatal growth, elevation, and extension; which are all indispensable for mammalian palatogenesis.

We analyzed the raw data of scanned GeneChip images using dChip, which normalizes and models the signal values of transcripts in order to accurately detect differential expression with FC values. We set the p-value <0.05 as a cut-off criterion of statistical significance for all raw data analysis. Although the raw signal values of hybridized transcripts, i.e., prior to statistical filtering, exhibited differential expression of more than 12,000 genes in all TGFβ-treated samples (data not shown), TGFβ3 group could not pass the more stringent statistical significance threshold. Therefore, only 8 genes were detected to be differentially expressed in response to TGFβ3 treatment in HEPM cells (Table 1), while expression of more than 500 genes was altered in TGFβ1 and TGFβ2 group. Based on these criteria of microarray data analysis, we concluded that TGFβ3 treatment did not result in significant changes in expression levels of transcriptome of HEPM cells.

There are a number of methods to compare the differential expression of genes in different treatment group using multiple microarray experiments. For several years, Venn diagrams have been shown to be the most effective graphical analysis tool to demonstrate similarities, differences and distinctions within gene sets (Kestler et al., 2005; Pirooznia et al., 2007). In our study, we used Venn diagrams to illustrate the number of altered genes within the HEPM cells in response to TGFβ isoform treatments as overlapping and individual gene numbers (Figure 1). The number of genes differentially expressed by TGFβ2 is much higher than TGFβ3, but less than TGFβ1. This, combined with the observation that TGFβ2 is expressed primarily in the mesenchyme region (Fitzpatrick et al., 1990), indicates its important role within the palatal mesenchyme. Since TGFβ1 is expressed in both the palatal mesenchyme and epithelium, it has been proposed to facilitate cross-talk between the cell types during palate development (Fitzpatrick et al., 1990). Therefore, it is expected that TGFβ1 would affect more genes compared to TGFβ2 and TGFβ3. Surprisingly, the number of genes significantly altered by TGFβ3 was minimal, which may suggest that TGFβ3 is less involved in the regulation of cellular functions of palatal mesenchyme. However, this does not rule out its major role within the MEE during palatogenesis. Meanwhile, the high number of differentially up- and downregulated genes within the intersection regions, particularly between TGFβ1 and TGFβ2, suggests that TGFβ isoforms may simultaneously coordinate the development of the palatal mesenchyme. We identified a strong correlation of gene expression patterns between microarray and qPCR data in terms of up- and downregulation of significantly differentially expressed genes, but not relative level of expression. Expression profiles of all 27 genes tested were found to be in agreement between the two methods (Figure 3), further confirming the reliability of the hybridization-based gene expression assay.

Visualization of the vast amount of microarray data is an essential step for the accurate interpretation of the acquired data. Unsupervised hierarchical clustering presented as a heatmap and dendrogram is a common approach for detecting coexpressed groups of genes (Eisen et al., 1998). The heatmaps of differentially expressed genes according to cellular biofunctions allowed us to detect progressively and distinctly up- or downregulated genes. The heatmaps also showed us that most of the genes altered significantly in response to TGFβ1 and TGFβ2 followed a very similar pattern, unlike TGFβ3, which effected gene expression to a lower extent in the HEPM cells (Figure 4).

The most significantly upregulated transcripts

There were 27 genes grouped as the most significantly altered, in which 13 of them were upregulated and 14 of them were downregulated, in response to TGFβ stimulation (Table 2). Although, not all of them presented a direct/indirect relation with TGFβ molecule according to IPA (Figure 2), here, we provide their potential interaction with TGFβ and/or cleft palate based on most recent literature review. Further detailed discussion of all 27 genes can be found at Table A3 in Appendix.

Table A3

Details of the most significantly altered genes within the palatal mesenchyme.

	Gene	Entrez Gene Name	Function in the cell
Most Upregulated	AMIGO2	Adhesion molecule with Ig-like domain 2	Cell–cell adhesion, apoptosis, adhesion, survival
	NOX4	NADPH oxidase 4	Proliferation, expression, apoptosis, morphogenesis, aging, growth, cell death, NADPH oxidase activity
	TNFAIP6	Tumor necrosis factor, alpha-induced factor 6	Signaling, adhesion, expansion, assembly
	DACT1	Dapper homolog 1	Multicellular organismal development; Wnt receptor signaling pathway
	ANKRD1	Ankyrin repeat domain 1	Apoptosis, colony formation
	PMEPA1	Prostate transmembrane protein, androgen induced 1	Apoptosis, growth; androgen receptor signaling pathway
	MSC	Musculin	Myogenesis; skeletal muscle development; palate development; regulation of transcription
	EDN1	Endothelin 1	Proliferation, migration, growth, apoptosis, gene expression, stimulation, pH, activation, invasion
	VDR	Vitamin D (1,25-dihydroxyvitamin D3) receptor	Apoptosis, development, proliferation, homeostasis, differentiation, osteoclastogenesis, transcription
	KCNH1	Potassium voltage-gated channel member 1	Growth, proliferation, fusion, hyperexcitation
	HAPLN1	Hyaluronan and proteoglycan link protein 1	Cell adhesion; hyaluronic acid binding
	NEDD9	Neural precursor cell expressed, developmentally downregulated 9	Migration, invasion, proliferation, apoptosis, quantity, cell division, cell spreading, survival, activation
	ITGB3	Integrin, beta 3	Adhesion, migration, aggregation, cell spreading, proliferation, activation, apoptosis, invasion, angiogenesis; negative regulation of cell death
Most Downregulated	GDF5	Growth differentiation factor 5	Differentiation, apoptosis, signaling, adhesion, proliferation, migration; cell–cell signaling; limb morphogenesis;TGFβR signaling pathway
	SECTM1	Secreted and transmembrane 1	Immune response; mesoderm development; positive regulation of IKB kinase/NFKB cascade; signal transduction
	NPTX1	Neuronal pentraxin I	Synaptic transmission, cell death, outgrowth, apoptosis, cell–cell contact, loss, chemotaxis; central nervous system development
	SMAD3	SMAD family member 3	Growth, proliferation, apoptosis, differentiation, chemotaxis; cell–cell junction organization; cell cycle arrest; developmental growth
	OAS1	2′–5′-oligoadenylate synthetase 1	Apoptosis; cellular response to interferon-alpha; cytokine-mediated signaling pathway; defense response to virus
	TBC1D8	TBC1 domain family, member 8	Proliferation, formation; blood circulation; cell proliferation; regulation of Rab GTPase activity
	TRERF1	Transcriptional regulating factor 1	Replication; homeostatic process; multicellular organismal development; regulation of transcription
	VSIG1	V-set and immunoglobulin domain containing 1	Member of the junctional adhesion molecule (JAM) family
	TGFBR3	Transforming growth factor, beta receptor III	Growth, apoptosis, proliferation, invasiveness, motility, movement, differentiation; epithelial to mesenchymal transition; TGFβ receptor signaling
	CA12	Carbonic anhydrase XII	Growth; one-carbon metabolic process; type I membrane protein
	ARHGAP28	Rho GTPase activating protein 28	Positive regulation of small GTPase activity; signal transduction
	HGF	Hepatocyte growth factor	Migration, scattering, proliferation, apoptosis, invasion, growth, expression, motility, branching, morphogenesis
	NR4A2	Nuclear receptor subfamily 4, group A, member 2	Differentiation, maturation, degeneration, quantity, migration, growth, survival; negative regulation of apoptosis
	IFIT2	Interferon-induced protein with tetratricopeptide repeats 2	Cellular response to interferon-alpha; cytokine-mediated signaling pathway; negative regulation of protein binding; type I interferon-mediated signaling pathway
Cleft Palate-related	EDN1	Endothelin 1	Proliferation, migration, growth, apoptosis, gene expression, stimulation, pH, activation, invasion
	INHBA	Inhibin, beta A	Apoptosis, growth, differentiation, cell cycle progression, proliferation, quantity, induction, cell death
	LHX8	LIM homeobox 8	Development, maturation, differentiation, quantity; female gonad development; forebrain neuron development
	PDGFC	Platelet derived growth factor C	Proliferation, migration, formation, mitogenesis, growth, chemotaxis, motility; activation of transmembrane receptor protein tyrosine kinase activity
	PIGA	Phosphatidylinositol glycan anchor biosynthesis, class A	Morphogenesis, colony formation, hemolysis, quantity; biosynthetic process; cellular protein metabolic process
	RUNX1	Runt-related transcription factor 1	Differentiation, apoptosis, growth, proliferation, transcription, activation, survival; skeletal system development; DNA-dependent transcription
	SMAD3	SMAD family member 3	Growth, proliferation, apoptosis, differentiation, chemotaxis; cell–cell junction organization; cell cycle arrest; developmental growth
	SNAI1	Snail homolog 1	Apoptosis, epithelial–mesenchymal transition, migration, cell cycle progression, binding, transcription, dissemination; osteoblast differentiation; palate development
	TGFβ1	Transforming growth factor, beta 1	Multifunctional peptide that regulates proliferation, differentiation, adhesion, migration
	TGFβ2	Transforming growth factor, beta 2	Proliferation, apoptosis, growth, differentiation, cell cycle progression, cell death, migration, cytostasis, development; angiogenesis; axon guidance; blood coagulation
	TGFβR1	Transforming growth factor, beta receptor 1	Apoptosis, differentiation, proliferation, growth, migration, motility, cell cycle progression, phosphorylation; signal transduction; skeletal system development

In our study, AMIGO2 was the most significantly upregulated gene in TGFβ-treated HEPM cells. In a recent study by Gimelli et al. (2011), it has been shown that a patient exhibiting cleft palate, mental retardation, and high myopia carries a de novo 12q13.11 microdeletion, which harbors 16 genes including AMIGO2. This genetic analysis shows that AMIGO2 might have a role in one of the symptoms mentioned above, including cleft palate.

Several studies showed that TGFβ1 regulates expression of NOX4 in airway smooth muscle cells (Sturrock et al., 2006), cardiac fibroblasts (Cucoranu et al., 2005), renal epithelial cells (Rhyu et al., 2005), and keratinocytes (Davies et al., 2005). Upregulated expression of Nox4 in palatal mesenchyme may assign a crucial role to Nox4 during palate development under TGFβ regulation.

In our study, TGFβ treatment induced high-level expression of TNFAIP6 and HAPLN1, both hyaluronan-related genes, in the palatal mesenchyme, which may suggest that TGFβ mediates hyaluronan production required for palatal shelf reorientation and extension.

Expression levels of two Wnt-mediated genes, DAPPER1 (DACT1) and ANKRD1 (CARP) genes were upregulated in response to TGFβ treatment. During embryonic development, TGFβ and Wnt signaling exhibit synergistic effects (Letamendia et al., 2001; Lei et al., 2004; Deng et al., 2008; Owens et al., 2008; Micalizzi et al., 2010), and both DAPPER1 and ANKRD1 are involved in both of these pathways (Waxman et al., 2004; Labbe et al., 2007; Su et al., 2007). Based on our detection of their upregulated expression, we suggest that these genes may serve as mediators of cross-talk between TGFβ and Wnt signaling during palate development.

PMEPA1 is a TGFβ-induced transmembrane protein (Brunschwig et al., 2003) and was upregulated following TGFβ stimulation. Although PMEPA1’s distinct role in development has not been studied yet, its upregulated expression and association with TGFβ make it an important candidate for the regulation of palatal mesenchyme development.

EDN1, which is a candidate for non-syndromic orofacial cleft disease based on evidence from linkage analysis and mouse model knockouts (Kurihara et al., 1994; Sato et al., 2008), was also upregulated. Thus, we suggest that TGFβ isoforms may regulate the expression of EDN1 and/or its transcriptional activators, which requires further elucidation.

The vitamin D receptor (VDR), which encodes the nuclear hormone receptor for vitamin D3, was also highly stimulated by TGFβ treatment (Table 2). The TGFβ and vitamin D signaling pathways are involved in synergistic and antagonistic cross-talk, in a tissue-specific manner, facilitated through the binding of SMAD3 and VDR to their associated binding sites (Yanagisawa et al., 1999; Subramaniam et al., 2001). Given its role as a mediator in normal physiology and its high-level expression in response to TGFβ treatment, we propose that VDR mediates the communication between TGFβ and Vitamin D during palatogenesis. Further verification of this intermediation may also shed light on the importance of Vitamin D as an environmental etiological factor of cleft palate.

It has been proposed that palatal shelf elevation can be driven by the growth (and movement) of the facial skeletal muscles inducing mechanotransduction pathways, which in turn regulate cytoskeleton remodeling, cell proliferation, and tissue differentiation (Innes, 1978; Farge, 2011). MSC (musculin), also known as MyoR, is a basic helix-loop-helix (bHLH)-type transcriptional regulator of myogenesis; and we detected its expression as upregulated with TGFβ treatment. Considering its capability of binding to the E-box element and muscular movement during palatogenesis, we postulate that MSC and TGFβ isoforms can modulate each other’s expression in a feedback loop during palate formation.

In a study conducted by Vogel et al. (2010), it is documented that TGFβ1 regulates expression of NEDD9, a scaffolding protein, to promote the differentiation of hippocampal and cortical cells into neurons. Similar to TGFβ isoforms, NEDD9 is associated with diverse biological processes including cell attachment, migration, and invasion, as well as apoptosis and cell cycle regulation. Considering its high-level expression observed in palatal mesenchyme, it is conceivable that NEDD9 may mediate crucial processes during palatogenesis in a close relation with TGFβ.

The most significantly downregulated transcripts

GDF5 is a member of the TGFβ superfamily and is involved in cell adhesion and differentiation of mesenchymal cells to chondrocytes (Takahara et al., 2004). Considering its reduced expression in the embryonic palatal mesenchyme, we hypothesize that TGFβ signaling inhibits GDF5 expression to block the mesenchymal differentiation into ossification centers, in favor of mesenchymal proliferation during palatal growth.

Based on the network analysis, which is developed using the IPA (data not shown), SECTM1 has an indirect relationship with SMAD3 and HGF, both of which were downregulated upon TGFβ treatment. Although, to our knowledge, there is no study available focusing on the role of SECTM1 during development, its downregulated expression in response to TGFβ may be associated with SMAD3’s decreased expression, which may suggest that SECTM1 is regulated under TGFβ signaling.

Human NPTX1 mRNA is exclusively localized to the nervous system during development (Omeis et al., 1996), and promotes apoptosis within the developing cortical neurons in culture (Abad et al., 2006). Although its function during craniofacial development has not been described yet, its significant downregulation in TGFβ-treated HEPM can be explained by requirement for inhibition of cell death during the morphogenesis of palatal mesenchyme.

Another apoptosis-related gene, OAS1 (Hale et al., 2008), was also downregulated with TGFβ treatment. Since the palatal mesenchyme is constantly proliferating and growing during the initial stages of palatogenesis, reduced expression of OAS1, resulting in inhibition of apoptosis, is a required step for the continuous development of the palate.

SMAD3 is one of the intracellular mediators that transduce signals from TGFβ receptors. It has been demonstrated that SMAD2/3 is involved in growth inhibition of mouse embryonic palate mesenchymal cells (MEMP) induced by all-trans retinoic acid (Wang et al., 2009). Our study showing downregulated expression of SMAD3 in response to TGFβ may suggest that SMAD3 is required to be inhibited during palatal development to maintain consistent growth of palatal shelves until attachment.

The VSIG1 gene encodes a member of the junctional adhesion molecule (JAM) family. Since the mesenchymal cells of the palate do not express adherens junction proteins (Nawshad et al., 2004), reduced expression of VSIG1 by TGFβ suggests that regulation of VSIG1 is required for integrity of mesenchymal cells of developing palate.

In epithelial cells, TGFBR3 inhibits TGFβ signaling by preventing the formation of the TGFBR1–TGFBR2 complex, thereby functioning as a potent TGFβ neutralizing agent (Vilchis-Landeros et al., 2001; Eickelberg et al., 2002). For the maintenance of proper communication of TGFβ signaling within the palatal mesenchyme, it is essential that expression levels of TGFBR3 are minimized; thus both TGFβ1 and TGFβ2 treatment facilitated downregulation of TGFBR3.

HGF has been demonstrated to cell growth (Boccaccio et al., 1998), cell motility (Hajjar and Nachman, 1996), and morphogenesis (Sunil et al., 2002) by activating a tyrosine kinase signaling cascade (Porter and Vaillancourt, 1998). In several types of cells, it has been shown that TGFβ, even at low concentrations, significantly reduces HGF production [87; 88]. Based on downregulated expression of HGF, we hypothesize that there is a negative feedback between TGFβ and HGF for the proper development of the palate.

Conclusion

In conclusion, we have determined that more than 2% (n = 566) of the human genes were differentially expressed in palatal mesenchymal cells in response to treatment with TGFβ isoforms. According to GO annotation clustering, 234 of these altered genes (41%) were associated with cellular biofunctions. Nevertheless, only 11 of these transcripts were implicated in craniofacial deformities, which suggest that other etiological factors are correspondingly involved in the development of palatal constituents. Our overall results suggest that both TGFβ1 and TGFβ2 isoforms regulate the expression of genes that govern cellular biofunctions of the palatal mesenchyme; and this regulation is crucial for normal palatogenesis. Further elucidation of the significantly up- and downregulated genes at the protein level will enhance our understanding of the mechanisms controlling palate development; thereby pave the way for prevention of cleft palate during development.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.