Identification of microRNAs and gene regulatory networks in cleft lip common in humans and mice.

The etiology of cleft lip with/without cleft palate (CL/P), one of the most frequent craniofacial birth defects worldwide, is complicated by contributions of both genetic and environmental factors. Understanding the etiology of these conditions is essential for developing preventive strategies. This study thus aims to identify regulatory networks of microRNAs (miRNAs), transcriptional factors (TFs) and non-TF genes associated with cleft lip (CL) that are conserved in humans and mice. Notably, we found that miR-27b, miR-133b, miR-205, miR-376b and miR-376c were involved in the regulation of CL-associated gene expression in both humans and mice. Among the candidate miRNAs, the overexpression of miR-27b, miR-133b and miR-205, but not miR-376b and miR-376c, significantly inhibited cell proliferation through suppression of CL-associated genes (miR-27b suppressed PAX9 and RARA; miR-133b suppressed FGFR1, PAX7, and SUMO1; and miR-205 suppressed PAX9 and RARA) in cultured human and mouse lip mesenchymal cells. Taken together, our results suggest that elevated expression of miR-27b, miR-133b and miR-205 may play a crucial role in CL through the suppression of genes associated with CL.

Introduction

Cleft lip with/without cleft palate (CL/P), one of the most common birth defects worldwide (1), constitutes a significant public health issue that negatively impacts esthetics, oral functions such as sucking, swallowing, and speaking, and quality of life (2). Approximately 70% of CL/P cases are non-syndromic, whereas the remaining 30% are syndromic (3,4). The etiology of non-syndromic CL/P is very complex, with contributions of both genetic and non-genetic factors, including environmental factors along with geographic, racial and ethnic influences, as well as maternal factors such as inadequate nutrition, smoking, and alcohol and drug consumption (3,4). Although novel gene mutations associated with CL/P are reported every year, it is unclear how these multiple genes are tightly regulated in a spatial and temporal manner.

The upper lip is formed through the growth and fusion of the maxillary process (MxP) and the medial nasal process (MNP) and lateral nasal processes (LNP) at the 4th–6th week of gestation in humans (5,6); any failure in these processes results in either unilateral or bilateral cleft lip (CL). Mouse models have been extensively used to study the pathological mechanisms of CL, since mouse lip formation is similar to that of humans, with a short time window, and the underlying molecular mechanism is well conserved. In mice, upper lip development begins at E9.5, when the stomodeum is bound laterally by a pair of the MxP, and caudally by a pair of mandibular processes (7–9). At E10.0, the nasal processes bulge within the ventrolateral region of the frontonasal prominence and outgrow around the periphery of the nasal placodes in a horseshoe shape, forming the nasal pit (6). At E10.5, the lateral growth of the MxP pushes the nasal pits toward the middle line of the face, and the ventral tip of the NP extends and fuses at the medial end of their boundary junction (6,10). By E11.5, the fusion of the nostrils and lip completes the disappearance of the medial edge epithelium while the MxP further grows and pushes the nostrils into the medial region (6,11). As in humans, any failure or delay in the outgrowth and fusion of these processes results in CL in mice.

Only ~2% of the human genome encodes genes that are transcribed and translated into proteins. The remaining 98% consists of so-called ‘noncoding sequences’, which include introns (26%), long terminal repeat retrotransposons (8%), short interspersed nuclear elements (13%), long interspersed nuclear elements (21%), segmental duplications (5%), simple sequence repeats (3%), DNA transposons (3%) and unique sequences (12). Although the role of non-coding DNA is not fully understood, it functions in various aspects of cellular physiology. Non-coding ribonucleic acids (RNAs), including transfer RNA, ribosomal RNA, microRNA, small nuclear RNA, circular RNA and long-non-coding RNA, are transcribed from either coding or non-coding regions and play essential roles in DNA transcription, protein translation and synthesis, and RNA processing (13–15). Introns, untranslated 5′ or 3′ regions, and regulatory DNA sequences such as the binding sites of transcription factors (TFs) or microRNAs are associated with the regulation of gene expression (16,17). Repetitive DNA sequences comprise about 50% of the human genome and are classified as tandem and interspersed repeats (18–21); for example, telomeres, which are regions of tandem repeats, protect chromosomes against degradation. The expansion of non-coding nucleotide repeats and trinucleotide repeats cause various diseases, including Huntington’s disease (22–24).

