Abstract

   Proper retinoic acid (RA) signaling is essential for normal
   craniofacial development. Both excessive RA and RA deficiency in early
   embryonic stage may lead to a variety of craniofacial malformations,
   for example, cleft palate, which have been investigated extensively.
   Dysregulated Wnt and Shh signaling were shown to underlie the
   pathogenesis of RA-induced craniofacial defects. In our present study,
   we showed a spatiotemporal-specific effect of RA signaling in
   regulating early development of facial prominences. Although inhibited
   Wnt activities was observed in E12.5/E13.5 mouse palatal shelves, early
   exposure of excessive RA induced Wnt signaling and Wnt-related gene
   expression in E11.5/E12.5 mouse embryonic frontonasal/maxillary
   processes. A conserved regulatory network of miR-484-Fzd5 was
   identified to play critical roles in RA-regulated craniofacial
   development using RNA-seq. In addition, subsequent
   osteogenic/chondrogenic differentiation were differentially regulated
   in discrete mouse embryonic facial prominences in response to early RA
   induction, demonstrated using both in vitro and in vivo analyses.

Introduction

   Retinoic acid, the metabolite of vitamin A, is the natural ligand of
   nuclear RA (retinoic acid) receptors and regulates the transcription of
   genes containing the RA response element (RARE). By recruiting either
   nuclear receptor coactivators or nuclear receptor corepressors, RA-RARs
   can directly activate/repress the transcription of key developmental
   genes such as Hox genes in hindbrain formation, Sox genes in
   neuromesodermal progenitor differentiation, and Fgf8 in developing
   heart ([34]Ghyselinck & Duester, 2019). Despite the discrepancies, both
   gain-of-function studies and genetic loss-of-function studies indicate
   that the dysregulated RA signaling has teratogenic effects and result
   in aberrant development of hindbrain, body axis formation, heart
   patterning, and eye morphogenesis ([35]Berenguer et al, 2022).

   Expression of retinaldehyde dehydrogenases, which is responsible for RA
   synthesis, starts at E7.5 mouse embryo in presomite mesoderm. The first
   embryonic RA signaling was detected and restricted within the
   developing eye around E8.5 ([36]Dubey et al, 2018). Excessive RA
   signaling during early embryonic development were shown to affect the
   fate of cranial neural crest cells (CNCCs) and resulted in craniofacial
   malformations including the cleft lip/palate ([37]Williams & Bohnsack,
   2019). Formation of the facial prominences appears around E8.5, upon
   receiving a great contribution of neural crest cells (NCCs) delaminated
   from between the dorsal neural tube and the overlying ectoderm
   ([38]Roth et al, 2021). The emigration and proliferation of NCCs, along
   with the paraxial mesoderm–originated mesenchymal cells, ensures the
   rapid growth of these facial prominences ([39]Murillo-Rincón & Kaucka,
   2020). By E11.5, the bilateral medial nasal prominences, lateral nasal
   prominences, and maxillary prominences are in close contacts and
   complete the fusion processes on each side of the developing face
   ([40]Ji et al, 2020). Ongoing convergence of these opposite structures
   in the midline form a continuous sheath of the front face and the
   primary palate ([41]Suzuki et al, 2018). Meantime, the mediolateral
   outgrowths of the maxillary form palatine processes starting around
   E11.5, growing vertically along the developing tongue in the mouth
   cavity ([42]Suzuki et al, 2018). The growth and elevation of the
   palatine processes lead to their horizontal fusion on top of the tongue
   in the midline and fusion to the primary palate just behind the
   incisive foramen by the time of E15.5 ([43]Martinelli et al, 2020). The
   mesenchyme derived from NCCs form ganglia, skeleton, and connective
   tissues, with musculature derived from paraxial mesoderm ([44]Roth et
   al, 2021).

   Underlying cellular events including cell proliferation, migration,
   apoptosis, and epithelial–mesenchymal transition of epithelium after
   fusion are well-synchronized processes, coordinated by key-signaling
   pathways, governing early morphogenesis. For instance, Wnt acts
   upstream of Bmp4 and induces its expression in overlapping domains
   ([45]Reynolds et al, 2020). The activation of Wnt and Bmp signaling is
   observed in the pre-fusion epithelia and the underlying mesenchyme of
   facial prominence to ensure proper mesenchymal growth, apoptosis, and
   epithelial fusion ([46]Jiang et al, 2006). However, excessive
   activation of Bmp signaling also causes craniofacial malformation.
   Knockout of noggin (Nog), an antagonist of Bmp pathway, displayed
   increased Bmp/Smad activity in E13.5 palate epithelium and resulted in
   a complete cleft palate in mice ([47]Reynolds et al, 2020). Although
   excessive intake of vitamin A and exogenous RA were shown to induce
   craniofacial malformation in humans and mice, the molecular
   pathogenesis associated with altered RA signaling was poorly understood
   ([48]Lammer et al, 1985; [49]Rothman et al, 1995; [50]Ackermans et al,
   2011). Inhibitory effects of RA on Tgf-β and Wnt signaling, whereas
   role of RA in maintaining Shh and Fgf8, were demonstrated in various
   cultured cell lines and in RA-induced CL/P in mouse and chick models
   ([51]Roa et al, 2019; [52]Reynolds et al, 2020; [53]Wu et al, 2022).
   However, controversial results showing RA-induced down-regulation of
   Shh expression were also presented ([54]El Shahawy et al, 2019;
   [55]Wang et al, 2019). The differences in animal models used and in the
   strategies of RA application may partially explain the discrepancies.

   In this work, we investigated the excessive RA-induced changes in gene
   expression profile in the developing craniofacial prominences using
   mouse model. Multiple dosages of all-trans RA (atRA) were applied to
   pregnant mice between early gestation (from E8.5 to E10.5), which was
   reported to induce full penetrance of CP ([56]Wang et al, 2019). The
   frontonasal prominence (FnP), a pair of maxillary processes (MxP) and a
   pair of mandibular processes (MdP) at E12.5 were collected and
   subjected to RNA-sequencing. The enrichment analysis of transcriptomic
   changes indicated the altered gene expression of 12 Wnt signaling
   components. A significantly down-regulated miRNA, miR-484, targeting
   five out of these 12 Wnt molecules, was predicted as a hub driver.
   Further screening for differentially expressed genes (DEGs) related to
   craniofacial malformations that were transcriptionally regulated by Wnt
   signaling revealed Foxn3 and Itgb1, a member of fork head family of
   transcription factors that regulate osteogenic gene expression and
   integrin β1 that regulate osteoblast differentiation. The gene
   expression of our bioinformatic analyses and transcriptome results were
   experimentally validated in prominence-specific fashion using
   quantitative qRT–PCR. The RA-induced different patterns of change in
   Wnt activities were assessed in craniofacial mesenchymal cells derived
   from various facial prominences at different stages using
   TCF/LEF-reporter TOP/FOP-flash assay. The subsequent influence of RA
   exposure on osteogenic differentiation was investigated both in vitro
   using mouse embryonic craniofacial mesenchymal cells derived from
   different facial prominences and in vivo using RA-induced mouse
   embryos. The results suggested a spatiotemporal-specific effect of
   early RA signaling in regulating Wnt activity and
   osteogenic/chondrogenic differentiation in discrete embryonic facial
   prominences.

