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

Abstract

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.

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.

Results

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).

Fig. 1

RAR Inhibition during early and late orofacial development. a Schematic of the RAR inhibitor treatment paradigm. b-g Representative images of frontal views of stage 43 (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). c,d Two different embryos treated from stage 24 to stage 29/30 (26–35 hpf) with 1 uM RAR Inhibitor. Black arrow indicates a triangular shaped mouth e Embryos treated with 1% DMSO from stage 29/30 to stage 40 (35–66 hpf). f,g Two different embryos treated with 1 uM RAR Inhibitor from stage 29/30 to stage 40 (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

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 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.

Table 1 Most dysregulated genes after early RAR inhibition. The top 20 decreased (left) and top 20 increased (right) genes after RAR inhibition during the early treatment window

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.

Table 2 Most dysregulated genes after late RAR inhibition. The top 20 decreased (left) and top 20 increased (right) genes after RAR inhibition during the early treatment window

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)

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).

Fig. 2

Functional classification of altered genes with RAR inhibition. (a). Functional categories of genes decreased (i) or increased (ii) after early RAR inhibition. (b). Functional categories of genes decreased (i) or increased (ii) after late RAR inhibition. For all graphs: numbers represent the percent of genes identified. (ci-ii): 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

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 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 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.

Fig. 3

Transcription regulation was altered in early RAR inhibition. a Functional network built in IPA software, utilizing DAVID pathway analysis. Blue genes are decreased relative to control; orange genes are increased relative to control. Legend denotes the type of gene product represented in the network. b-c Representative images of Control (b) and CHD1 morphants (c). Facial structures are labeled in b. d-e Representative images of Control (1% DMSO, d) and TSA treated embryos (e). In all representative images the mouth is outlined in red dots

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.

Fig. 4

Transcription regulation was largely decreased in late RAR inhibition. Functional network build in IPA software, utilizing IPA and DAVID pathway analysis. Green genes are decreased relative to control; purple genes are increased relative to control. Key denotes the type of gene product represented in the network

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)

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).

Fig. 5

Calcium signaling in Organ Development and Differentiation were altered with early RAR inhibition. a Functional network built in IPA software utilizing IPA and DAVID pathway analysis. Blue genes are decreased relative to control; orange genes are increased relative to control. Key denotes type of gene product represented in diagram. b-c Representative images of control (b, 1% DMSO) and CAMKII inhibitor treated embryos (c, 100 μM KN93). d-e Representative images of control (1% DMSO, d) or calcineurin inhibitor treated embryos (e, 100 μM FK506). Mouths are outlined in red dots

Table 5 Genes dysregulated after early RAR inhibition associated with organ development and differentiation. Genes shown in Fig. 5 along with their general function (adapted from genecards.org)

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)

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).

The antagonist FK506, a previously validated inhibitor, was used to inhibit Calcineurin during early orofacial development in Xenopus [34]. Embryos were treated with 100 μM of FK506 from stage 24–35 (26–50 hpf, n = 20, 2 biological replicates) and facial morphology was assessed at stage 43 (87–92 hpf). FK506 treated embryos had a narrower midface, protruding lenses, and a square shape mouth (Fig. 5d-e).

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 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.

Fig. 6

Wnt and BMP signaling were altered with late RAR inhibition. a Functional network built in IPA software utilizing IPA and DAVID functional analyses. Green genes are decreased and purple genes are increased relative to control. Key denotes type of gene product represented in diagram. b-c Representative images of control (b, 1% DMSO) and Wnt inhibitor treated embryos (c, IWR-1 10 μM). d-e Representative images of control (d, 1% DMSO) and BMPR Inhibitor (e, LDN193189, 100 μM). Mouths are outlined in red dots

The pharmacological inhibitor IWR-1, was next utilized to block the WNT-induced accumulation of β-catenin through stabilization of the destruction complex member AXIN2 [37,38,39,40]. Embryos were treated with 10 μM IWR-1 from stage 29/30–40 (Fig. 6b-c, 35–66 hpf, n = 20 from 3 biological replicates). At stage 43 (87–92 hpf) facial characteristics were assessed. IWR-1 treated embryos showed lens and eye defects, a narrower face and smaller misshapen mouth (Fig. 6c).

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 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.

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

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)

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. 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, SMARCAD1, 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²⁺/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 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 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 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.

Methods

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 and Faber [94]. Briefly, eggs were collected in high salt 1X Modified Barth’s Saline (MBS, pH 7.8) from at least two different females and in vitro fertilized using sperm obtained from one male. Embryos were cultured in 0.1X MBS, refreshed daily, until used for experiments. Embryos were housed in a 23 °C or 15 °C incubator (Torrey Pines Scientific, Cat. No. IN30). After the experiments were completed all embryos were given a lethal dose of anesthetic (10% tricaine).

Pharmacological inhibitor treatments and antisense Morpholinos

Two periods of development were chosen for our analysis based on craniofacial developmental events. An early period or phase, when the neural crest is migrating and facial prominences are specified. A later phase, when the face is growing and the jaw muscle and cartilage just begin to form. At the beginning of each phase (early = stage 24 or late = 30), embryos were arrayed in 60 mm culture dishes containing 10 ml of 0.1X MBS with the desired chemical and DMSO (final concentration of DMSO was 1%). Chemicals used for analysis were: 1 μM RAR antagonist (BMS453, 1 mM stock solution in DMSO, Tocris, Cat. No. 3409), 100 μM CaMKII antagonist (KN93, 5 mM stock solution in water, Sigma-Aldrich Cat. No. 422711-M), 125 μM HDAC inhibitor (Trichostatin A, 5 mM stock solution in DMSO, Sigma-Aldrich, Cat. No T8852), 25 μM Dyrk1a inhibitor (INDY, 100 mM stock solution in DMSO, Tocris, Cat. No. 4997), 40 μM Wnt1 inhibitor (10 mM stock solution in DMSO, Sigma-Aldrich, Cat. No. I0161), and 25 μM BMPR1 Inhibitor (LDN-193189, 10 mM stock solution dissolved in water, Calbiochem, Cat. No. 509882). During pharmacological treatments, embryos were maintained at 23 °C in an Echotherm incubator (Torrey Pines Scientific, Cat. No. IN30). Antisense morpholinos (MOs) with fluorescein tags were purchased from Genetools. 60 ng of a splice blocking Chd1 MO (CTAGTCACTTTACCTCTGCTCCATC) was injected into the 1 cell stage. The equivalent standard control morpholino was used as a control in this experiment. In all experiments at least 40 embryos were treated in 2 biological replicates.

Intercanthal distance measurements and Morphometrics

Embryos were fixed in 4% paraformaldehyde (PFA) overnight at 4 °C. A sterile disposable No. 15 scalpel (VWR, cat. no: 82029–856) and Dumont No. 5 forceps (Fisher, cat. no: NC9404145) 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. Landmarks 18–27 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 GeneChip 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.