Transcriptional analysis of cleft palate in TGFβ3 mutant mice

Cleft palate (CP) is one of the most common craniofacial birth defects, impacting about 1 in 800 births in the USA. Tgf-β3 plays a critical role in regulating murine palate development, and Tgf-β3 null mutants develop cleft palate with 100% penetrance. In this study, we compared global palatal transcriptomes of wild type (WT) and Tgf-β3 −/− homozygous (HM) mouse embryos at the crucial palatogenesis stages of E14.5, and E16.5, using RNA-seq data. We found 1,809 and 2,127 differentially expressed genes at E16.5 vs. E14.5 in the WT and HM groups, respectively (adjusted p < 0.05; |fold change|> 2.0). We focused on the genes that were uniquely up/downregulated in WT or HM at E16.5 vs. E14.5 to identify genes associated with CP. Systems biology analysis relating to cell behaviors and function of WT and HM specific genes identified functional non-Smad pathways and preference of apoptosis to epithelial-mesenchymal transition. We identified 24 HM specific and 11 WT specific genes that are CP-related and/or involved in Tgf-β3 signaling. We validated the expression of 29 of the 35 genes using qRT-PCR and the trend of mRNA expression is similar to that of RNA-seq data . Our results enrich our understanding of genes associated with CP that are directly or indirectly regulated via TGF-β.