Results

Early prenatal RA exposure induce developmental defects including midfacial
anomalies

   Prenatal RA exposure affects embryonic development in a both
   dose-dependent and embryonic stage–dependent manner. Oral
   administration of 25 μg/g body weight atRA consecutively at E8.5, E9.5,
   and E10.5 gestational timepoints was performed according to the
   description from [57]Wang et al (2019). Full penetrance of CP was
   observed, in line with the observation reported by [58]Wang et al
   (2019). Consistent with previous studies ([59]Berenguer et al, 2018),
   gross morphological analysis between E12.5–E15.5 revealed a range of
   malformations, resulting from defective early RA signaling, including
   microphthalmos, limb reduction, and superficial hemorrhages, with
   various incidences ([60]Fig 1A–G′).

Figure 1. Abnormalities at different developmental stages induced by early RA
exposure.

   [61]Figure 1.
   [62]Open in a new tab

   Representative images of some of the malformations resulted from early
   RA exposure. (A, A′, B, B′, D, D′, E, E′, F, F′, G, G′, H, H′)
   Whole-mount nuclear fluorescent imaging of control (A, B, D, E, F, G,
   H) and RA-treated (A′, B′, D′, E′, F′, G′, H′) embryos between E12.5
   and E15.5, displaying hypoplasia of frontonasal prominence (FnP) (A′),
   reduced growth of medial nasal prominence (B′), cleft of upper lip
   (D′), cleft palate (E′), and cleft in mandible (F′). Scale bars: 1 mm;
   PS, palatal shelf. (G, G′) Gross morphology of E15.5 control (G) and
   RA-treated (G′) mouse embryos. * indicate microphthalmos and
   superficial hemorrhages. Scale bars: 1 mm. (H, H′) H&E staining the
   coronal sections through the palatal regions of E15.5 embryos, showing
   failed fusion of PS in RA-treated mice, indicated with an arrowhead.
   Scale bars: 200 μm; NS, nasal septum; PS, palatal shelf; T, tongue. (C)
   Schematic of E12.5 mouse embryonic head. Dotted line illustrated the
   FnP, bilateral maxillary processes and mandibular processes that were
   dissected for transcriptomic profiling. FnP, frontonasal prominence;
   MxP, maxillary processes; MdP, mandibular processes.

   For transcriptomic profiling, the FnP, bilateral MxP and MdP were
   micro-dissected from fetuses at E12.5, at which stage the craniofacial
   abnormalities induced by RA could be easily verified ([63]Fig 1C).
   E12.5 is also the turning point from the completion of facial
   development to the beginning of palate formation, with the vertical
   palatal shelves (PS) being easily distinguished on both sides of the
   MxP. To minimize the chance of technical variance in sampling, each
   fetus was collected from an independent conception with (n = 4) or
   without RA exposure (n = 4). The dysregulated craniofacial development
   was confirmed by morphological analysis before sampling.