Recent studies indicate that environmental and epigenetic factors contribute to alternation between active/suppressed gene expression through the regulation of non-coding RNAs (25,26). miRNAs, which are small RNAs comprising 21–22 nucleotides, are typically transcribed as pri-miRNAs ranging from hundreds to thousands of nucleotides in length, and then processed into a stem-loop structure with approximately 70 nucleotides long in the nucleus by the microprocessor complex (MC) comprising DROSHA and DGCR8 (27,28). The pri-miRNAs are bound to EXPOTIN5, which exports pre-miRNAs into the cytoplasm, and are further processed into mature double-strand miRNAs by RNase III enzyme DICER. The mature miRNAs are then loaded into the ARGONAUTE protein guide-strand channel to form an RNA-induced silencing complex (RISC) (29,30). With a guide of a 5′ end of seed sequence of mature miRNAs, the miRNA–RISC complex binds their mRNA targets at the 3′ UTR, resulting in inhibition of translation or destabilization of mRNA (31). Thus, miRNAs negatively regulate the expression of target mRNAs at the post-transcriptional level.

An increasing body of evidence points to a crucial role for miRNAs in embryonic development and various physiological and pathological conditions (32), and recent studies show that a large number of miRNAs are expressed in the developing lip and palate (33–36) and contribute to the etiology of CL/P in humans (33,37,38). For example, a human genetic association study showed that single-nucleotide polymorphisms (SNPs) in DROSHA are associated with susceptibility to non-syndromic CL/P (39). A deficiency of Drosha in mice results in embryonic lethality by E7.5 (40), while a conditional null mutation in craniofacial tissues has not been reported yet. Similarly, Dgcr8-deficient mice die by E5, whereas heterozygous Dgcr8 mutant mice can survive until adulthood, although they display defective adult neurogenesis (41,42). As other examples, mice with a conditional deletion of Dgcr8 in cranial neural crest (CNC) cells, from which the majority of craniofacial mesenchymal cells derive (Wnt1-Cre;Dgcr8 mice), exhibit cardiac defects and severe craniofacial deformities, including agnathia and midfacial hypoplasia (43). Similarly, mice with a deficiency of Dicer1 in CNC cells (Wnt1-Cre;Dicer1 mice) exhibit severe craniofacial deformities, including severe midfacial hypoplasia, cleft palate of the secondary palate, lack of maxillary bony components and agnathia (44–46). In addition, Pax2-Cre;Dicer1 mice, in which Pax2 is expressed at E7.5 in the prospective midbrain-hindbrain region where CNC cells migrate into the first pharyngeal arch, exhibit severe craniofacial deformities, including midfacial hypoplasia, cleft palate of the secondary palate, lack of the palatal process of the maxilla and micrognathia (47). By contrast, mice with an epithelium-specific deletion of Dgcr8 (K14-Cre;Dgcr8 mice) display no craniofacial developmental defects but exhibit defects in the hair germ and epidermis (48). Moreover, mice with conditional deletions of Dicer1 in the epithelium, i.e. K14-Cre;Dicer1 and Shh-Cre;Dicer1, display no cleft in the lip and palate but present tooth developmental defects such as abnormal molar shape and enamel formation defects, and supernumerary upper incisors, respectively (49,50). Taken together, these observations point to a crucial role for miRNAs in CNC cells, but not in K14-positive epithelial cells, during lip and palate development. However, the specific miRNAs involved in lip formation and pathogenesis of CL, and how they regulate other genes, remain elusive.

TFs and miRNAs are the largest families of molecules that share a common regulatory pattern (51), with sets of coexpressed TFs and miRNAs precisely regulating gene expression in each cell type. By binding to discrete sequence elements, individual TFs (binding to promoter regions of genes, including miRNA genes) and miRNAs (binding to the 3′ UTR regions) can regulate multiple target genes and bind to their cognate sequences, either alone, cooperatively or even TF-miRNA composites (52,53); the cooperative regulatory mechanism through feed-forward loops (FFLs) and feedback loops (FBLs) has been reported in various diseases (54–56). Since the expression of a large number of genes is altered in CL, the etiology may be explained by the disruption of miRNA–gene regulatory networks commonly involved in both human and mouse CL.

Results

Prediction for regulatory pairs and loops of genes, TFs and miRNAs associated with CL

