Transcriptome analysis of Xenopus orofacial tissues deficient in retinoic acid receptor function

Background Development of the face and mouth is orchestrated by a large number of transcription factors, signaling pathways and epigenetic regulators. While we know many of these regulators, our understanding of how they interact with each other and implement changes in gene expression during orofacial development is still in its infancy. Therefore, this study focuses on uncovering potential cooperation between transcriptional regulators and one important signaling pathway, retinoic acid, during development of the midface. Results Transcriptome analyses was performed on facial tissues deficient for retinoic acid receptor function at two time points in development; early (35 hpf) just after the neural crest migrates and facial tissues are specified and later (60 hpf) when the mouth has formed and facial structures begin to differentiate. Functional and network analyses revealed that retinoic acid signaling could cooperate with novel epigenetic factors and calcium-NFAT signaling during early orofacial development. At the later stage, retinoic acid may work with WNT and BMP and regulate homeobox containing transcription factors. Finally, there is an overlap in genes dysregulated in Xenopus embryos with median clefts with human genes associated with similar orofacial defects. Conclusions This study uncovers novel signaling pathways required for orofacial development as well as pathways that could interact with retinoic acid signaling during the formation of the face. We show that frog faces are an important tool for studying orofacial development and birth defects. Electronic supplementary material The online version of this article (10.1186/s12864-018-5186-8) contains supplementary material, which is available to authorized users.


Background
The orofacial region, the area between the eyes from the bridge of the nose through the chin, is an integral center for recognition, communication and ingestion. Birth defects affecting this region such as cleft lip and palate are among the most common worldwide. A complete understanding of how the face forms in the embryo is critical to developing methods to prevent these types of birth defects. The face develops through precise coordination of growth and convergence of craniofacial prominences, migration of neural crest cells, and interactions between mesoderm derived tissues and the surrounding ectoderm (reviewed in [1]). We know many of the key players responsible for the orchestration of orofacial development such as WNTs, bone morphogenetic proteins (BMPs), fibroblast growth factors, sonic hedgehog proteins, and retinoic acid (reviewed in [2][3][4][5][6]). However, our understanding of how such signaling pathways interact with each other and implement changes in gene expression is rudimentary. Therefore, this study focuses on uncovering potential cooperation between transcriptional regulators and signaling pathways, in particular retinoic acid, during orofacial development.
Perturbation of retinoic acid (RA) signaling can lead to the development of a cleft lip/palate in humans (reviewed in [7,8]). Work in animal models has further demonstrated that retinoic acid signals are specifically required for orofacial development [6,7,[9][10][11][12][13]. In particular, Xenopus retinoic acid pathway components are expressed in the developing midface and embryos exposed to an retinoic acid receptor (RAR) antagonist during early orofacial development form a median orofacial cleft [14]. RA ligand binds to a heterodimer of two nuclear receptors most often consisting of RXRs and RARs [15]. These receptors bind to specific enhancer regions in the DNA called retinoic acid response elements. Upon RA binding to RAR/RXR, complexes of coactivators and epigenetic regulators are recruited. These then in turn modify the chromatin structure, allowing the transcriptional machinery to access the DNA and transcription can proceed. Without RA ligand, the receptors are bound by corepressors and repressive epigenetic regulators that stabilize the nucleosome structure so that the DNA is inaccessible to the transcriptional machinery (reviewed in [16,17]). This balance of RAR activation and repression is integral in regulating gene expression during embryonic development [18].
We now know that RA can modulate the expression of hundreds of genes during development and the expression of such genes can differ widely across developmental events (for examples compare [19][20][21]). Thus, to gain a more complete understanding of the role of RA during midface development we; 1) examined global gene expression changes in embryos where retinoic acid signals are perturbed and 2) specifically analyzed expression changes in the orofacial tissues during two different phases of its development. By doing so, this work provides a comprehensive picture of how RA is required during orofacial development, independent of its roles in earlier whole body development. Further, we have identified novel signaling and transcriptional regulators that may coordinate with RA during the development of the face. Finally, our work reveals that many of the genes altered in Xenopus embryos with a median cleft are also implicated in humans with similar orofacial defects. As a whole, this work furthers our understanding of RA signaling during orofacial development and exhibits that frog faces are an ideal tool for craniofacial research, especially to formulate a more comprehensive understanding of the complex network of signals and transcriptional regulators of this region.