An overview of the transcriptome and quality assessment

   94.9% of the mRNA-seq datasets and 98.7% miRNA-seq datasets,
   respectively, remained after the initial filtering. The medians of
   log-expression (log-count-per-million) of each sample were presented
   post filtering using threshold of count-per-million >0.8 ([64]Fig 2A
   and B, Table S1). Genes were included in the analysis if their
   expression were identified in at least four samples with reads count
   >0, detected by filterByExpr function of edgeR (version 3.14.0). In the
   mRNA expression datasets, >97% of the paired-end raw reads were mapped
   using the mm (10) reference genome (GRCm38.p4). Of the 34,754 total
   genes, 17,024 (49.0%) gene expressions were detected in all samples. In
   the miRNA expression datasets, >80% of the paired-end raw reads were
   mapped using the reference miRBase (version 20.0). Of the 1,907 total
   genes, 835 (43.8%) gene expressions were detected in all samples. The
   expression of Zic3 was below the threshold in the mRNA expression
   datasets of all samples, excluding the potential contamination from
   adjacent forebrain tissue ([65]Inoue et al, 2007).

Figure 2. Early RA exposure induced altered gene expression in developing
mouse face.

   [66]Figure 2.
   [67]Open in a new tab

   (A, B) Boxplots of log-count-per-million values showing expression
   distributions for normalized data of mRNA-seq and miRNA-seq. (C, D)
   Principal component analysis of the facial gene expression of control
   (grey) and RA-treated (orange) E12.5 mouse embryos. (E, F) Volcano plot
   of the differentially expressed genes and differentially expressed
   miRNAs between control and RA-treated mouse embryos. The horizontal
   dashed line represents a P-value of 0.05, and the vertical dashed lines
   represent fold changes of −1.5 and 1.5, respectively. (G, H) Heatmap of
   hierarchical clustering dendrograms of top 50 differentially expressed
   genes and differentially expressed miRNAs in control and RA-treated
   mouse embryos. Control samples. Both mRNA and miRNA of two groups were
   clustered with approximately unbiased (AU) P-values and bootstrap
   probabilities more than 95%. Color corresponds to the relative
   expression levels, red for up-regulated and blue for down-regulated
   genes, respectively.

   Table S1. [68]Filtered gene expression values.^ (1.6MB, xls)

   The RA induced divergence in overall gene expression was demonstrated
   by the clustering of gene expression data among treatments in principal
   component analysis ([69]Fig 2C and D). We identified a total of 1,047
   DEGs (525 up-regulated and 522 down-regulated) and a total of 76
   differentially expressed miRNAs (DEMs) (54 up-regulated and 22
   down-regulated) with statistical significance (one-way ANOVA, adjusted
   P-value <0.05; [70]Fig 2E and F, Tables S2, and S3). The agglomerative
   hierarchical clustering of the DEG and DEM using Pearson correlation
   coefficient and average distance also supported the dominant
   differences between treatments. Cluster analysis with approximately
   unbiased (AU) P-values and bootstrap probability (BP) values showed
   that two main groups were identified, significant at least at the 95%
   CI (P <0.05, [71]Fig 2G and H).

   Table S2. [72]Differentially expressed genes.^ (192KB, xls)

   Table S3. [73]Differentially expressed miRNAs.^ (71.5KB, xls)

   The RA-synthesizing enzymes, Rdh10, Aldh1a2, and Aldh1a3 were
   abundantly expressed, and Cyp26b1 turned out to be the dominant subtype
   of enzyme for RA degradation, in agreement with the expression pattern
   of these genes in previous transcriptomic analysis of mouse embryonic
   craniofacial tissues ([74]GSE62214
   and [75]GSE7759) ([76]Fig S1). In addition, our results, along with the
   others, showed that mouse embryonic face expressed all isoforms of RA
   receptors including nuclear RA receptors (Rarα, Rarβ, Rarγ) and
   retinoid X receptors (Rxrα, Rxrβ, Rxrγ) ([77]Fig S1). Rxrγ was
   expressed at a lower level in the epithelium and mesenchyme of all
   facial prominences between E10.5 and E12.5, although relatively
   abundant Rxrγ was found in MdP ([78]GSE7759), primarily in the
   mandibular mesenchyme ([79]GSE62214) ([80]Feng et al, 2009; [81]Hooper
   et al, 2017). No significant differences in the expression of these
   genes were noticed between the control mouse embryos and RA-induced
   mouse embryos ([82]Fig S1A).