A total of 173 human genes and 131 mouse genes associated with CL were collected through systematic reviews (57–60) (Supplementary Material, Tables S1 and S2). The human CL-associated genes were grouped into 31 TF genes and 142 non-TF genes. In addition, 140 miRNAs regulating these CL-related genes were collected through a literature search (Fig. 1A). Using CL-associated gene sets common in humans and mice, we conducted bioinformatic analyses as shown in the overall analysis framework (Fig. 2). We estimated possible molecular interactions and found 241 miRNA–gene interactions, 102 miRNA–TF interactions, 579 TF–gene interactions and 330 TF–miRNA interactions (Fig. 1B and Table 1). In addition, we found that these interactions formed three types of motifs (Fig. 2): 214 type I motifs (TF regulates miRNA and gene; miRNA represses gene), 140 type II motifs (miRNA represses TF and gene; TF regulates gene) and 14 type III motifs (TF regulates miRNA and gene; miRNA represses TF and gene) (Table 2, Supplementary Material, Tables S3–S5). In mice, 131 genes were subcategorized into 37 TF genes, 94 non-TF genes and 141 miRNAs (Fig. 2A). Among these molecules, we found that there were 452 miRNA–gene interactions, 410 miRNA–TF interactions, 406 TF–gene interactions and 385 TF–miRNA interactions (Table 3). The number of types of motifs was 450 for type I, 328 for type II and 88 for type III (Table 4,  Supplementary Material, Tables S6–S8).

Figure 1

Flowchart of the miRNA–TF regulatory network analysis for CL in humans and mice. (A) Lists of miRNAs, TFs and non-TF target genes in humans and mice manually curated from the literature. (B) Different types of interactions between TFs, miRNAs and non-TF target genes. (C) Three types of regulatory motifs (Type I, Type II and Type III). (D) The integrated miRNA–TF–non-TF target gene regulatory network. (E) Hub gene selection, subnetwork identification and gene set enrichment analysis of network genes. (F) The consensus miRNAs/TFs/non-TF genes, subnetworks, motifs, pathways and GO terms shared between humans and mice.

Figure 2

Distribution and intersections of miRNAs, TFs and non-TF target genes, and their pairwise interactions between humans and mice. (A) Distribution and intersections of miRNAs, TFs and non-TF target genes between humans and mice. (B) Distribution and intersections of TF–non-TF gene, TF–miRNA, miRNA–non-TF gene and miRNA–TF between humans and mice. (C) Motifs with unique or consensus interactions in humans and mice.

Table 1

Summary of the integrated regulatory relationships among human CL-related miRNAs, TFs and non-TF target genes

Relationship	#pairs	#miRNAs	#genes	#TFs	Methods
miRNA–Gene	241	27	68	/	A
miRNA–TF	102	23	/	24	A
TF–Gene	579	/	132	23	B
TF–miRNA	330	119	/	20	B
Total	1252	126	134	28

sTwo methods are used. Method A: regulation relationships were supported by at least two of the following four databases: miRanda, miRTarBase, PITA and TargetScan. Method B: FIMO was used for the identification of TF–target–non-TF target gene relationships.

Table 2

Summary of human CL-associated 3-node feed-forward loops (FFLs) and feedback loops (FBLs)

Type		#nodes			#pairs
	Motifs	TFs	miRNAs	Genes	miRNA–Gene	miRNA–TF	TF–Gene	TF–miRNA
I	214	10	19	54	49	/	131	35
II	140	16	20	35	90	56	75	/
III	14	4	4	13	14	4	14	4
Total	368	18	24	54	160	56	176	49

Table 3

Summary of the integrated regulatory relationships among mouse CL-related miRNAs, TFs, and genes

Relationship	#pairs	#miRNAs	#genes	#TFs	Methods
miRNA–Gene	452	99	59	/	A
miRNA–TF	410	104	/	34	A
TF–Gene	406	/	84	24	B
TF–miRNA	385	111	/	22	B
Total	1653	134	94	36

sTwo methods are used. Method A: Regulation relationships were supported by at least two of the following four databases: miRanda, miRTarBase, PITA and TargetScan. Method B: FIMO was used for the identification of TF–target relationships.

Table 4

Summary of CL-associated mouse 3-node feed-forward loops (FFLs) and feedback loops (FBLs)

Type		#nodes			#pairs
	Motifs	TFs	miRNAs	Genes	miRNA–Gene	miRNA–TF	TF–Gene	TF–miRNA
I	450	15	71	45	247	/	145	176
II	328	19	74	41	231	155	101	/
III	88	10	25	20	77	33	49	33
Total	866	20	92	48	325	155	168	176