Inhibition of retinoic acid signaling during two phases of orofacial development showed overlapping and distinct phenotypes
To better understand the evolving role of retinoic acid during orofacial development, we perturbed RAR function over two distinct phases. Treatment 1 consisted of RAR antagonist administration during the early phase of facial development, from stage 24-30, (26-35 hpf ). At this time the neural crest is migrating and facial prominences are being specified. Treatment 2 consisted of RAR antagonist administration over a later phase from stage 29/30-40, (35-66 hpf ; Fig. 1a). During this time the face is growing and facial structures such as jaw cartilage are specified. 100% of the embryos treated with the RAR inhibitor during the early treatment phase developed a median cleft whereas 91% of embryos developed a median cleft with RAR inhibition during the later treatment phase (Fig. 1b-g; n = 22/treatment, 2 biological replicates). The median cleft produced by each treatment was assessed by investigating the shape of the mouth. The mouth phenotypes existed as a spectrum, at one end of this spectrum was a median cleft that came to a point, resulting in a triangular shaped mouth (eg Fig. 1c arrow). This phenotype was more prevalent in embryos after treatment 1 (68% of treated embryos). The other shape was where the midline of the dorsal aspect of the mouth was flattened (eg Fig. 1g arrow). This latter phenotype was more common in embryos after treatment 2 (64% of treated embryos).
In order to quantify these differences in the mouth shape, as well as the surrounding midface, we performed a geometric morphometric analysis of embryos treated with the RAR inhibitor. Thirty-three standardized landmarks were used to label the orofacial region for morphometric analysis (Additional file 1: Figure S1B, [14,22]). Landmark data were aligned by Procrustes fit to eliminate information about size and analyzed using multivariate statistical analyses. First, a discriminant function analysis was utilized to visually assess the changes in landmark position. The results are represented as vector diagrams, such that the flat end of the vector represents the average landmark position in control embryos, and the closed circle represents the average landmark position in the treated embryos. The length of the vectors represents the magnitude of shape change. (Treatment 1 Fig. 1h, Treatment 2 Fig. 1i). These vectors are presented on a transformation grid; warping of the grid highlights regions with the most significant changes. The transformation grids from both the early and late treatments revealed significant shifts that were consistent with a narrowing of the face, especially just above the mouth (Fig. 1h, i). The dramatic warping of the grid in the midface region also exhibited a narrowing of the face. Despite the overall similarity between the changes in landmark position in each treatment phase, there were also important differences. The direction and length of the landmarks just above the mouth were different in embryos after the two different treatments (Fig. 1h, i red arrowheads). This difference could reflect changes in the amount and distribution of primary palate tissue. The landmarks around the mouth were angled differently in the embryos from the two treatment phases, indicative of the spectrum of mouth shape changes that was qualitatively observed (Fig. 1h, i).
Statistical relationships across the RAR antagonist treatment and control groups were visualized simultaneously using a canonical variate analysis (CVA). Data from the CVA is represented as a scatter plot; where each dot represents the average landmark position of a single embryo in a given treatment (Fig. 1j). Procrustes distances (PDs) measure how different one treatment group is from another, with significance indicated by p-values. The CVA confirmed the findings of each individual DFA; there were significant differences between the control and RAR antagonist treated groups (Treatment 1 Procrustes distance: 0.2967, p < 0.0001, Treatment 2 Procrustes distance: 0.2719, p < 0.0001). There was no difference between the two control groups, indicated by the overlap of the circles encompassing average landmark position of the embryos in each condition (Fig. 1j, light blue = control for treatment 1, light green = control for treatment 2;PD = 0.0485, p = 0.028). The two different RAR antagonist treatment groups (T1 RAR Inhib, dark blue and T2 RAR Inhib, dark green) appeared separated on the CVA graph (Fig. 1j, PD = 0.0749, p = 0.004). However, the Procrustes distance is low and p-value is high reflecting only subtle differences in the shape of the orofacial region after early and late RAR antagonist treatments.
Morphometric analysis provided an assessment of the orofacial shape changes but did not address changes in orofacial size. To investigate changes in size, we measured the intercanthal distance, the distance between the eyes. Intercanthal distance was significantly decreased when retinoic acid signaling was inhibited compared to  (87-92 hpf) embryos. The mouths are outlined in red dots and the shape of the mouth is shown in red in the left corners. b) Embryos treated with 1% DMSO from stage 24 to stage 29/30 (26-35 hpf (35-66 hpf). Black arrow indicates a mouth shape with a flat dorsal aspect at the midline h-i Transformation grids exhibiting changes in orofacial landmarks during early (h) and late (i) RAR Inhibition. Flat end of the vector represents average position of control landmark; closed circle represents average position of landmark in treated embryos. j CVA analysis scatter plot. k Intercanthal distance relative to control. For all experiments n = 22/treatment, 2 biological replicates controls during both treatment windows (Fig. 1k, p < 0.0001, n = 22/treatment with 2 biological replicates). The intercanthal distance was not different between the early and late phase RAR antagonist treated groups (Fig. 1k, p = 0.19).
Together our analyses suggested a continuous requirement for retinoic acid during early and late phases of orofacial development. The spectrum of the mouth shape changes shifted from early to late RAR inhibition. This change was reflected in the subtle quantitative shape differences while the size of the face was the same. Further analyses are required to determine how these changes reflect the evolving role of RAR during orofacial tissue development.
Large scale transcriptome analyses revealed highly altered genes after early and late RAR inhibition To better understand the underlying mechanisms driving the changes in the shape of the face when embryos have decreased RAR function, we performed two global unbiased expression analyses. For each analysis, the facial region was dissected at the conclusion of the treatment window (Fig. 1a, Additional file 1: Figure S1A). After early phase RAR inhibition, a microarray expression platform assessed relative changes in RNA (Fig. 1a, n = 2 biological replicates). An RNA-seq analysis was used to assess relative changes in RNA expression derived from microdissected face tissues after late phase RAR inhibition (reviewed in Fig. 1a, n = 3 biological replicates). A random subset of genes were validated by quantitative RT-PCR (Additional file 2: Figure S2); others have been validated in the course of other work in our laboratory [14]. It is important to note that two different methods to evaluate gene expression were utilized, during different years, and therefore the specific genes from each list are not directly comparable.

Top genetic alterations after early inhibition of RAR
The most highly altered genes after early RAR inhibitor treatment (Table 1) were first examined. Genes that encode transcription regulators, RNA processing proteins, and calcium signaling mediators, such as CAMKIIG, were among the most decreased transcripts. Genes that were most enriched included those encoding proteins involved in heart and pancreas function and development, suggesting that RA signals inhibit more posterior tissues from developing in the face.

Top genetic alterations after late inhibition of RAR
The most highly altered genes after the late RAR inhibitor treatment (Tables 2 and 3) were also examined. Transcription factors and TGF-beta/BMP family members were among the top decreased genes. In addition, genes in the retinoic acid signaling pathway were also identified in this analysis. Genes with increased expression after late inhibition of RAR included those encoding proteins associated with cells of the eye and the vasculature and heart. Again this result suggests that retinoic acid signals are required to restrict the expansion of surrounding or posterior tissues.

Functional classification of altered genes after early and late RAR inhibition
We utilized the online analysis software DAVID to categorize the genes perturbed by RAR inhibition into functional processes (version 6.8, NIH [23,24]. Gene lists were divided into those that were decreased and those that were increased relative to the controls after early and late phase RAR antagonist treatment. After early RAR antagonist treatment, the top functional categories from depleted gene list were "Transcription" and "Cell Proliferation and Growth" (Fig. 2 Ai). Categories important for the developing embryo "Organ Development" and "Cell Differentiation" were also highly represented in this list. Other functional categories included: "RNA splicing/Processing", "Intracellular Protein Transport" and "Chromatin Remodeling". The top functional categories from the enriched gene list after early RAR antagonist treatment were also "Transcription" and "Cell proliferation and Growth" as well as "Cell Differentiation" (Fig. 2 Aii). After late RAR antagonist treatment, the top functional categories from depleted gene list were "Transcription", "Cell Proliferation and Growth" and "Organ Development" (Fig. 2 Bi). Other functional categories included: "Extracellular Matrix/Adhesion", "Metabolic Pathways" and "Cartilage Development". The top functional categories from the enriched gene list after late RAR antagonist treatment were "Cell Differentiation", "Organ Development" and "Signal Transduction" (Fig. 2 Bii).
Next, we sought to compare the functional categories identified after each treatment phase. Note, that we comparing the functional groups identified, not the individual genes that were altered. In this analysis we focused on the functional categories representing more than 10% of the genes that were altered after RAR inhibition. This comparison revealed that "Transcription" and "Cell Proliferation/Growth" were both decreased after early and late RAR inhibition (Fig. 2 Ci). This overlap indicated that retinoic acid may be required for these processes throughout orofacial development. Additionally, this finding could offer an explanation for the phenotypic similarities observed after the two treatments ( Fig.  1) For example, RAR antagonist treatment during both phases resulted in smaller faces that were not statistically different. Such similar perturbation in the size of the face could be explained by the common effect on genes Table 3 Genes dysregulated after early RAR inhibition associated with transcriptional regulation. Genes shown in Fig. 3 along with their general function (adapted from genecards.org)

Gene Name General Protein Function
Decreased ATRX α-thalassemia mental retardation X-linked protein A SWI/SNF family member protein involved in chromatin remodeling.

ESR1
Estrogen Receptor 1 A nuclear receptor for the estrogen ligand.

EZH2
Enhancer Of Zeste 2 Polycomb Repressive Complex 2 Subunit A polycomb group family member that methylates Lys-9 and Lys-27 of histone H3, leading to transcriptional repression.

HDAC1
Histone Deacetylase 1 A histone deacetylase that targets lysine residues on the N-terminal part of the core histones (H2A, H2B, H3 and H4) which is necessary for transcriptional repression.

SMARCAD1
SWI/SNF-Related, Matrix-Associated Actin-Dependent Regulator Of Chromatin, Subfamily A, Containing DEAD/H Box 1 A SNF subfamily of helicase proteins that mediates histone H3/H4 deacetylation that is necessary for restoration of heterochromatin organization and propagation of epigenetic patterns following DNA replication.

SUPT6H
SPT6 Homolog, Histone Chaperone A transcription elongation factor which binds to histone H3 and plays a key role in the regulation of transcriptional elongation and mRNA processing.

SUZ12 Polycomb Repressive Complex 2 Subunit
A Polycomb group protein that methylates Lys-9 and Lys-27 of histone H3, leading to transcriptional repression.

UHRF1
Ubiquitin Like With PHD And Ring Finger Domains 1 A protein that bridges DNA methylation and chromatin modification.
BCL2 Associated Athanogene 6 A nuclear protein that forms a complex with E1A binding protein p300 and is required for the acetylation of p53 in response to DNA damage.

BRD4
Bromodomain Containing 4 A chromatin reader protein that recognizes and binds acetylated histones and plays a key role in transmission of epigenetic memory across cell divisions and transcription regulation.

CEBPB
CCAAT/Enhancer Binding Protein Beta A transcription factor that contains a basic leucine zipper (bZIP) domain.

FOS
Fos Proto-Oncogene, AP-1 Transcription Factor Subunit A leucine zipper protein that can dimerize JUN and together form part of the AP-1 transcription factor complex.

GATA3
GATA Binding Protein 3 A protein that belongs to the GATA family of transcription factors.

JUN
Jun Proto-Oncogene, AP-1 Transcription Factor Subunit A protein that can dimerize FOS and together form part of the AP-1 transcription factor complex.

MTA2
Metastasis Associated 1 Family Member 2 A component of NuRD complex which is a deacetylase that is required for remodeling nucleosomes.

NCOA3
Nuclear Receptor Coactivator 3 A histone acetyltransferase that recruits p300/CBP as well as CREB binding proteins to form a transcriptional coactivation complex.

NCOR1
Nuclear Receptor Corepressor 1 A protein that is associated with thyroid-hormone and retinoic-acid receptors and promotes chromatin condensation and transcriptional repression.

OTUB1
OTU Deubiquitinase, Ubiquitin Aldehyde Binding 1 A hydrolase that can specifically remove Lys-48-linked conjugated ubiquitin from proteins and thereby prevents protein degradation.

PAX7
Paired Box 7 A member of the paired box (PAX) family of transcription factors.

PBRM1
Polybromo 1 A protein necessary for ligand-dependent transcriptional activation by nuclear hormone receptors.

WBP2
WW Domain Binding Protein 2 A WW domain binding protein that is a transcriptional coactivator of the estrogen and progesterone receptor. regulating "Cell Proliferation". On the other hand, several functional categories did not overlap; this could offer insight into the differences in the phenotype after each treatment. "Cell Death" was a category of genes enriched after the late RAR inhibitor treatment but not the early treatment ( Fig. 2 Ci). We could hypothesize that deficient RAR function during the later time points could result in excess apoptosis which could in turn cause a larger gap in the palate region. These analyses provide new information about the potential roles of RA signaling during craniofacial development and allow for the formation of new hypotheses that can be tested in future work.

Network analysis centered around the functional classification and highly dysregulated genes
While DAVID analysis identified functional categories perturbed after RAR inhibition, it did not provide any information about how signals, pathways, transcription factors and epigenetic modifiers might be integrated in orofacial tissues. Therefore, we utilized Ingenuity Pathway Analysis (version 01-07, Qiagen) to better dissect networks of interacting proteins that were perturbed upon RAR inhibition. Ingenuity Pathway Analysis is a platform that synthesizes literature from human, mouse, and cell culture research to build a database of functional interactions. We utilized this database to build networks between genes that were dysregulated during early or late RAR inhibition, with a focus on two of the top functional categories identified by DAVID analysis Transcription and Organ Development.
Networks of genes were identified after early and late RAR inhibition that contained genes from the functional category of "transcription" Transcription was identified as the largest DAVID functional category of altered genes, after both early and late RAR inhibition. We utilized the previously characterized transcription pathways in Ingenuity to build networks around the most highly altered genes in the transcriptional functional categories after each treatment. For Venn diagrams representing the overlap between decreased genes (i) and increased genes (ii) in early and late RAR inhibition. Blue and orange represent functional categories altered with early RAR inhibition; green and purple represent functional categories altered with late RAR inhibition each network, we limited the number of genes to 30 and focused on genes that were altered, minimizing the inclusion of non-affected intermediaries.
Transcriptional network altered in early orofacial development contained chromatin regulators and modulators of RAR Among the most highly depleted transcriptional regulators altered the early RAR inhibitor treatment was Chromodomain helicase DNA binding protein 1 (CHD1) (see Table 1). CHD1 belongs to a family of proteins that modify the chromatin and regulate transcription during development [25]. This network linked CHD1 to several other epigenetic regulators also altered by RAR inhibition, such as HDAC1, EZH2, ATRX and CHD2 (Fig. 3a, Table 3). A subset of the genes from this network encodes proteins that are coactivators or repressors of retinoic acid receptors. For example, NCOR1, EZH2 (also called PRC2) and HDAC1 encode proteins that belong to complexes that have been shown to repress RAR transcription, while NCOA encodes a protein that is a co-activator of retinoic acid [26][27][28][29][30][31].
Overall, this analysis revealed that the transcriptional regulators that were altered after early RAR inhibition are modulators of chromatin and RAR function.
Functional validation: Chromatin modifiers CHD1 and HDACs are required for craniofacial development in Xenopus embryos We sought determine whether the chromatin remodelers identified in our analyses, CHD1 and HDAC1 are required for craniofacial development in Xenopus. To accomplish this, we decreased the function of these proteins and examined the effects on facial morphology. A custom designed antisense oligonucleotide morpholino was used to investigate the function of CHD1 in orofacial development. Embryos were injected at the 1 cell stage with a CHD1 morpholino or a standard control morpholino; craniofacial phenotypes were assessed at stage 43 (87-92 hpf ). CHD1 morpholino injected embryos had craniofacial defects including midface hypoplasia, misshapen mouth and protruding lenses (n = 33, 3 biological replicates, Fig. 3b-c). Further validation and analysis of this protein is provided in another study (Wyatt and Dickinson, in preparation).
To investigate the role of HDACs during early Xenopus orofacial development a well-characterized pharmacologic inhibitor, Trichostatin A was utilized [32]. Embryos were exposed to 100 μM Trichostatin A from stage 24-30, (26-35 hpf), transferred to normal culture media after treatment, and assessed for facial phenotype at stage 43 (87-92 hpf n = 20/treatment, 3 biological replicates). HDAC inhibitor treated embryos had craniofacial phenotypes including: small eyes, lens defects, midfacial hypoplasia, and a small, round mouth when compared to controls (Fig. 3d-e).
In summary, inhibition of either CHD1 or HDACs during the early phase of orofacial development resulted in defects including a narrower midface and smaller mouth, similar to what was observed in embryos with reduced RAR function. Together, these data indicate the possibility that retinoic acid signaling may work with CHD1 and HDACs during early orofacial development.
Transcription network identified with late orofacial development included transcription factors Genes identified by DAVID analyses as transcriptional regulators after late RAR inhibition were used to create a network with Ingenuity Pathway Analysis software ( Fig. 4 and Table 4). The majority of genes identified in this network were transcription factors. Genes that were depleted included members of the Forkhead (FOXP4, FOXD1 and FOXC1) and Aristaless families (ALX1 and ALX4), as well as HOXA4, SIX2 and GATA2. Other transcription factors were enriched in tissues deficient in RAR function, including DLX2 and LHX1. These results reveal that RA signaling is required for the proper expression of transcription factors in the developing face. Networks of genes were identified after early and late RAR inhibition that contained genes from the functional category of "organ development" Organ development was a functional category highly represented in the DAVID analysis (Fig. 2). We chose these categories as our next target for Ingenuity Pathway Analysis because genes in this category are integral to developmental processes including facial morphogenesis. Custom networks were built, focusing on the most highly altered genes in the categories of organ development and differentiation. The number of genes in each network was limited to 30 and unaffected intermediaries were minimized.
Calcium signaling genes were dysregulated after early RAR inhibition Calcium/calmodulin dependent protein kinase 2-isoform gamma (CAMKIIg) was highly depleted after early RAR inhibition and was present in the functional category of "Organ Development" in the DAVID analysis. This gene is known to be important for development, which led us to construct an Ingenuity network around CAMKIIg (Fig. 5). Predictably, this network identified other proteins involved in calcium signaling such as the calcium channel, Inositol 1,4,5-Trisphosphate Receptor (ITPR) and Calcineurin A (PP3CA). This network also included NFAT and DYRK1A, which are misregulated in Trisomy 21 ( Fig. 5a and Tables 5).
Functional validation: CAMKII and CALCINEURIN are required for craniofacial development in Xenopus embryos We validated the findings of this network by examining altered levels of two proteins identified in the Ingenuity analysis, CAMKIIg and Calcineurin A (also called PPP3CA). To determine if these proteins were required during the early orofacial development phase, pharmacological inhibitors were utilized to diminish their function.
The antagonist KN93 was used to disrupt CaMKII function. KN93 prevents the association between calmodulin and CaMKII and has been widely used in experimental models, including zebrafish [33]. Embryos were exposed to 100 μM KN93 from stage 24-30 (26-35 hpf ), the same time window as our early RAR inhibitor treatment (reviewed in Fig. 1a, n = 20, 2 biological replicates). After treatment, embryos were raised until stage 43 (87-92 hpf ) and then facial phenotypes were assessed. KN93 treated embryos had a narrower midface, close-set eyes and protruding lenses (Fig. 5b-c). The shape of the mouth resembled those treated with the RAR inhibitor with a median cleft that either came to a point (not shown), or formed a gap (Fig. 5c).
In summary, inhibition of either CaMKII or Calcineurin during the early phase of orofacial development resulted in defects similar to RAR inhibitor treated Xenopus embryos. These results suggest that RA signaling may interact with calcium signaling during orofacial development.
BMP and WNT family members were dysregulated after late RAR inhibition BMP5 and the Wnt  Table 4 Genes dysregulated after late RAR inhibition associated with transcriptional regulation. Genes shown in Fig. 4 along with their general function (adapted from genecards.org)

Gene
Gene Name General Protein Function

Decreased ALDH1A2
Aldehyde Dehydrogenase 1 Family Member A2 A member of the aldehyde dehydrogenase family of enzymes that catalyzes the conversion of retinaldehyde to retinoic acid, the active derivative of vitamin A.

ALX1
ALX Homeobox 1 Unknown in humans, in rodents it is necessary for survival of forebrain mesenchyme and possibly neural tube development.

ALX4
ALX Homeobox 4 Homeodomain transcription factor expressed during development of bones, limbs, hair, teeth, and mammary tissue.

BMP2
Bone Morphogenetic Protein 2 Secreted ligand of the TGF-beta family that is important in bone and cartilage development.

FOXP4
Forkhead Box P4 A member of the forkhead box family of transcription factors that is a transcription repressor.

GATA2
GATA Binding Protein 2 Zinc-finger transcription factor that regulates genes involved in development and proliferation of hematopoietic and endocrine cell lineages.

HAS2
Hyaluronan Synthase 2 Enzyme involved in the synthesis of hyaluronic acid.
HMGA2 High Mobility Group AT-Hook 2 A member of the non-histone chromosomal high mobility group protein family that acts as a transcription regulator.

HOXA4
Homeobox A4 A member of the homeobox family of proteins that are transcription factors, important for gene expression, morphogenesis, and differentiation.

Phenylalanine Hydroxylase
Enzyme that catalyzes the rate-limiting step of phenylalanine catabolism.

Retinoic Acid Receptor Beta
Ligand-dependent transcription regulator important for embryonic morphogenesis, cell growth, and differentiation.

RARG
Retinoic Acid Receptor Gamma Ligand-dependent transcription regulator implicated in embryonic growth, limb bud development, and differentiation.

SFRP1
Secreted Frizzled Related Protein 1 Soluble modulator of Wnt signaling; important for cell growth and differentiation.

SIM2
Single-Minded Family BHLH Transcription Factor 2 Transcription factor that serves as a master regulator of neurogenesis, located within Down's Syndrome Critical Region.

SIX2
SIX Homeobox 2 A member of the sine oculis family of proteins that are transcription factors involved in limb and eye development.

SSTR2
Somatostatin Receptor 2 Receptor that mediates the inhibitory effect of somatostatin on hormones and secretory proteins. Also involved in neurodevelopment.

WNT1
Wnt Family Member 1 Soluble transcription regulator that binds to frizzled receptors and plays several roles in development.

Increased DLX2
Distal-Less Homeobox 2 Transcription factor expressed in the developing head and limbs; important for the terminal differentiation of neurons and craniofacial patterning.

LHX1
LIM Homeobox 1 Transcription factor important for development and differentiation of the neural crest, ectoderm, and urogenital system.

VSX1
Visual System Homeobox 1 Transcription factor that regulates opsin gene expression.

UNCHANGED
CSNK2A2 Casein Kinase 2 Alpha 2 Protein kinase that phosphorylates acidic C-terminal residues to regulate cell cycle progression, apoptosis, transcription, and viral infection. antagonist FRZB were among the most highly changed genes after late RAR inhibition. These genes were also identified in the functional category of "Organ Development". An Ingenuity network was constructed around these genes and identified additional members of the WNT and BMP pathways, including WNT1, BMP7, as well as TGFB1 a close relative of BMP. This network also connected genes associated with cartilage (COL9A, and BGN) and transcription factors ALX1, SIX2 and DLX2 (Fig. 6a, Table 6).
BMP and WNT signals are required for craniofacial development in Xenopus We functionally tested whether BMP and WNTs are required for craniofacial development during the late orofacial development phase. To accomplish this, we utilized pharmacological inhibitors that could specifically test the role of these morphogens during this phase of orofacial development, without perturbing earlier events in development.
To disrupt BMP signaling, Type I BMP receptors, which are activated by BMP2 and BMP5 were blocked utilizing an established pharmacological antagonist [35,36]. Embryos were treated with 100 μM of LDN193189 from stage 30-40 (35-66 hpf, n = 20, 3 biological replicates, Fig. 6d-e). Facial characteristics were assessed at stage 43 (87-92 hpf ). The midface appeared narrow, and the mouth was both smaller and rounder (Fig. 6e). These results demonstrate that BMP is required during late phase orofacial development in Xenopus.
In summary, these experiments revealed that both BMP and Wnt are required over a very specific period of orofacial development.

Genes dysregulated in Xenopus after RAR inhibitor treatments are associated with human craniofacial birth defects
In Xenopus, RAR inhibition during specific windows of development result in median clefts in the primary palate [14]. While this phenotype resembles the human median cleft, it remained to be determined if the phenotypes had similar genetic origins. We explored this question by investigating whether genes altered in our Xenopus cleft phenotype were associated with a human craniofacial defect. Using the Online Mendelian Table 4 Genes dysregulated after late RAR inhibition associated with transcriptional regulation. Genes shown in Fig. 4

along with their general function (adapted from genecards.org) (Continued) Gene
Gene Name General Protein Function expression, morphogenesis, and differentiation.

HOXD11 Homeobox D11
A member of the homeobox family of proteins that are transcription factors; important for morphogenesis; limb and genital morphogenesis in particular.

RARA
Retinoic Acid Receptor Alpha Ligand-dependent transcription regulator implicated in development, differentiation, apoptosis, and transcription of clock genes.

SRA1
Steroid Receptor RNA Activator 1 Transcriptional repressor that binds to non-coding RNAs and has roles in metabolism, adipogenesis, and chromatin organization.

Radixin
Cytoskeletal protein that important in linking actin to the plasma membrane.

SPTBN2
Spectrin Beta, Non-Erythrocytic 2 Component of spectrin, important in securing glutamate transporter EEAT4 to the plasma membrane and regulating glutamatergic signaling.

Increased AK1
Adenylate Kinase 1 Catalyzes the reversible transfer of the terminal phosphate group between ATP and AMP BET Bromodomain containing, BET Family of genes that contain 110 amino acid domains that recognize acetylated lysine residues and translate the acetylation signal into various phenotypes depending on the family member's location and expression.

Basic Helix-Loop-Helix Family Member E40
Basic helix-loop-helix protein expressed in various tissues and is believed to be involved in the control of circadian rhythm and cell differentiation.

BRD4
Bromodomain Containing 4 Protein with a conserved sequence motif involved in chromatin targeting and binds acetylated lysine residues.

CKMT1A/ CKMT1B
Creatine Kinase, Mitochondrial 1A, Creatine Kinase, Mitochondrial 1B Serine/threonine kinases that play a central role in energy transduction in tissues with large, fluctuating energy demands, such as skeletal muscle, heart, brain and spermatozoa.

COL2A1
Collagen Type II Alpha 1 Chain Alpha-1 chain of type II collagen, a fibrillar collagen found in cartilage

DYRK1A
Dual Specificity Tyrosine Phosphorylation Regulated Kinase 1A Located on the Down's Syndrome Critical Region, this gene product is involved in cell survival, proliferation, quiescence, and differentiation.

Fos Proto-Oncogene, AP-1 Transcription Factor Subunit
Leucine zipper protein that dimerizes with proteins of the JUN family to form the transcription factor AP-1; regulates cell proliferation, differentiation, and transformation.

FOXD1
Forkhead Box D1 Member of the forkhead family of transcription factors, functions in the differentiation of nephron progenitors and the development of the retina.

GATA5
GATA Binding Protein 5 Transcription factor containing 2 GATA-type zinc fingers, involved in establishment of cardiac smooth muscle diversity.

BRD4-interacting chromatin domain remodeling complex associated protein
Plays a role in BRD4-mediated transcription.

GPSM1
G Protein Signaling Modulator 1 A receptor-independent activator of G protein signaling, regulates spindle orientation and cortical neuronal cell fate. May also be involved in macroautophagy in the intestine.

HOXA11
Homeobox A11 Transcription factor that provides cells with specific positional identities on the anterior-posterior axis.

HOXD11
Homeobox D11 Transcription factor that provides cells with specific positional identities on the anterior-posterior axis.

ID2
Inhibitor Of DNA Binding 2 Basic helix-loop-helix transcription regulator that binds to other HLH transcription regulators and inhibits their activity; important for glial differentiation and chondrocyte development.

PP2A (PPP2CA)
Protein phosphatase 2 catalytic subunit alpha Catalytic subunit of serine/threonine phosphatase involved in the negative regulation of cell growth.

Peroxisome Proliferator Activated Receptor Gamma
Forms heterodimers with retinoid X receptors (RXRs) to regulate transcription.
Inheritance of Man database (OMIM.org), we searched for 6 keyword terms that applied to symptoms of a defect in the midface region (listed in Fig. 7). 406 human genes were associated with these keywords (Additional file 3: Table S1). We then compared the list of human genes to the significantly altered genes after both early and late RAR inhibition (Additional file 3: Table S1). Of the 415 genes that were significantly altered after early RAR inhibition, 19 were associated with a median orofacial defect (4.6%, Fig. 7, Table 7A). After late RAR inhibition, 9 of the 93 genes that were significantly altered were also associated with a median orofacial defect (9.7%, Fig. 7, Table 7B). These results suggest that retinoic acid could also regulate a number of genes that can cause midface birth defects in humans.

Discussion
Here we combine transcriptional expression analysis, functional validation and comparison with human databases to provide a global view of the effects of retinoic acid perturbation and clefting in a model vertebrate.

Retinoic acid requirement during early orofacial development
We defined "early" orofacial development to include the time when cranial neural crest cells migrate and facial prominences are specified. Also during this period, the first steps in embryonic mouth occur and the cranial placodes are formed.

Epigenetic regulation
Generally during development, retinoic acid is an important morphogen that directs gene expression that is in turn necessary for specifying cell fates. Therefore, it would not be surprising that RA is also required for maintaining the expression of transcriptional machinery, epigenetic factors and cofactors required for gene expression during orofacial specification. Consistent with this hypothesis, we identified a network of genes dysregulated in facial tissues deficient for RAR function, which included chromatin readers and remodelers. For example, two Chromodomain helicase DNA-binding (CHD), CHD1 and CHD2, were decreased in Xenopus embryos exposed to a RAR antagonist. Further, functional validation revealed that CHD1 is required for orofacial development in Xenopus. The CHD family of enzymes belong to the SNF2 superfamily of ATP-dependent chromatin remodelers and are associated with active transcription (reviewed in [41]). Complete loss of CHD1 in mice is devastating; embryos die before craniofacial development can occur [42]. Other CHD family members, specifically CHD7 and CHD8 are essential for craniofacial development in mice and humans [43][44][45]. For example, CHD7 is associated with CHARGE syndrome, which can include craniofacial defects such as a small mouth, and cleft lip and/or cleft palate and are the result of problems in neural crest development [44,45] In addition to CHD family members, our analysis uncovered other chromatin readers and remodelers including ATRX [46], BRD4 (MGI) and MTA2 [47] that are also associated with craniofacial birth defects. Current work in the lab is investigating the role of CHD1 and other epigenetic regulators during Xenopus craniofacial development. RAR and other retinoic acid pathway members are expressed in the face throughout both early and later orofacial development in Xenopus. We might expect, therefore, that during early orofacial development it is necessary to maintain and perpetuate retinoic acid signaling. To accomplish this, RAR function could be critical for the expression of its partners. Indeed our analysis identified a network that included chromatin remodelers and cofactors involved in RAR receptor mediated transcriptional regulation. Traditionally, upon retinoic acid binding, histone modifying complexes are recruited to target genes to promote their transcription. Carbohydrate binding protein subunit of the protein phosphatase 1 (PP1) complex. Involved in the maintenance of glycogen levels by activating/inhibiting synthesis and breakdown enzymes.

PPP2C
Protein Phosphatase 2 Catalytic Subunit Catalytic subunit of serine/threonine phosphatase involved in the negative regulation of cell growth.

PRKAR1A
Protein Kinase cAMP-Dependent Type I Regulatory Subunit Alpha Regulatory subunit of the cAMP-dependent protein kinases involved in cAMP signaling in cells including GPCR signal transduction.

Serum Response Factor
Ubiquitous transcription factor that stimulates proliferation and differentiation.

TAL BHLH Transcription Factor 1, Erthyroid Differentiation Factor
Basic helix-loop-helix transcription factor that plays a role in hematopoietic and erthyroid differentiation.

Inositol 1,4,5 triphosphate receptor family
Family of membrane glycoproteins that act as calcium channels to regulate several cellular processes including proliferation, learning, and memory.
Notably, we determined that one such histone modifying complex member, NCOA3, was dysregulated in the face tissues of Xenopus embryos treated with and RAR antagonist. On the other hand, in the absence of retinoic acid, members of several repressive complexes are recruited to target genes to promote transcriptional repression. Several such complex members including HDAC1 EZH2, SMAR-CAD1, SUZ12, and NCOR1, were dysregulated in embryos with deficient RAR function. Functional validation demonstrated a distinct role for HDACs during the early phase of orofacial development, which is supported by work in other model vertebrates [48][49][50].
Together our results suggest that RAR function is required for maintaining the proper expression levels of a number of epigenetic regulators necessary for retinoic acid mediated gene expression. A role for retinoic acid in regulating DNA methylation, histone modification, and miRNA expression has also been demonstrated during development, stem cell maintenance and in cancer (reviewed in [51]). Thus, while a role for retinoic acid in regulating epigenetic factors is not unexpected, a better understanding of the coordination of RAR and chromatin remodelers will provide a more complete view of how this important morphogen gradient is utilized during the development of the face.

Calcium-NFAT pathway
Network analysis of genes dysregulated in the orofacial tissues after early RAR inhibition identified several connected genes encoding proteins associated with calcium signaling such as the calcium channel ITPR, CAMKIIg and Calcineurin A. Functional validation demonstrated that CaMKII and Calcineurin are indeed required during early orofacial development in Xenopus. Calcium levels in the cell mediate a multitude of processes and events during development (reviewed in [52]). The importance of proper calcium regulation insures that there are several modulators of calcium levels, depending on the tissue type and process. For example, calcium can be released from internal stores upon binding of inositol triphosphate to a receptor that acts as a calcium channel (ITPR). Calcium levels then modulate calcium sensitive kinases such as Ca 2+ /calmodulin-dependent protein kinases (CaMK) and phosphatases, such as calcineurin. These enzymes regulate the activity of a wide range of signaling and transcription factors in the embryo (see [52]). One such transcription factor, identified in our network analysis, is nuclear factor of activated T cells (NFAT) [53]. Calcineurin dephosphorylates this protein, which results in its translocation to the nucleus. Once in the nucleus, NFAT can activate the transcription of many different genes that in turn influences diverse cellular functions including diversification of neural crest cells and craniofacial development [54][55][56][57] . In Xenopus, nfatc3 is expressed in tissues consistent with the neural crest during early orofacial development (Xenbase). Another regulator of NFAT activity that was identified in our network analysis of genes dysregulated after RAR inhibition was DYRK1A. This protein phosphorylates NFAT triggering its movement into the cytoplasm and thus preventing its transcriptional activity. Haploinsufficiency of DYRK1A results in craniofacial defects that can include cleft lip/palate in humans [58,59]. Further, DYRK1A is located on the Down's Syndrome Critical Region of chromosome 21 and has been linked to the characteristic facial features associated with this disorder [60]. Moreover, the craniofacial features associated with increased DYRK1A dosage closely resemble those of NFAT and Calcineurin loss of function in mice [57]. The Calcineurin-NFAT-DYRK1A axis has also been shown to interact with CAMKIIg. CAMKII inhibits NFAT nuclear translocation by phosphorylation and inhibition of Calcineurin in cardiac myocytes [61]. While we know that CAMKII is integral to embryonic development, including signaling in cilia [62][63][64], a role in craniofacial development has never been demonstrated. Therefore, a novel avenue for future work will be to investigate the roles of CAMKII, Calcium signaling and NFAT regulation and how these interact with retinoic acid signaling during orofacial development.
In summary, we demonstrate that retinoic acid is required for the expression of epigenetic factors and members of Calcium-NFAT signaling which in turn could influence events such as neural crest migration and specification of the facial prominences during early orofacial development.
The requirement for retinoic acid during late orofacial development We defined late orofacial development to include a time when facial prominences grow and converge as well as the differentiation of the facial tissues. During this time Table 6 Genes dysregulated after late RAR inhibition associated with organ development and differentiation. Genes shown in Fig. 6 along with their general function (adapted from genecards.org)

Gene
Gene Name General Protein Function

Decreased ALX1
Alx Homeobox 1 Unknown in humans, in rodents it is necessary for survival of forebrain mesenchyme and possibly neural tube development.

Biglycan
Small, leucine-rich proteoglycan that plays a role in bone growth, muscle development, and collagen fibril assembly.

BMP2
Bone Morphogenetic Protein 2 Secreted ligand of the TGF-beta family that is important in bone and cartilage development. A member of the non-histone chromosomal high mobility group protein family that acts as a transcription regulator.

RARB
Retinoic Acid Receptor Beta Ligand-dependent transcription regulator important for embryonic morphogenesis, cell growth, and differentiation.

Retinoic Acid Receptor Gamma
Ligand-dependent transcription regulator implicated in embryonic growth, limb bud development, and differentiation.

SIX2
SIX Homeobox 2 A member of the sine oculis family of proteins that are transcription factors involved in limb and eye development.

TGFB1
Transforming Growth Factor Beta 1 Ligand for TGF-beta superfamily of signaling molecules, activates Smad signaling to regulate cellular proliferation and differentiation.

WNT1
Wnt Family Member 1 Soluble transcription regulator that binds to frizzled receptors and plays several roles in development.

Increased DLX2
Distal-Less Homeobox 2 Transcription factor expressed in the developing head and limbs; important for the terminal differentiation of neurons and craniofacial patterning. Chordin family member that interacts with TGF-beta signaling and may play a role in myoblast and osteoblast differentiation.

SIX5
SIX Homeobox 5 Transcription factor involved in the regulation of organogenesis.

SNAI1
Snail Family Transcriptional Repressor 1 Zinc finger transcriptional repressor that downregulates the expression of ectodermal genes in the mesoderm.

TGFB2
Transforming Growth Factor Beta 2 Ligand for TGF-beta superfamily of signaling molecules, activates Smad signaling. the embryonic mouth is formed and jaw cartilage and muscle begin to form.

WNT/beta-catenin and BMP
Our work suggests that during the later period of orofacial development, retinoic acid signaling cooperates with well-known regulators of craniofacial development, WNT and TGF-beta family members including BMP. Ingenuity pathway analysis identified a network of dysregulated genes connecting WNT1, FRZB, BMP2 and BMP5. Further, functional validation using pharmacological WNT and BMPR antagonists over this late phase of orofacial development resulted in Xenopus embryos with defects in the midface. Specifically, embryos treated with a WNT antagonist developed a large round embryonic mouth and narrower midface. This is consistent with previous work Xenopus that demonstrated that overexpression of Wnt inhibitors frzb-1, crescent and dkk similarly resulted in a larger mouths but narrower faces [65]. Together, work in various vertebrate models indicate that Wnt/beta-catenin inhibitors define the region that will form the embryonic mouth while Wnt/ beta-catenin signals are required in the surrounding facial prominences to maintain proliferative regions and direct jaw cartilage development [66][67][68][69][70]. Thus, disruption of the balance of WNT activation and inhibition can alter the relative facial dimensions in the developing embryo [71]. Next we determined that Xenopus embryos exposed to a BMPR antagonist over the late phase of orofacial development resulted in both a smaller mouth and midface. Consistently, mice deficient for BMP4 during orofacial development have decreased facial growth and reductions in the jaw skeleton (reviewed in [72]{Shuman, 2007 #3380} [73]. In humans, mutations in members of either WNT or BMP pathways are associated with craniofacial defects [74][75][76][77][78] Thus, a role for WNT/beta-catenin and BMP signaling during orofacial development of Xenopus is not surprising [79].
Interactions between WNT, BMP and retinoic acid have also been previously demonstrated in craniofacial development in model organisms and humans. For example, exogenous retinoic acid suppresses canonical Wnt in the developing palate [80,81]. Retinoic acid interacts with BMPs to determine the identity of the frontonasal mass and maxillary prominences in developing chick embryos [82,83]. Additionally, teratogenic levels of retinoic acid and BMP signaling in humans contribute to the severity of holoprosencephaly, a disorder with characteristic midline facial defects [84]. Thus, our study demonstrates that Xenopus, like other vertebrates, likely requires the coordination of retinoic acid, WNT and BMP signals to direct orofacial development.

Transcription factor network
Retinoic acid signaling is well known for its role in regulating the expression of transcription factors, especially those with a homeobox {Langston, 1994 #3084} [18,85]. Therefore, we were not surprised to uncover a network of homeobox containing transcription factors disrupted in orofacial tissues with deficient RAR function. Several of these transcription factors (ALX1, ALX4, FOXC, FOXD1 and SIX2) are required for craniofacial [86][87][88][89][90][91]. Moreover, at least of subset of such transcription factors could be direct targets of retinoic acid signaling. The promoter regions of two genes in our analysis, FOXP4 and HOXA4, were bound by RAR in embryonic stem cells. These genes were also activated upon exogenous retinoic acid application [92]. Additional experiments would be necessary to determine whether a subset of the transcription factors identified in this work are indeed direct targets of retinoic acid Fig. 7 Overlap of altered gene expression in RAR deficient face and human orofacial defects. Venn diagram exhibiting the overlap among genes altered in our treatment paradigm in Xenopus and human median orofacial defects. Blue circle: total number of Xenopus genes with human orthologs (415) that were significantly altered with early RAR inhibition. Green circle: total number of Xenopus genes with human orthologs (93) that were significantly altered with late RAR inhibition. Grey circle: Human genes associated with one of the six craniofacial defects listed (325 total). Overlap of circles denotes genes that appeared in both lists signaling during orofacial development. Regardless, these results present the possibility that retinoic acid signals coordinate either directly or indirectly a host of transcription factors critical for the development and differentiation of the facial tissues.

XENOPUS as a tool for understanding human orofacial birth defects
Twenty-six genes were dysregulated in orofacial tissues after RAR inhibition are associated with specific birth defects affecting the midface (Table 7, Fig. 7). Humans Table 7 Genes associated with human medial craniofacial defects that were dysregulated with RAR inhibition in X. laevis. A-B: Human orthologs of Xenopus genes significantly altered with early RAR inhibition (A) or late RAR inhibition (B) that were associated with a craniofacial defect (compiled with the use of OMIM.org) with these defects present with orofacial abnormalities that share features observed in Xenopus embryos deficient in retinoic acid signaling. For example, birth defects including Stickler Syndrome and Frontonasal Dysplasia include midline or median clefts, a hallmark of RAR deficiency in Xenopus. Therefore, in humans just as in Xenopus, retinoic acid signals may coordinate the expression of many of the genes required for the formation of the midface and associated tissue types such as cartilage and bone. This could mean that in cases where human midface defects are a result of changes in the dosage of RA-regulated genes, RAR agonists and vitamin A supplementation could potentially serve as future therapies. Our work demonstrates that Xenopus could be an ideal organism in which to test the therapeutic potential of retinoic acid analogs for preventing or ameliorating craniofacial birth defects.

Conclusions
This study utilized Xenopus laevis to study retinoic acid-mediated mechanisms of orofacial development.
Through this work, we confirmed the efficacy and usefulness of Xenopus as a craniofacial model. Further, we determined signaling pathways and networks affected by RAR inhibition and provided functional validation of members of those networks in orofacial development. We identified several genes that are altered with RAR inhibition in Xenopus that correspond to genes altered in human disorders. This work provides evidence for how RA regulates orofacial development and identifies several targets for further research into the etiology of facial development in Xenopus, which may provide insight for human facial development.

Obtaining Xenopus embryos
Xenopus laevis adults were created in our breeding colony and purchased from Nasco. All procedures were approved by the VCU Institutional Animal Care and Use Committee (IACUC protocol number 5 AD20261). Embryos were collected using standard procedures [93] and were staged according to Nieuwkoop  ) were used to make two cuts to isolate the head: first-at the posterior end of the gut and then second caudal to the gills. Isolated heads were mounted in small depressions in a clay-lined dish containing Phosphate Buffered Saline with 0.1% Tween (PBT) [95]. The faces were imaged using a Discovery.V8 stereoscope fitted with an Axiovision digital camera (Zeiss). AxioVision40 V.4.8.1.0 software (Zeiss) was used to measure the intercanthal distance, the distance between the eyes. Statistical Analysis: Student's t-tests assuming unequal variance between each treatment group and corresponding control were performed to determine statistical relationships. Error bars representing standard error were calculated in Excel. Morphometrics analysis was performed as described previously with some modifications [95]. A total of 33 landmarks were placed on tadpole faces using the Pointpicker plugin of ImageJ (NIH, Additional file 1: Figure S1B). Landmarks 1-5 were placed on the midline of the face between the top of the cement gland and the top of the eyes. Landmarks 6-17 outlined each eye.  completed the outline of the midface region from the center of the eyes to the top of the cement gland, eliminating variation due to forebrain or cement gland size. Landmarks 28-33 outlined the mouth opening. Landmark coordinates were analyzed in MorphoJ 1.05f with Procrustes superimposition to remove information concerning scale, size, position, and orientation [96]. Statistical Analysis: To study statistical relationships among the treatment paradigms and corresponding controls, a canonical variate analysis (CVA) of Procrustes fit landmark data was performed and visualized as a bivariate plot. Discriminate function analysis (DFA) was also performed to examine statistical relationships between each treatment and its respective control. These were represented as transformation grids, in which the flat end of the vector represented average landmark position in control embryos, and the round end of the vector represented average landmark position in RAR-inhibitor exposed embryos. Statistical relationships are represented by p-values and all p-values that are less than 0.0001 are shown as p < 0.0001 as designated by MorphoJ.

Orofacial microdissections and RNA extraction
The orofacial region of control and RAR-inhibitor treated embryos were isolated by manual dissection at 35 or 66 hpf. Briefly, the embryo was placed in a chilled dish of 0.1X MBS (15C) and the head was isolated. The target area surrounding the stomodeum or mouth (Additional file 1: Figure S1A), excluding most of the brain and eyes was removed using a disposable scalpel with a 15 T blade Sklar Instruments, Cat. No 06-3120). Dissected tissues were added to TriReagent (Sigma, Cat. No. T9424) and frozen at − 80°C until RNA Isolation. 1 Biological replicate consisted of approximately 100 embryos from 2 to 3 pooled experiments. RNA extraction was conducted using the TriReagent protocol and lithium chloride precipitation (Ambion, AM9480). RNA quantity was assessed by Quawell UV-Vis spectrophotometer Q5000 and accompanying software.

Microarray and RNAseq
A microarray expression analysis was performed by the genome core at Virginia Commonwealth University (VCU) using the Xenopus laevis Genome 2.0 Array Gen-eChip containing probes to over 32,000 probe sets which represented more than 29,000 X. laevis transcripts (# 901214, Affymetrix). The VCU core performed RNA quality analysis, cDNA construction, and hybridization to the GeneChips, expression and statistical analysis. Statistical Analysis: Statistics and relative gene expression was performed with the Affymetrix statistical expression algorithm provided in the gene Chip Operating Software package. Data was compiled into an Excel file for further analysis in the Dickinson lab. X laevis genes that were significantly altered by 1.75 fold or greater (415 unique genes) were converted to the human ortholog for further functional and network analyses. Genetic similarity was determined to be > 90% identity and similarity. The files generated by this analysis are publically available at NCBI Gene Expression Omnibus [97] accession number GSE116819. Genes identified and used in this analysis are provided in table S2 (Additional file 4).
RNAseq was performed by Hudson Alpha Institute for Biotechnology (Huntsville, AL) where a HiSeq v4 PE100 250 mil reads 50 GB max output with Teir 1 analysis performed. The resultant data was analyzed by Dr. Stephen Turner and colleagues at the Bioinformatics Core at the University of Virginia (Charlottesville, VA). Reads were aligned using Xenopus laevis reference genome and STAR alignment software, and then counted using featureCounts software in the Subread package (http://www.xenbase.org/, RRID:SCR_003280) [98][99][100]. Statistical Analysis: The DESeq2 Bioconductor package in the R statistical computing environment to normalize count data, estimate dispersion, and fit a negative binomial model for each gene [101][102][103]. R Code available upon request. Data was then exported to Excel for further analysis in the Dickinson lab. Xenopus genes that were significantly altered (p < 0.05) were converted to the human ortholog for functional and network analyses (93 unique genes). Genetic similarity was determined to be > 90% identity and similarity. The files generated by this analysis are publically available at NCBI Gene Expression Omnibus [97], accession number GSE116819.

RT-PCR validation
RNA extraction was performed as described above. 1 μg of RNA was used to synthesize cDNA using the Qiagen quantitect reverse transcription kit (Qiagen 205,311) or the Biosystems high capacity cDNA reverse transcription kit (Thermofisher 4,368,814). RT-PCR was performed with SensiFAST SYBR No Rox (Bioline, BIO-98002) on a Bio-Rad CFX96 real time PCR system. Primer sequences are available upon request. Delta-delta CT method was used to calculate changes in expression, relative to expression levels of actin, in the PCR products.

DAVID and ingenuity network analysis
The Human Genome Organization (HUGO) Human Gene Nomenclature Committee (HGNC) database was utilized to convert Xenopus probes and reads to human nomenclature (HUGO Gene Nomenclature Committee at the European Bioinformatics Institute, genenames.org). The converted Xenopus gene names were used for all analyses and used within the text to maintain consistency and readability.
The significantly altered genes from each dataset were imported in DAVID Bioinformatics Resources [23,104]. Genes that were significantly altered were assessed for their molecular function using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database resources [105] and Gene Ontology (GO) Molecular Function from the Gene Ontology Consortium [106]. These functions were compiled and sorted according to genetic alteration (increase vs. decrease, relative to control). To further examine the functions of these genes, they were imported into Ingenuity Pathway Analysis software and the resultant networks examined (IPA version 31813283, Qiagen, Hilden, Germany). Custom pathways were built using the My Pathway tool.
In the microarray analysis, 415 unique genes with human orthologs were altered by 1.75 fold or greater with treatment 1 RAR inhibition. In the RNAseq Analysis, 93 unique genes with human orthologs were altered significantly with RAR inhibition during the late treatment window. For each network generated in IPA, the maximum number of altered genes included in the network was set at 30 genes. This number was chosen to ensure interpretability of the network. Occasionally genes that were altered by less than 1.75 fold in treatment 1 were included in the network if their function was integral to that network.

OMIM analysis
A keyword search was performed on the Online Mendelian Inheritance of Man (OMIM.org) to identify genes that were related to a craniofacial defect. Our search was limited to the keywords shown in Fig. 7. Three hundred and twenty five genes were identified, associated with 307 unique disorders (see Additional file 4: Table S2). These genes were compared to the genes that were significantly altered with early RAR inhibition (415 genes) and late RAR inhibition (93 genes) to identify the human disorders associated with genes altered with RAR inhibition.