Figure S1. The expression pattern of RA metabolism–related genes in
transcriptomic analysis of mouse embryonic craniofacial tissues.

   [83]Figure S1.
   [84]Open in a new tab

   (A) The expression of RA metabolism–related genes between the control
   mouse embryos and RA-induced mouse embryos at E12.5. (B, C, D, E, F, G,
   H, I, J, K) The expression pattern of RA metabolism–related genes in
   previous transcriptomic analysis of mouse embryonic craniofacial
   tissues, in [85]GSE62214
   (B, C, D, E, F) and [86]GSE7759 (G, H, I, J, K), respectively.

Functional annotations and enrichment analysis revealed anomalous Wnt-related
gene expression in RA-treated embryonic craniofacial prominences

   Given the heterogeneity of the samples with respect to the tissue
   compositions and sex differences (with three females and one male in
   the control group whereas one female and three males in the RA-induced
   group), small fold changes might still be biologically important.
   Therefore, all DEGs with statistical significance were functionally
   annotated and undertaken enrichment analysis with no fold change
   cutoffs. Considering our particular interest in the early embryonic
   signaling, the Gene Ontology biological process terms (GO-BP terms) and
   Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enriched in
   genes differentially expressed were analyzed in up-regulated and
   down-regulated genes induced by RA exposure, respectively. The enriched
   GO-BP terms were identified with DAVID functional annotation (DAVID
   v2022q2) which include DAVID Gene Functional Classification Tool and
   DAVID Functional Annotation Clustering Tool. The resulting enriched
   terms were ranked by gene count and the Wnt-signaling–associated terms
   were the only specific developmental pathway recognized in up-regulated
   DEGs with adjusted P <0.05 ([87]Fig 3A, Table S4). No over-represented
   categories of GO-BP terms was identified in the down-regulated DEGs
   ([88]Fig S2A). Furthermore, KEGG pathway analysis also recognized the
   over-represented genes involving Wnt signaling in the up-regulated
   DEGs, with 2.79-fold enrichment at adjusted P <0.05 ([89]Fig 3B, Table
   S5). Down-regulated DEGs lacked enrichment in any categories of
   specific developmental signaling pathways ([90]Fig S2B).