Results
Our previously acquired RNA-seq data 12 involving wild type and Tgf-β3 homozygous knockout mouse models on E14.5 and E16.5 represented four sample groups (denoted as "WT14.5, " "WT16.5, " "HM14.5, " and "HM16.5"), each represented by two biological replicates (denoted as "a, " or "b") for a total of eight samples. The average read count for the raw RNA-seq data was ~ 65.5 M paired-end reads (i.e., ~ 130.1 M total reads) per sample providing a high coverage of the transcriptome. The re-analysis of this data involved data trimming and filtering, normalization, expression quantification, clustering, differential analysis, comparative analysis, systems biology analysis, and qRT-PCR validation.
RNA-seq data processing and clustering. After trimming and filtering, the number of average total reads per sample came down to ~ 128.1 M (Supplementary Fig. 1A). The average read length was 101 bp in the raw data, which decreased to 95.86 bp following trimming and filtering ( Supplementary Fig. 1B). On the other hand, the average Phred read quality score increased to 36.79 from 35.84 ( Supplementary Fig. 1C); and the percentage of high-quality bases (bases with a quality score > 20) per sample increased to 99.20% from 96.28% ( Supplementary Fig. 1D) following trimming and filtering. Therefore, both the total number of reads and the average read length parameters showed small changes in quantity after trimming and filtering; but this resulted in significant data quality improvement.
RNA-seq analysis generated expression data for 103,215 transcripts. Transcripts that showed a TPM value of less than 1 in both samples in all of the four groups were eliminated from downstream analysis leaving 52,475 transcripts. Biological replicates showed high degree of correlation (r > 0.98). In Supplementary Fig. 2, we show the hierarchical clustering of the samples using all 52,475 transcripts. This global unsupervised grouping reveals that the samples were separated clearly by time as there are two main clades where one clade only consists of E14.5 samples and the other clade only consists of E16.5 samples. The effect of genotype on the transcriptional profiling was subtle because in the clustering tree we do not see a grouping based on the genotype neither globally, nor within a time point. Hence, there is a need for supervised analysis methods to identify the differences in gene expression due to genotypic variance. Furthermore, the height of the branching points in Supplementary  Fig. 2 both for the E14.5 and E16.5 sample groups implied that the similarity between the WT and HM groups was higher at E14.5 than it was at E16.5. This is because of the shorter branch lengths among samples in the E14.5 clade compared to the branch lengths among samples in the E16.5 clade. Therefore, the samples at E16.5 showed a more divergent transcriptional profile and thus the effects of Tgf-β3 knockout were more pronounced at E16.5. Differential expression analysis. Differential expression analysis also showed the stark temporal difference in gene expression between samples (Table 1). There were 4,115 and 5,304 significantly differentially expressed genes (SDEG) between the E16.5 and E14.5 time points for the WT and HM samples, respectively indicating that the temporal transcriptional change in HM samples was more widespread than in normal controls. Furthermore, in both genotypes, the number of upregulated genes with time (2,421 for WT and 3,153 for HM) was more than the number of downregulated genes ( www.nature.com/scientificreports/ We observed a very subtle difference between the two genotypes at a given time point. At days E14.5 and E16.5, there were only 13 and 38 SDEGs between the WT and HM samples, respectively. In concordance with the unsupervised hierarchical clustering results, we saw a greater difference between the genotypes at E16.5. This was further strengthened by observing 13 SDEGs between the two genotypes with an |FC|> 2.0 at E16.5, while there were no SDEGs between the WT and HM samples with an |FC|> 2.0 at E14.5. WT specific and HM specific gene identification. Since a direct comparison between the HM and WT groups yielded a very subtle difference, we defined the effect due to Tgf-β3 knockout by comparing the temporal SDEGs in the two groups. Following our adjusted p < 0.05 and |FC|> 2.0 cutoffs, we compared the 1,675 and 1,936 SDEGs that were upregulated in E16.5 vs. E14.5 in the WT and HM groups, respectively. We identified 429 genes that were upregulated in the WT group, but not in the HM group; similarly, we identified 690 genes that were upregulated in the HM group, but not in the WT group. Conversely, we compared the 134 SDEGs downregulated in E16.5 vs. E14.5 in the WT group with the 191 SDEGs downregulated in E16.5 vs. E14.5 in the HM group. We identified 72 genes that were downregulated in the WT group, but not in the HM group and similarly we identified 129 genes that were downregulated in the HM group, but not in the WT group. These results are summarized in Fig. 2A. We called the 501 (429 + 72) SDEGs uniquely up/downregulated in the WT group "WT specific"; and similarly, we called the 819 (690 + 129) SDEGs uniquely up/downregulated in the HM group "HM specific. " In Fig. 2B, we show the hierarchical clustering of the 1,320 (501 + 819) WT and HM specific SDEGs across the four sample groups. This visual representation of expression values (heatmap) clearly shows how WT specific genes are significantly up/downregulated between WT E14.5 and WT E16.5 samples, whereas the expression of these genes remains "not significantly altered" between the HM E14.5 and HM E16.5 samples. A similar (but complementary) observation holds for the HM-specific genes. Figure 2B also shows the sample similarity based on genes specific to the genotypes, which clustered the biological replicates together and grouped the samples in E14.5 under the same clade. This implies that the temporal difference in the dataset was more dominant than the genotypic difference as we do not see samples from the same genotype in the two different time points group together. The HM samples at E16.5 stood out as a distinct group, which is reasonable as the phenotypic difference between the genotypes was not visible at E14.5 but emerged clearly at E16.5, making the HM samples at E16.5 a distinctly separate phenotypic group among the four. www.nature.com/scientificreports/ Functional analysis. In order to understand the genes and corresponding signaling network and functional mechanisms that would explain the differences observed between HM and WT samples across E16.5 and E14.5 time points, we highlighted the GO functional categories and KEGG pathways that were statistically significantly enriched in the WT and HM specific gene lists (Fig. 2C). A complete list of SDEGs, enriched GO categories and KEGG pathways can be found in the Supplementary Data. In Table 2, we list 24 HM-specific and 11 WT-specific genes based on their functional relevance by any one or more of the following criteria: (1) falling under the TGF-β signaling pathway, (2) contributing to palatogenesis or cleft palate in mice, (3) being involved in palatal cellular processes, such as EMT, apoptosis, differentiation, proliferation that are functional in palatal and craniofacial morphogenesis, and (4) being directly or indirectly regulated by TGF-β signaling.
Systems biology analysis. Using IPA, we analyzed the signaling networks and crosstalks amongst the downstream molecules that could be uniquely functional in WT and HM palates. Our network analysis results suggest that the TGF-β, ERK/MAPK, p38MAPK and PI3K/AKT pathways directly regulate p38MAPK in WT ( Fig. 3A and Supplementary Fig. 3), which is not the case in HM ( Fig. 3B and Supplementary Fig. 7). The canonical pathway analysis using IPA further strengthened this observation when we specifically considered the TGF-β signaling pathway separately for WT specific genes (Fig. 3A) and for HM specific genes (Fig. 3B). Furthermore, we showed that Goosecoid (Gsc) is downregulated by TGF-β signaling in WT (Fig. 3A), whereas, in the HM, TGF-β signaling downregulates transcription factor Tlx2 (Fig. 3B). Finally, our results suggest that, in both in WT and HM, TGF-β signaling regulates expression of neither Smads nor of transcripts encoding proteins that act upstream or downstream of Smads (Fig. 3A,B; Supplementary Figs. 3 and 7). We identified WT specific genes that are also TGF-β downstream target genes. When we overlaid the Epithelial Mesenchymal Transition (EMT) functional category on these genes, we observed that only two genes were selected ( Supplementary Fig. 4a). When we similarly overlaid the "cell death" functional category on the same gene list, we observed nine selected genes ( Supplementary Fig. 5a). These results potentially demonstrate that in WT, TGF-β signals to regulate "cell morphology" and "differentiation" that facilitate apoptosis over EMT as the number of TGF-β downstream target genes are predominantly pro-apoptotic genes but not EMT. We further analyzed the mechanistic networks generated by IPA that involve TGF-β signaling. These networks show the interaction among upstream regulators that best explain the changes observed in the SDEG list. We performed this analysis for WT specific genes and found that these upstream regulators were predominantly pro apoptotic, (12 out of 14, Supplementary Fig. 5b) not EMT (none out of 14, Supplementary Fig. 4b) genes.

Discussion
All isoforms of TGF-β are involved in both palatal epithelial and mesenchymal cell proliferation, differentiation, transformation, and apoptosis 16 . TGF-β3 plays a distinctive role in periderm removal, palatal fusion and seam disintegration chronologically that involves crosstalk and reciprocal signaling between the two cell populations by epithelial-mesenchymal interaction (EMI) 17,18 , which cannot be achieved by other isoforms 10,11 . Therefore, in this study, we analyzed Tgf-β3 HM whole palates to accurately explain the global TGF-β3 signaling during palate development. Usage of both epithelial and mesenchymal cells is a standard practice as previously done 19,20 . In one of the most elegant studies of the palatal transcriptome analysis by NIH Facebase consortium, the authors also used regions of palate with a combination of both palatal epithelia and mesenchyme 13 .
Our differential expression analysis showed a very subtle difference between the WT and HM palates-there were only 13 and 38 SDEGs between the WT and HM samples, respectively, at days E14.5 and E16.5 (see Supplementary Data). In agreement with the unsupervised hierarchical clustering data, we noticed a greater difference between the genotypes at E16.5. This difference was not only in numbers but also in the degree of differential expression. When the |FC|> 2.0 cutoff was employed, there were no SDEGs at E14.5 and only 13 SDEGs at E16.5 between the HM and WT groups (Table 1). This observation may be attributable to the obvious phenotypical difference, a palatal cleft at E16.5 and no difference at E14.5. Interestingly, 9 out of the 13 SDEGs in HM vs. WT at E16.5 were epithelial-specific genes with a majority of them being involved in palate deformities and cleft palate (Table 3). Henceforth, the differences between these two genotypes mainly impacted palatal epithelial cell-specific function, which is in accordance with the expression of TGF-β3 in WT or lack of it in HM.
The subtle difference in gene expression at a given time point is likely to trigger larger downstream changes through direct and/or indirect regulatory signaling that results in a phenotypic difference. Such a large effect would be observed through variations in the temporal gene expression patterns for the two genotypes, which also accounts for any subtle cell population differences as direct comparisons between the genotypes are not involved. Therefore, we focused on genes that show temporal changes in WT (from E14.5 to E16.5) due to the presence of Tgf-β3, which are not identified as changed in HM (from E14.5 to E16.5) due to the absence of Tgf-β3 Scientific RepoRtS | (2020) 10:14940 | https://doi.org/10.1038/s41598-020-71636-0 www.nature.com/scientificreports/ and vice versa. In order to identify genes that show temporal difference uniquely in the WT or HM groups, we identified WT specific and HM specific gene lists ( Fig. 2A). The functional analysis of these gene lists (Fig. 2C, Supplementary Data) render GO categories and KEGG pathways that are statistically significantly enriched in the WT or HM, such as the "apoptotic process, " "cell adhesion molecules (CAMs), " and "focal adhesion. " These lists show specific signaling pathways and TGF-β downstream molecules modulating both the epithelial and the mesenchymal cellular functions that are unique to WT and HM, resulting in immaculate palatogenesis in WT but cleft in HM.
In the Supplementary Data, we also list the 1,308 genes commonly up/down regulated between E16.5 and E14.5 both in WT and HM groups. The rationale behind the table is to describe the common genes that are responsible for phases that are common in both WT and HM, such as palatal cell proliferation, differentiation, transformation and apoptosis. These cellular changes are ongoing in both genotypes, hence, the genes listed in this Supplementary Data represent a common palatal cellular behavior and function that may not necessarily imply genes that can trigger cleft palate. Of note among commonly upregulated genes are Sntn, which participates in cell growth/maintenance and is involved in cell communication process in the nasopharynx 21 , and Tmem212, which codes for a transmembrane protein known to interact with transcription factor Tcf12, an important paralog of Tcf4 22 and Forkhead box k2 (Frk2). Frk2 is an important regulator palatal cell proliferation via activation of the phosphoinositide 3-kinase/AKT pathway 23 . Other genes to names are, Ppp1r32, which is a substrate for cGMP-dependent protein kinase and is involved in central nervous system development and intracellular signal transduction 24 , is upregulated in both WT and HM palates. It implements protein serine/threonine phosphatase inhibitor activity and inhibits phosphatase activities of protein phosphatase 1 (Pp1) and protein phosphatase 2A (Pp2A) complexes. Pp1 and Pp2A were reported by Weston et al. 25 to account for virtually all detectable serine/ threonine protein phosphatase activity during the development of embryonic palate. We also observed common upregulation of Tspan1 and 33 in both WT and HM going from E14.5 to E16.5. Tspan1have been shown to play crucial roles in biologic processes including cell adhesion, proliferation, differentiation, and migration 26 . Our data suggest that the role of Tspans may be limited to cellular proliferation and differentiation via Smad pathways. Another commonly upregulated transcript is Tnfrsf11a (RANK), which is a key regulator of bone homeostasis 27 . As early as E14.5, mesenchymal condensations undergo chondrogenesis initially and ultimately membranous ossification that give rise to the hard palate 28 . We expect the palatal cells express RANK to regulate osteoclast function in palatal bone formation in both genotypes. In terms of epithelial markers, Krt13, which has been shown to be expressed in the suprabasal layer of stratified palatal keratinized epithelia 29 , is also among the commonly upregulated transcripts. Krt4, a type II cytokeratin, is specifically expressed in differentiated layers of all of oral mucosal epithelia along with family member Krt13 30 . Similar to Krt13, Krt4 is also upregulated in both WT and HM. Showing a similar expression pattern to Krt13 among commonly upregulated transcripts is Sprr3 (aka Loricrin), a marker for terminally differentiated keratinized and non-keratinized oral mucosa 31 . This is potentially due to the fact that palatal epithelia undergo significant stratification with terminally differentiated keratinized epithelia covering the oral side of the palatal lining in both WT and HM palates.
Using IPA, we identified several novel findings which may suggest, and potentially imply, new signaling network molecules as well as cellular functions. Our data suggests that TGF-β signaling may induce palatogenesis through regulating p38MAPK in WT palates (Fig. 3A,B, Supplementary Figs.3 and 7). These findings are in agreement with previous work [32][33][34] showing p38MAPK activation by TGF-β signaling during palatogenesis. It has been previously demonstrated that in palatogenesis TGF-β3 signals through the SMAD pathway [35][36][37][38] . Our results indicate no change in SMADs and its down-and up-stream genes at the mRNA levels both for WT and HM palates (Fig. 3A,B, Supplementary Figs. 3 and 7). This observation does not necessarily imply an inactivation of the SMAD pathway as transcriptional level activity does not always imply protein level functionality. On the other hand, these results may imply non-Smad pathways to be also at play both in the WT and HM while in HM, unlike WT, this signaling cascade does not include p38MAPK (Supplementary Figs. 3 and 7).
Uniquely in the WT palate, Goosecoid (Gsc), a homeobox-containing gene, is downregulated by TGF-β signaling in palatogenesis (Fig. 3A). Gsc mutant mice display defects in the pharyngeal muscles and the pharyngeal mucosa 39 . TGF-β signaling is known to promote EMT by regulating the Gsc gene during embryonic Spemann's organizer formation 40 as well as breast cancer metastasis 41 . Since GSC is known to be a homeobox transcription factor that promotes EMT, it is likely that apoptosis is favored for seam disintegration since an EMT gene, such  (Fig. 3B). The Tlx2 gene encodes transcription factors essential in the development of neural-crest-derived cells suggesting a physiological role in the transcription-factor cascade underlying the differentiation of neuronal lineages during embryogenesis 42 . We did not find any direct relationship between TGF-β and Tlx2 but based on the fact that Tlx2 is crucial for neural-crest-derived cell development, it is therefore likely that palatal epithelia (which is an ectoderm derived cell) has limited or no role for Tlx2 and therefore it is downregulated. Our upstream regulator analysis identified targets of Tgf-β in the WT specific gene list. Exploring the involvement of these genes in the EMT and cell death mechanisms revealed limited involvement (Cdh1 and Lef1) of the EMT pathway ( Supplementary Fig. 4a) whereas nine apoptotic genes (Ace, Cdh1, Dcn, Dlx2, Krt18, Pparg, Rasgrp, Sema7a, and Tnfrsf11b) were directly regulated by TGF-β (Supplementary Fig. 5a). We also identified mechanistic networks, which are interaction networks of upstream regulators that explain the changes observed in the WT (or HM) specific genes. When we explored the functional characteristics of these interconnected upstream regulators, none of the genes were involved in the EMT (Supplementary Fig. 4b); and 12 out of 14 genes were involved in "cell death" (Supplementary Fig. 5b).
We identified genes regulated by Tgf-β3 in WT specific genes and showed that among these Tgf-β3 targets, Cdh1 is upregulated and Lef1 is downregulated (Supplementary Fig. 6). Loss of Cdh1 is a key marker of EMT 43 ; and the expression of Cdh1 can be repressed by the transcription factor Lef1 in palatal EMT 22,44 . Additionally, anti-EMT cell-cell adhesion genes (Cdh1, Ocln, Tspan2) are upregulated, which may suggest a potential relationship between Tgf-β3 signaling and suppression of palatal EMT. These findings may indicate two possible outcomes: (a) since palatal seam disintegration is complete at E16.5, the EMT markers are no longer expressed, or (b) EMT may not be a mechanism of palatal seam disintegration. Our findings, based on mechanistic network and upstream regulator analysis (Supplementary Figs. 4a, 4b, 4a, 4b, and 4), suggest that the latter is more probable as the number of apoptotic regulatory network genes are significantly more than the number of EMT genes.

Conclusion
Identifying transcripts that play key roles in regulating palatal development in critical stages has been a powerful approach to understanding how Tgf-β3 controls normal palatogenesis and how the lack of signaling (and its downstream signaling partners) is associated in induction of cleft palate. This study identifies potential CP-related genes based on differential expression between genotypes and gestational ages. The data presented in this work provide a strengthened understanding of the complex genetic mechanism of Tgf-β3-regulated palatogenesis. In addition, we discussed the variations in gene expression in the absence of Tgf-β3 in HM implicated in cleft palate. Our results represent a comprehensive analysis of the gene profile in murine cleft palate due to the absence of Tgf-β3. Further elucidation of the significantly up/downregulated genes will enhance our understanding of the mechanisms controlling palate development, thereby paving the way for prevention of cleft palate during development.

Methods
Animal selection and breeding. Tgf-β3 heterozygous (+/−) C57BL/6J breeder mice were obtained from Tom Doetschman (BIO5 Institute, University of Arizona, AZ). The reproduction and genotyping of Tgf-β3 −/− mice was conducted as previously described 9 . Mice were accommodated and subject to procedures at the University of Nebraska Medical Center (UNMC) College of Dentistry Animal Facilities under the approval of the UNMC Institution Animal Care and Use Committee (IACUC # 06-064). Null mutant embryos were generated by intercrossing Tgf-β3 heterozygous male and female mice in a Mendelian fashion.
Genomic DNA purification and genotyping. Palatal tissues were dissected under the NIKON SMZ1000 stereo microscope system (NIKON, Tokyo, Japan) from embryos collected on embryonic day (E) 14.5 and E16.5 following the identification of vaginal plugs, which are considered to be E0.5. Palatal samples were stored in RNAlater Stabilization Reagent (QIAGEN, Hilgen, Germany) to preserve the gene expression profile and individually labeled and matched with the corresponding tail tissue used for genotyping as detailed in our previous study 12 . Extraction of RNA, construction of small RNA libraries, and RNA-Seq. Two biological replicates from each genotype and gestational stage were designed to ensure reproducibility and rule out the possibility of differences caused by technical procedures. Palatal shelves were harvested in pairs from eight fetuses out of four litters. Each sample consisted of two pairs of palatal shelves dissected from fetuses of the same genotype from the same litter. The total RNA was purified using Arcturus PicoPure RNA Isolation Kit (THERMOFISHER   RNA-seq analysis. RNA-seq data was obtained for four groups: two genotypes, Tgf-β3 −/− HM and Tgf-β3 +/+ WT samples, profiled at two time points, embryonic days E14.5 and E16.5. Each group was represented by two biological replicates, resulting in eight samples. Each sample was run in two lanes on the ILLUMINA HISEQ2000 next generation sequencer using the 2 × 101 bp paired-end mode. Due to the high correlation coefficient between them (r > 0.995) and in order to increase the coverage per biological sample and reduce the lane effect 45 , lane data for each sample were pooled at the read level. Raw reads were analyzed with FASTQC (V. 0.11.5) for quality control 46 . Overrepresented (e.g., adapter and similar technical) sequences remaining in the raw reads were assessed and subsequently removed using TRIM-MOMATIC (V. 0.36) in the palindrome mode based on default alignment detection and scoring parameters 47 . Trimmomatic was also used for low quality base filtering. Maximum information quality filtering was employed with a minimum average read quality threshold of 25. Following technical sequence and low-quality base removal, reads that were shorter than 36 bp were filtered out.
Transcript quantification was done based on the GRCm38.p5 reference genome using Salmon (v. 0.8.2) with default parameters 48 . Salmon uses sample-specific models, such as correction for GC-content bias, that improve the accuracy of transcription abundance estimates. We used transcripts per million (TPM) in Salmon's output as the normalized relative abundance measure employed in our downstream analysis. Differential gene expression analysis was done using DESeq2 (Love, Huber et al. 2014), which has been shown to perform well in experimental designs with few replicates 49 . The RNA-seq data used in this paper were deposited under NCBI's Gene Expression Omnibus (GEO) database (Accession No.: GSE109838).
Clustering of samples and/or genes was done using the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) method (also known as hierarchical clustering) with Pearson's correlation as the distance measure 50 . The expression data matrix was row-normalized prior to the application of average linkage clustering. The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 51 was used for functional analysis of the gene lists, interrogating the Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) Gene Ontology (GO) categories 52 , and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways 53 . Biologically relevant categories that are overrepresented in the gene set and, therefore, may be of further interest were assessed using the Expression Analysis Systematic Explorer (EASE) score in the DAVID tool. The EASE score is the upper bound of the distribution of jackknife iterative resampling of Fisher's exact probabilities with Bonferroni multiple testing correction. Categories containing low numbers of genes are underweighted so that the EASE score is more robust than the Fisher exact test. The EASE score is a significance level, with smaller EASE scores indicating increasing confidence in overrepresentation. We selected GO categories that have EASE scores of 0.05 or lower as significantly overrepresented.
We further analyzed the differentially expressed gene lists using the Ingenuity Pathway Analysis (IPA; QIA-GEN Inc., https ://www.qiage nbioi nform atics .com/produ cts/ingen uity-pathw ay-analy sis) software. IPA is based on the manual curation of scientific literature to identify pathways, networks, and functional categories that are significantly represented in the input gene list 54 . The computational analysis methods used in IPA are based on enrichment approaches 55,56 where pathways or functional groups in which the input gene lists are overrepresented are identified. By the same token, IPA identifies upstream regulators (e.g., transcription factors, microRNAs, kinases, compounds, or drugs) and generates interaction networks (based on known interactions identified in the literature) that best explain the transcriptional changes observed in the input gene list. www.nature.com/scientificreports/ Confirmation of differentially expressed genes with qRT-PCR. To verify differentially expressed genes in Tgf-β3 WT and HM samples, qRT-PCR was undertake as previously described 19,57 . Embryonic palates were extracted from at E14.5 and E16.5, and RNA extraction was conducted as previously described 19,57 using Arcturus PicoPure RNA Isolation Kit (THERMOFISHER SCIENTIFIC, San Francisco, CA) to reliably extract high-quality RNA from a few cells. RNA (500 ng) was converted to cDNA using Invitrogen Superscript IV VILO Master Mix (THERMOFISHER SCIENTIFIC, San Francisco, CA) that generated a significant cDNA yield at high temperatures in less time. An additional preamplification step was performed using TaqMan PreAmp Master Mix (THERMOFISHER SCIENTIFIC,

Data availability
The RNA-seq data used in this paper were deposited under NCBI's Gene Expression Omnibus (GEO) database (Accession No.: GSE109838).