To determine conserved mechanisms in humans and mice, we extracted and analyzed the genes and miRNAs using miRbase (61) for miRNAs and the NCBI HomoloGene database (https://www.ncbi.nlm.nih.gov/homologene, Release 68) for genes (both TFs and non-TFs) (62). We identified nine miRNAs (miR-124, miR-133b, miR-205, miR-24, miR-27b, miR-369, miR-376a, miR-376b and miR-376c), eight TF genes (GLI2, MSX1, MSX2, PAX3, PAX7, PAX9, RARA and SATB2) and eight non-TF genes (BMP4, FGFR1, KIF7, PTCH1, SPRY1, SUMO1, WNT5A and WNT9B) as CL-associated molecules conserved in humans and mice (Figs 1 and 2A). A total of 11.0% (19/173) of human CL-associated genes and 14.5% (19/131) of mouse CL-associated genes were commonly identified, suggesting that these genes are more likely to be crucial for lip formation.

Next, we performed gene set enrichment analysis of the consensus TF and non-TF genes (n = 19) using the DAVID online tool (63) and Gene Ontology (GO) annotation terms. The significantly enriched GO terms included GO:0000122 ‘negative regulation of transcription from RNA polymerase II promoter’ [human, false discovery rate (FDR) = 1.79 × 10−3; mouse, FDR = 4.96 × 10−5], GO:0045892 ‘negative regulation of transcription, DNA-templated’ (human, FDR = 7.73 × 10−2; mouse, FDR = 2.53 × 10−4), GO:0045944 ‘positive regulation of transcription from RNA polymerase II promoter’ (human, FDR = 8.52 × 10−4; mouse, FDR = 2.44 × 10−5), and GO:0060021 ‘palate development’ (human, FDR = 7.39 × 10−4; mouse, FDR = 8.80 × 10−4) (Supplementary Material, Tables S9 and S10). Interestingly, among these enriched signaling pathways, all except for the Rap1 pathway have been reported in CL mouse models (64–68).

Next, we examined the consensus pairwise interactions within the tri-molecule motifs (e.g. TF–gene, TF–miRNA, miRNA–gene and miRNA–TF) between humans and mice and found two consensus TF–gene interactions, two consensuses TF–miRNA interactions, five consensus miRNA–gene interactions and nine consensus miRNA–TF interactions (Figure 2B,  Supplementary Material, Table S11). We compared the tri-molecule motifs that were common in humans and mice and identified one consensus motif (a miRNA–TF–target gene motif): miR-376b–GLI2–FGFR1 (Fig. 2C).

Regulatory networks and features of human CL-related miRNAs, TFs and genes

To construct the human-specific miRNA–TF–gene regulatory network, we integrated all three types of human motifs and found 96 nodes (24 miRNAs, 18 TFs and 54 genes) and 441 edges (Supplementary Material, Fig. S1). Using the maximal clique centrality (MCC) (69) score obtained through the cytoHubba plugin (version 0.1), we identified 10 hub nodes from this network, including five miRNAs (miR-27b, MCC = 15; miR-30b, MCC = 17; miR-374a, MCC = 19; miR-374b, MCC = 18; and miR-381, MCC = 15) and five TFs (ESRRG, MCC = 8; FOXE1, MCC = 38; RARA, MCC = 18; RUNX2, MCC = 13; and TFAP2A, MCC = 27). The target genes (Supplementary Material, Table S12) were subjected to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the WebGestalt online tool (70). The target genes of two hub TFs (FOXE1 and RUNX2) and two hub miRNAs (miR-374a and miR-381) were significantly enriched in several pathways closely related to the CL-specific biological processes (Fig. 3A). For example, the target genes (n = 38) of TF FOXE1, known to be associated with CL (71), were significantly enriched in the signaling pathways regulating pluripotency of stem cells (FDR = 1.38 × 10−2) and the Hippo signaling pathway (FDR = 2.13 × 10−2). The target genes of miR-374a were enriched in the signaling pathways regulating pluripotency of stem cells (FDR = 8.82 × 10−2) using an FDR cut-off value of 0.1.

Figure 3

KEGG pathway enrichment analysis. (A) KEGG pathway enrichment analysis of TF targets and miRNAs related to human CL. (B) KEGG pathway enrichment analysis of targets of TF and miRNAs related to mouse CL.

TFs and miRNAs crucial for the Wnt pathway in human CL

Previous studies show that the dysregulation of Wnt pathway plays an essential role in CL (72); for instance, hypermethylation of WNT3A is associated with risk for non-syndromic CL/P (73). The KEGG pathway database (74) identified that, among the human CL-associated genes, 11 genes were related to the Wnt pathway (AXIN2, DVL2, FZD6, WNT10A, WNT11, WNT3, WNT3A, WNT5A, WNT5B, WNT6 and WNT9B). Next, we constructed a subnetwork by integrating all motifs for the 11 genes identified and found 24 nodes (nine miRNAs, nine TFs and six non-TF genes) and 62 interactions involved in the subnetwork (Fig. 4A); for instance, FOXE1, GLI2 and TP63 were direct upstream regulatory molecules of DVL2, WNT3A, WNT5A and WNT5B.

Figure 4

miRNA–TF–non-TF target gene regulatory subnetwork in humans and mice. (A) Network consisting of miRNA–TF–non-TF target gene motifs related to the Wnt signaling pathway. (B) Extracted network modules from the miRNA–TF–non-TF target gene network in humans. (C) Network consisting of miRNA–TF–non-TF target gene motifs related to the Hedgehog signaling pathway. (D) Extracted network modules from the miRNA–TF–non-TF target gene network in mice.

Next, using the Markov cluster (MCL) algorithm (75) with the default setting (inflation = 2.5), we identified three network modules from the corresponding miRNA–TF–gene regulatory network (Fig. 4B). Each module comprised a set of miRNAs, TFs and target genes that were adjacent in terms of network topology. Among these modules, four TFs (E2F1, FOXE1, GLI2 and TP63) were considered to be the most promising CL nodes because the network interactions’ direction stated that these three regulatory molecules could modulate the others at the transcriptional or post-transcriptional level (Fig. 4A).

Regulatory networks and features of mouse CL-related miRNAs, TFs and genes

To compare the regulatory networks between humans and mice, we constructed the mouse-specific miRNA–TF–gene regulatory network by integrating all three types of mouse motifs and found that this network consisted of 160 nodes (92 miRNAs, 20 TFs and 48 non-TF genes) and 824 interactions (edges) (Supplementary Material, Fig. S2). Using MCC scores, six hub nodes (six TFs: FOXD3, SP8, PAX9, PBX2, PBX3 and ZEB1) were extracted, and the molecular targets were identified for each hub TF (Supplementary Material, Table S13). Next, we conducted KEGG pathway enrichment analysis of these hub TFs (FOXD3, SP8, PAX9, PBX2, PBX3 and ZEB1) and found that a set of three hub TFs (FOXD3, PAX9 and SP8) was significantly enriched in various pathways that are closely associated with CL-specific biological processes (Benjamini-Hochberg’s FDR < 0.05) (Fig. 3A). For instance, targets of PAX9 were enriched in the RAP1 signaling pathway and adherens junction (FDRs = 2.76 × 10−2 and 7.19 × 10−3, respectively), and targets of FOXD3 were enriched in the adherens junction (FDR = 2.39 × 10−2) (Fig. 3B). All of these pathways have been previously reported to be associated with CL (76).

Identification of potential mouse CL-associated TFs and miRNAs in the Hedgehog signaling pathway

The Hedgehog signaling pathway plays critical roles in the formation of the lip and secondary palate (77). Therefore, we investigated the association of candidate miRNAs, TFs, and non-TF genes with the Hedgehog signaling pathway in mouse CL (Fig. 4C) and found that 15 miRNAs (miR-129, miR-132, miR-144, miR-195, miR-203, miR-214, miR-30e, miR-320, miR-322, miR-370, miR-376b, miR-410, miR-493, miR-539 and miR-9) directly regulated CL-associated genes. For example, in module 1 (Fig. 4D), TF GLI2 regulated miR-376b and FGFR1 (two consensus interaction pairs).

Conserved miRNA–TF–non-TF gene regulatory network in human and mouse CL

We identified a consensus subnetwork through comparison of the miRNA–TF–non-TF target gene regulatory network between humans and mice (Fig. 5). This overlapped network consisted of five miRNAs, eight TFs, and one gene. Among the five miRNAs in the overlapped subnetwork, miR-27b regulated three TFs (PAX3, PAX9 and SATB2) and two non-TF genes (RARA and SUMO1), miR-133b regulated one TF (PAX7) and two non-TF genes (FGFR1 and SUMO1), miR-205 regulated two TFs (PAX9 and SATB2) and one non-TF genes (RARA), miR-376b regulated one TF (GLI2) and one non-TF gene (FGFR1), and miR-376c regulated one TF (GLI2) and one non-TF gene (SATB2) (Fig. 5).

Figure 5

Consensus subnetwork between humans and mice. This subnetwork was extracted from the whole human miRNA–TF–non-TF target gene regulatory network and the mouse miRNA–TF–non-TF target gene regulatory network.

Experimental validation of miRNA–gene regulatory networks for miR-27b, miR-133b and miR-205 in human and mouse lip mesenchymal cells

Our bioinformatic analyses highlighted miR-27b, miR-133b, miR-205, miR-376b and miR-376c as potentially associated with CL through the suppression of CL-associated genes. To test the function of these miRNAs in the regulation of cell proliferation, we conducted cell proliferation assays with a mimic for either miR-27b, miR-133b, miR-205, miR-376b or miR-376c in cultured human embryonic lip mesenchymal (HELM) cells and mouse embryonic lip mesenchymal (MELM) cells. We found that treatment with miR-27b, miR-133b and miR-205 mimic suppressed cell proliferation in HELM cells (Fig. 6A), and overexpression of miR-27b, miR-133b, miR-205 and miR-376b suppressed cell proliferation in MELM cells (Fig. 6B). To assess the anticorrelation between miRNA and its target genes, we conducted quantitative RT-PCR for the predicted target genes of miR-27b, miR-133b or miR-205. We found that treatment with miR-27b mimic downregulated expression of PAX3, PAX9, RARA and SUMO1 in HELM cells (Fig. 7A) and Pax9, Rara and Satb2 in MELM cells (Fig. 7B), whereas miR-133b mimic downregulated FGFR1, PAX7 and SUMO1 in HELM cells (Fig. 7C) and Fgfr1, Pax7 and Sumo1 in MELM cells (Fig. 7D). In addition, miR-205 mimic downregulated PAX9 and RARA in HELM cells (Fig. 7E) and Pax9, Rara and Satb2 in MELM cells (Fig. 7F). Taken together, our results indicate that the miR-27b–PAX9/RARA, miR-133b–PAX7/FGFR1/SUMO1 and miR-205–PAX9/RARA networks are conserved in humans and mice.

Figure 6

Experimental validation for effects on cell proliferation through overexpression of miR-27b-3p/5p, miR-133b, miR-205-3p/5p, miR376c-3p/5p and miR-376b-3p/5p in HELM and MELM cells. (A) Cell proliferation assays using HELM cells treated with the indicated miRNA mimic. (B) Cell proliferation assays using MELM cells from E11.5 maxillary process treated with the indicated miRNA mimic. *P < 0.05, **P < 0.01, ***P < 0.001. Each treatment group was compared to the controls.

Figure 7

Experimental validation for target gene regulation by miR-27b-3p/5p, miR-133b and miR-205-3p/5p in HELM and MELM cells. (A–F) Quantitative RT-PCR for the indicated genes after treatment with miR-27b-3p/5p mimic or control in HELM cells (A) and MELM cells (B), after treatment with miR-133b mimic or control in HELM cells (C) and MELM cells (D), and after treatment with miR-205-3p/5p mimic or control in HELM cells (E) and MELM cells (F). **Adjusted P < 0.01, *** adjusted P < 0.001 versus control (n = 6).

Discussion

In this study, we constructed miRNA–TF–non-TF target gene networks using potential CL-associated TFs, miRNAs and non-TF target genes conserved in humans and mice and found that five miRNAs (miR-27b, miR-133b, miR-205, miR-376b and miR-376c) and their consensus regulatory loops may play essential roles in CL pathogenesis through bioinformatic analyses. In addition, we further evaluated the functional relevance of five miRNAs in human and mouse palatal mesenchymal cells. Among them, we found that three miRNAs (miR-27b, miR-133b and miR-205) and their networks were involved in cell proliferation through downregulation of target genes. Our bioinformatic analyses and experimental validation suggest that miRNA–TF–non-TF gene networks can explain one of the pathological mechanisms of CL.

Our previous study suggested that hsa-miR-27b might activate the Wnt signaling pathway, leading to human CL (58) and human and mouse cleft palate (78). In addition, miR-27b expression is upregulated at E11.5, compared to E10.0, in the mouse nasal process and maxillary process during embryogenesis (36). In this study, we found that miR-27b suppressed cell proliferation by suppressing PAX9 and RARA expression, which are associated with the Wnt pathway (79). Therefore, miR-27b-regulated Wnt signaling may play a role in CL.

miR-133b plays a vital role during craniofacial development. Expression of miR-133b in zebrafish results in clefts of the ethmoid plate (80), whereas its expression in the maxillary process in mice at E12.5 is weak (36,80). Overexpression of miR-133b suppresses cell proliferation through downregulation of the expression of the FGFR1, GCH1, PAX7, SMC2 and SUMO1 genes in human palatal mesenchymal cells (81). In this study, we found that miR-133b suppresses cell proliferation in lip mesenchymal cells through downregulation of FGFR1, PAX7 and SUMO1. These results indicate that there is high similarity in target genes between the lip and palate. FGFR1, a receptor tyrosine kinase whose ligands are specific members of the fibroblast growth factor family, regulates cell proliferation, cell migration and cell differentiation (82). In mice, the CNC-specific deletion of Fgfr1 (Wnt1-Cre;Fgfr1 cKO mice) results in cleft palate (83). In humans, SNPs or mutations in the FGFR1 gene are associated with non-syndromic CL/P (84,85). PAX7, a TF containing a paired box domain, an octapeptide, and a paired-type homeodomain, plays roles in craniofacial development (86), and PAX7 nucleotide variants are associated with non-syndromic CL/P (87). SUMO1, an ubiquitin-like protein, plays roles in various cellular processes, such as nuclear transport, transcriptional regulation, and apoptosis, through sumoylation of many genes. For example, CL/P-associated genes TBX22, MSX1, SATB2, P63, PAX9, TRPS1 and EYA1 are essential for sumoylation leading to their functional activation (88). There are SNPs in SUMO1 associated with non-syndromic CL/P in humans (89–92). Genetically engineered mice generated using a gene-trap strategy in RRQ016 ES cells exhibit a mutation in Sumo1 and cleft palate (93). The detailed analysis in these mice revealed that there was unexpected chromatin rearrangement in intron 1 of Sumo1, as well as a 3-kbp duplication of the preintegration site of β-galactosidase. Interestingly, despite of the rearrangement in Sumo1, expression of the SUMO1 protein was not affected. This suggests that cleft palate in the RRQ016-derived mouse line is caused by mutations in other genes or in additional genes associated with CL and palate (94). This is further supported by recent findings that Sumo1 null mice (95), and mice generated through a gene-trap strategy in XA024 ES cells (94), exhibit no developmental defects and no alternation in sumoylation with complete depletion of SUMO1. In agreement with these findings, recent studies showed that sumoylation by SUMO1 is compensated by homologs SUMO2 and SUMO3 (94,95).

miR-205 expression is upregulated at E11.5, compared to E10.0, in the mouse nasal process and maxillary process during development (36) and is associated with the TGF-β signaling pathway (96). In humans, expression of hsa-miR-205-5p is higher in individuals with non-syndromic CL/P compared to healthy people (33). Since both miR-205 and miR-27b can suppress cell proliferation through downregulation of PAX9 and RARA in human and mouse lip mesenchymal cells, the suppression of these genes through miR-205 and miR-27b might be associated with the pathogenesis of CL.

Although our approach is useful to identify causative molecular mechanisms, it presents several limitations. On one hand, we have data limitations. We used the CL-associated gene set to build the networks, but this set included highly heterogeneous data. On the other hand, there is a statistical power limitation. Although we used stringent criteria to develop the miRNA–TF-non-TF target gene networks, some of the interactions might include false positive or negative results. However, we systematically examined the miRNA–TF–non-TF target gene regulation in CL using all the available data we could collect. Therefore, we believe this work is applicable, as a general procedure, for constructing a conserved regulatory network for other diseases.

Materials and Methods

Collection of CL-associated genes and miRNAs

The list of CL-associated genes was identified in a systematic review and Mouse Genome Informatics database search of genes associated with CL in mice (59) and humans (58,60). A list of 173 human genes (31 TFs and 142 non-TF genes) and another 131 mouse genes (37 TFs and 94 non-TF genes) was obtained after manual curation. Through an extensive literature search, 46 human miRNAs (33,97–103) and 30 mouse miRNAs (34,67,101,104–107) were collected using the following conditions: (1) miRNAs that target CL-associated genes [obtained from genome-wide association studies], (2) those that are differentially expressed between CL and control samples and (3) those that are differentially expressed at the different stages of lip development.

Bioinformatic analyses

PITA (version 6), miRanda (Release August 2010), miRTarBase (version 7) and TargetScan (Release 7.1) (108–111) were used for the prediction of target genes of the miRNAs. The miRNA–target pairs identified in at least two databases mentioned above were retained for further analysis. The UCSC Table Browser (112) was used to extract the promoter sequences (−1000 to +200 bp of the transcription start site) of human and mouse GENCODE protein-coding genes and precursor miRNAs (wgEncodeGencodeBasic V27 for human; wgEncodeGencodeBasic VM16 for mouse), as previously described (56). TRANSFAC® Professional version (Release 2016.1) (113) was used for the identification of TF motifs, as previously described (78). Significant hits with P-value <1 × 10−5 that occurred within the promoter of GENCODE genes were retained.

Identification of motifs (FFLs and FBLs)

Our in-house MotifPredictor, available at https://www.uth.edu/bioinfo/software.htm and https://github.com/emanlee/MotifPredictor, was used to extract three types of motifs (Fig. 1C). Briefly, in the type I motif, a TF regulates a miRNA and a non-TF gene at the transcriptional level, whereas the miRNA regulates the TF and the non-TF gene at the post-transcriptional level. In the type II motif, a TF regulates a non-TF gene at the transcriptional level, whereas a miRNA regulates the non-TF and the TF gene at the post-transcriptional level. In the type III motif, both a TF and a miRNA regulate their consensus target non-TF gene, whereas the miRNA and TF mutually regulate each other (Fig. 1C). A regulatory network was built by integrating all these motifs. Network visualization and analysis were conducted by using Cytoscape software version 3.6 (http://cytoscape.org/) (114). After constructing the network, centrality measure and MCC were applied through the cytoHubba plugin (version 0.1) implemented in Cytoscape (version 3.6) (69). MCL (75) was applied for extracting modules from the miRNA–TF–target gene regulatory networks. To identify pathways enriched in the target genes and TFs, the WebGestalt webtool (70) was used for gene set enrichment analysis with KEGG pathway annotations; the DAVID online tool was used for GO term annotations (63).

Cell proliferation assay

HELM cells, a cell line originated from human lip cancer tissue, were purchased from JCRB Cell Bank (#JCRB9103KD) and maintained in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% fetal bovine serum (FBS), L-glutamine and penicillin/streptomycin. HELM cells were used by passage 5–10. Primary MELM cells were isolated from the maxillary process, a developing lip region, of E11.5 C57BL6/J mice, as previously described (59). MELM cells were maintained in DMEM supplemented with 10% FBS, penicillin/streptomycin, L-glutamine, β-mercaptoethanol and non-essential amino acids. MELM cells were used by passage 2. Cells were plated at a density of 10 000 (HELM cells) or 5000 (MELM cells) per well on 96-well plates and transfected with a mimic for the (negative control, miR-27b-3p/5p, miR-133b, miR-205-3p/5p, miR-376b-3p/5p and miR-376c-3p/5p (mirVana miRNA mimic, Thermo Fisher Scientific) using the lipofectamine RNAiMAX transfection reagent (Thermo Fisher Scientific). Cell proliferation assays were conducted using the Cell Counting Kit 8 (Dojindo Molecular Technology) (n = 6 per group), as previously described (59,60).

Quantitative RT-PCR

Total RNAs were isolated from cultured HELM and MELM cells, which were transfected with a mimic for the target miRNAs for 24 h (n = 6 per group), using the QIAshredder and RNeasy mini extraction kit (QIAGEN), as previously described (59,60). The PCR primers used are listed in Supplementary Material, Tables S14  and  S15.

Statistical analysis

The statistical significance for multiple groups of two was evaluated using a multiple t-test with Bonferroni adjustment. Multiple comparisons in cell proliferation assays were evaluated using a two-way analysis of variance with Dunnett’s test. A P-value <0.05 was considered statistically significant. The FDR was calculated to adjust the significance of the multiple tests. For all groups, data are represented as mean ± standard deviation.

Funding

National Institute of Dental and Craniofacial Research, National Institutes of Health (R03DE028103, R03DE027393, and R01LM012806 to ZZ; and R01DE026767, R01DE029818, R03DE026509, and R03DE028340 to JI).

Authors’ Contributions

ZZ and JI designed the project. AL conducted the bioinformatic analyses. HY and SSR performed the experiments. AS collected genes related to CL through a database search. HY, AL, AS, ZZ and JI wrote the manuscript. All authors read and approved the final manuscript.

Conflict of Interest statement. The authors have declared that no competing interests exist.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.