Figure 3. The top 20 Gene Ontology biological process terms and top 20 Kyoto
Encyclopedia of Genes and Genomes pathways in enrichment analysis of
up-regulated differentially expressed genes.

   [91]Figure 3.
   [92]Open in a new tab

   (A, B) The colors of dots (A) and bars (B) indicate the adjusted
   P-value of the Fisher’s exact test. The dot size indicates the number
   of genes, and the X-axis represents the percentage of differentially
   expressed genes identified in each term.

Figure S2. Enriched functional annotation and pathway analysis of
down-regulated differentially expressed genes.

   [93]Figure S2.
   [94]Open in a new tab

   Table S4. [95]Gene Ontology biological process terms enriched in genes
   differentially expressed by prominence.^ (70.5KB, doc)

   Table S5. [96]Kyoto Encyclopedia of Genes and Genomes pathway
   enrichment analysis in genes differentially expressed by prominence.^
   (79KB, doc)

Identification of mouse-specific miRNA–gene interactions targeting the
up-regulated genes in Wnt pathway in DEMs

   First, we searched TargetScan 8.0 and miRWalk for potential miRNAs–gene
   interactions for the 12 up-regulated genes enriched in Wnt signaling in
   KEGG pathway analysis. miRNAs were only considered as potential
   regulators of these genes when predicted in both databases, and as a
   result, 1,584 miRNAs and 5,473 miRNA–gene interactions were initially
   identified (Table S6). Next, the overlapping miRNAs between these
   potential miRNA regulators and the RA-induced down-regulated DEMs were
   selected ([97]Fig S3A). 26 nodes (14 DEMs and 12 Wnt genes) and 38
   edges with 47 MREs in the mRNA of these 12 Wnt-related genes were
   identified ([98]Fig 4A, Table S7). According to the maximal clique
   centrality (MCC) score calculated using cytoHubba plugin (version 0.1),
   miR-484 (MCC = 5) were identified as a hub miRNA. Four miR-484 binding
   sites on Fzd5 3′UTR were predicted by TargetScan ([99]Fig S3B) and the
   miR-484: Fzd5 mRNA association was further characterized using
   luciferase assay. 515 bp murine Fzd5 3′UTR containing either three
   clustered WT seed sequences of miR-484 or sense mutations of all these
   three seed sequence was cloned into Sac I/Xba I site of pmirGLO vector,
   respectively ([100]Fig S3C). Reduced luciferase activities were
   observed in primary mouse embryonic mesenchymal cells and in HEK-293T
   cells co-transfected with pmirGLO-WT-Fzd5 3′UTR and miR-484 mimics,
   whereas no significant change was observed in cells co-transfected with
   pmirGLO-MUT-Fzd5 3′ UTR and miR-484 ([101]Figs 4B and [102]S3D).
   Scramble miRNA transfections were used as controls.

Figure S3. miRNAs in down-regulated differentially expressed miRNAs with
potential targets of the up-regulated Wnt-related genes in differentially
expressed genes.

   [103]Figure S3.
   [104]Open in a new tab

   (A) Venn diagram showed the overlapping of down-regulated
   differentially expressed miRNAs and the potential miRNA regulators of
   the up-regulated 12 Wnt genes predicted by TargetScan and miRWalk. (B)
   miR-484–binding sites on 3′UTR of Fzd5 predicted by TargetScan. (C) The
   515 bp Fzd5 3′UTR containing mouse WT or sense mutation in the three
   clustered seed sequences of the miR-484–binding sites. (D) Luciferase
   assay showed the direct association of miR-484 and the 3′UTR of Fzd5 in
   cultured 293T cells. The data were presented as means ± SD (n = 3).
   One-way ANOVA, Turkey’s test.

Figure 4. miR-484 was recognized as the hub regulator.

   [105]Figure 4.
   [106]Open in a new tab

   (A) miRNA–mRNA regulatory network of the down-regulated differentially
   expressed miRNAs and the up-regulated Wnt-related genes induced by RA.
   The size of miRNA nodes indicated the number of their target genes. (B)
   Fzd5 was verified to be the direct target miR-484 using luciferase
   assay. Luciferase reporter construct contained three clustered WT seed
   sequences of miR-484 or sense mutations of all three seed sequences
   were co-transfected with miR-484 mimics or scramble RNA into primary
   mouse mesenchymal cells. Luciferase activities were measured 48 h post
   transfection. The relative luciferase activities were normalized to
   blank controls of mock transfection. The data were presented as means ±
   SD (n = 3). **P < 0.01, one-way ANOVA, Turkey’s test. (C) The relative
   expression of miR-484 in each facial prominences of E12.5 embryos of
   RA-induced mice and controls. (D) The relative expression of miR-484 in
   cultured primary embryonic mesenchymal cells derived from E12.5
   frontonasal prominence/maxillary processes/mandibular processes/PS,
   induced with or without 1 μM atRA. (E) The sequence of the miR-484
   promoter containing WT or mutant retinoic acid response element, boxed
   and highlighted. The transcriptional start site was indicated. (F) Fold
   change of the expression of luciferase driven by the WT and the mutated
   promoter of miR-484 in response to 1 μM atRA induction. Values (C, D,
   F) were presented in means ± SD (n = 3). *P < 0.05, **P < 0.01, ***P <
   0.001, two-tailed, unpaired t test.

   Table S6. [107]miRNA–gene interactions.^ (189KB, xls)

   Table S7. [108]miRNA–gene network.^ (63KB, xls)

   The RA-induced change of expression of miR-484 was subsequently
   verified in vivo using qRT-PCR. Interestingly, distinctive regulation
   of miR-484 expression was observed in different facial prominences,
   with RA-induced down-regulation in E12.5 FnP/MxP/MdP and RA-induced
   up-regulation in E12.5 palate shelves ([109]Fig 4C). We further
   verified the prominence-specific regulation of miR-484 in vitro using
   cultured primary mesenchymal cells collected from each individual
   prominence. Reduced miR-484 expression in mesenchymal cells derived
   from E12.5 FnP/MxP/MdP whereas elevated expression of miR-484 in
   mesenchymal cells derived from E12.5 palate shelved were observed upon
   1 μM atRA induction ([110]Fig 4D). 1,442 bp sequence of the
   proximal-promoter region from −1,437 to +5 upstream with respect to the
   transcription start site of miR-484 gene was used to identify the
   potential presence of RARE. Two putative RAR-binding motifs were found
   at −40 to −29 and −190 to −179 ([111]Fig 4E). A 500 bp
   promoter–enhancer region containing the WT or mutant sequence of RAREs
   was cloned into the vector pEZX-FR01, which were transiently expressed
   in HEK-293T cells and then induced with or without 1 μM atRA,
   respectively. The RA-induced fold change of the relative expression of
   luciferase driven by promoter containing mutated proximal RA
   cis-regulatory element (mut#1) was significantly reduced, compared with
   that driven by the WT promoter ([112]Fig 4F). Mutation of both RAREs
   (mut#2) in the promoter resulted in further decreases in the RA-induced
   luciferase reporter gene expression. The results confirmed the
   predicted cis-acting motif in the proximal-promoter region of miR-484
   gene per se was RA responsive.

RA enhanced Wnt signaling and induced Wnt-related gene expression during
early craniofacial development

   The Wnt activities were determined in E12.5 mouse primary craniofacial
   mesenchymal cells and in E12.5/E13.5 mouse embryonic palatal
   mesenchymal (MEPM) cells using Wnt reporter TOP/FOP-flash. The
   mesenchymal cells were isolated from E12.5 MxP, E12.5 MdP, E12.5, and
   E13.5 PS, respectively, from the embryos of control and RA-treated
   pregnant mice. The MEPM cells collected from E13.5 mouse PS were used
   as references for palatal development that were widely used by others.