Integrated assessment of differentially expressed plasma microRNAs in subtypes of nonsyndromic orofacial clefts.

BACKGROUND: Orofacial clefts include cleft lip only (CLO), cleft palate only (CPO), and cleft lip with palate (CLP). Previously, we reported the expression profile of plasma microRNAs in CLO, CPO, and CLP, respectively. However, the interaction of each subtype remains poorly investigated.
METHODS: In this study, we integrated the expression profiles of plasma miRNAs in these 3 subtypes, and assessed the distinct and overlapping dysregulated miRNAs using Venn diagrams. Their respective target genes reported in the literature were further analyzed using pathway analysis. RESULTS AND
CONCLUSION: The results showed that distinct or overlapping signaling pathways were involved in CLO, CPO, and CLP. The common key gene targets reflected functional relationships to the Wnt, Notch, TGF-beta, and Hedgehog signaling pathways. Further studies should examine the mechanism of the potential target genes, which may provide new avenues for future clinical prevention and therapy.

Introduction

As common, complex birth defects with an overall birth prevalence of 1:800 live births worldwide,[ orofacial clefts (OFCs) are categorized as cleft lip only (CLO) as usual, cleft palate only (CPO), and cleft lip with palate (CLP). The overall prevalence of OFCs in China was 1.4 per 1000 live births.[ Around 70% of the whole of clefts were categorized as nonsyndromic, with no identifiable structural imperfection except for the clefts.[ Children born with OFCs may encounter many difficulties in adapting to future social life, including feeding, language communication, dental development, hearing, and so on. Those influenced families will face tremendous financial burden and psychological distress from the community. The genetic causes of nonsyndromic clefts have been clarified gradually with the development of various ways including linkage, association studies, GWAS, exome sequencing, copy number variation, and whole-genome sequencing.[ In addition, epigenetic data are able to assist in recognizing risk elements for OFCs or causal effects of clefts by disclosing epigenetic mechanisms, or by catching information about causal genetic or environmental points.[

As epigenetic factors, plasma miRNAs have emerged as attractive diagnostic and prognostic biomarkers for cancer, diabetes, allergic asthma, and oral cancers.[ The identification of the plasma miRNA and potential targets’ gene ontology and KEGG pathway enrichment is important to understand the regulatory mechanisms of CLO,[ CPO, and CLP.[ However, the interaction of each subtype remains poorly investigated. Recently, mounting evidence suggests there may be unique underlying pathophysiology and/or genetic modifiers influencing expression of these 3 subtypes of OFCs including CLO, CPO, and CLP.[ Therefore, we explored whether children with different cleft subtypes showed distinct plasma miRNAs’ profiles. In this study, we analyzed the expression profiles of plasma microRNAs in nonsyndromic CLO, CPO, and CLP as a whole using our previously reported human miRNA microarray data.[ Based on KEGG pathway analysis, we found that there were distinct biological processes involving in the physiology of CLO, CPO, and CLP.

Materials and methods

Sample collection

Institutional Review Board of Nanjing Medical University (2014-10-16) has approved this research. Nonsyndromic cleft lip (CLO), cleft palate (CPO), and cleft lip with palate (CLP) patients’ plasma samples were collected using methods as we previously described.[ EDTAK2 tubes (regular type) were utilized to collect the peripheral blood. The supernatant plasma was centrifuged at 12,000 g for 10 minutes. The supernatant was aliquoted and stored at −80°C.

Human miRNA microarray analysis

Twelve plasma samples were divided into 4 groups including CLO (n = 3), CPO (n = 3), CLP (n = 3), and control (n = 3). These 4 groups were analyzed using human miRNA microarray chips (8×60K) v21.0 as we previously described.[ miRNAs that show expression changes above twofold were selected. We draw Venn diagrams on this website (http://bioinformatics.psb.ugent.be/webtools/Venn/).

KEGG pathway enrichment analyses

We use miRTarbase and Tarbase databases through the ShanghaiBio Corporation (SBC) online (http://sas.ebioservice.com) to predict the target genes of these differentially expressed plasma miRNAs. KEGG pathway enrichment analyses were performed as we previously described.[ Fisher's exact test was utilized to conduct the calculation of the enrichment P-values of the pathway enrichment analyses, which was rectified applying enrichment q-values (the false discovery rate) that were computed utilizing the method of John Storey.

Results

miRNA microarray and bioinformatics analysis

The 4 groups CLO, CPO, CLP, and control of plasma miRNA microarray analysis were displayed using hierarchical clustering (Fig. 1A). Scatter plot visualization showed miRNAs expression variation between the 2 groups of control and CLO, control and CPO, control and CLP, CLO, and CPO, CLO, and CLP, CPO and CLP (Fig. 1B). Red dots represent a twofold or greater upregulation between the 2 groups, green dots refer to a twofold or greater downregulation between the 2 groups (Fig. 1B). Upregulated or downregulated miRNAs’ numbers were shown in Figure 1B. Correlation analysis of these expressed miRNAs in the control, CLO, CPO, and CLP plasma samples revealed that CLO was closely related to CLP (Fig. 1C). The correlation coefficient between CLO and CLP was 0.9749, meaning that CLO was highly correlated with CLP.

Figure 1

Plasma miRNA microarray expression data among the patients with nonsyndromic CLO, CPO, CLP, and healthy children. (A) Hierarchical clustering reveals the miRNA expression profile. (B) Scatter plot represents the differentially expressed miRNAs between indicated groups. The miRNAs numbers are listed above the data dot. (C) Correlation analysis shows the relationship between indicated groups. (D) Venn diagrams show the number of distinct and overlapping upregulated or downregulated miRNAs among CLO, CPO, and CLP. CLO = cleft lip only, CLP = cleft lip with palate, CPO = cleft palate only, miRNAs = microRNAs.

We selected differentially expressed plasma miRNAs that showed expression changes above twofold (compared to the control group) in the CLO, CPO, and CLP groups. Overlapping data sets were performed using Venn diagrams (Fig. 1D). Totally, 57 miRNAs were co-upregulated among the CLO, CPO, and CLP plasma samples (Fig. 1D). Table 1 shows all those co-upregulated miRNAs (fold change≥2). Among these CLO, CPO, and CLP plasma samples, 45 miRNAs were co-downregulated (Fig. 1D). Table 2 shows all these co-downregulated miRNAs (fold change≥2).

Table 1

List of 57 miRNAs that were co-overexpressed in CLO, CPO, and CLP plasma samples.

Table 2

List of 45 miRNAs that were co-downregulated in CLO, CPO, and CLP plasma samples.

Previously, we have demonstrated the expression changes of several miRNAs in larger samples of CLO, CPO, and CLP using BulgeLoop qRT-PCR, and confirmed that the miRNA microarray data was confidential.[ Therefore, we performed bioinformatics analysis as follows using these differentially expressed miRNAs in respective groups.

Functional analysis of potential target genes of co-expressed dysregulated miRNAs

We want to examine whether there are overlapping biological processes involved among CLO, CPO, and CLP groups. Therefore, we performed KEGG pathway analysis. Around 57 co-upregulated and 45 co-downregulated miRNAs among CLO, CPO, and CLP groups could produce 13,846 nonduplicated target genes reported in the literature using the miRTarbase and Tarbase (Supplemental file 1). The top 30 enriched pathways for these 13,846 nonduplicated target genes are listed in Figure 2. Overall, these 13,846 target genes were involved in Wnt, p53, TGF-beta, and Hedgehog signaling pathways (Fig. 3).

Figure 2

KEGG pathway analyses show the associated function of the target genes of the miRNAs dysregulated among CLO, CPO and CLP. Top 30 enriched pathways for the 13846 target genes of the 57 upregulated miRNAs and 45 downregulated miRNAs among the CLO, CPO and CLP. Fisher's exact test is utilized to calculate the enrichment P-value. Based on the first letter of the path name, the terms/paths on the vertical axis are sorted in descending order. Enrichment factors are plotted on the horizontal axis, that is, the number of dysregulated genes/total number of dysregulated genes in the pathway/(number of genes in the pathway in the database/total number of genes in the database). The first 30 enrichment approaches are based on enrichment factors. Pathway in the number of genes ≥ 4, P < .05 is the selection criteria. P-values are expressed in different colors from green to red. The different sizes of the round shapes represent the gene count number in a pathway. CLO = cleft lip only, CLP = cleft lip with palate, CPO = cleft palate only, KEGG = Kyoto Encyclopedia of Genes and Genomes, miRNAs = microRNAs.

Figure 3

KEGG pathway analyses show the associated function of the target genes of the miRNAs dysregulated in CLO. The top 30 enriched pathways for the 14114 target genes of the 94 upregulated miRNAs and 22 downregulated miRNAs in CLO. The enrichment P-values were calculated using Fisher's exact test. The term/pathway on the vertical axis was drawn according to the first letter of the pathway name in descending order. The horizontal axis represents the enrichment factor, that is, (the number of dysregulated genes in a pathway/the total number of dysregulated genes)/(the number of genes in a pathway in the database/the total number of genes in the database). The top 30 enriched pathways were selected according to the enrichment factor value. The selection standards were the number of genes in a pathway ≥4 and P < .05. The different colours from green to red represent the P-value. The different sizes of the round shapes represent the gene count number in a pathway. CLO = cleft lip only, miRNAs = microRNAs, KEGG = Kyoto Encyclopedia of Genes and Genomes.

Function analysis of the potential target genes of the miRNAs dysregulated in CLO, CPO, or CLP, respectively

To explore whether distinct biological processes regulated CLO, CPO, and CLP, we performed KEGG pathway analysis of the forecasted targets of the distinctively expressed miRNAs in CLO, including 94 upregulated miRNAs and 22 downregulated miRNAs yielding 14,114 nonduplicated potential target genes (Supplemental file 2). The top 30 enriched KEGG pathways for the 14,114 potential target genes of the dysregulated miRNAs in CLO showed that those genes were mainly involved in Wnt, Notch, TGF-beta, and Hedgehog signaling pathways (Fig. 4). Then, we performed KEGG analysis of the forecasted targets of the distinctively expressed miRNAs in CPO, including 10 upregulated miRNAs and 3 downregulated miRNAs producing 9330 nonduplicated possible target genes (Supplemental file 3). The top 30 enriched KEGG pathways are listed in Figure 5. The analysis revealed that the potential target genes of the dysregulated miRNAs in CPO were associated with the Wnt, p53, Notch, and Hedgehog signaling pathways (Fig. 5). Additionally, we performed KEGG analysis of the forecasted targets of the distinctively expressed miRNAs in CLP, including 23 upregulated miRNAs and 6 downregulated miRNAs creating 11,232 nonduplicated possible target genes (Supplemental file 4). The top 30 enriched KEGG pathways revealed that the potential target genes of the dysregulated miRNAs in CLP were associated with the Wnt, TGF-beta, Notch, and Hedgehog signaling pathways (Fig. 2).

Figure 4

KEGG pathway analyses show the associated function of the target genes of the miRNAs dysregulated in CPO. Top 30 enriched pathways for the 9330 target genes of the 10 upregulated miRNAs and 3 downregulated miRNAs in CPO. Fisher's exact test is utilized to calculate the enrichment P-value. Based on the first letter of the path name, the terms/paths on the vertical axis are sorted in descending order. Enrichment factors are plotted on the horizontal axis, that is, the number of dysregulated genes/total number of dysregulated genes in the pathway/(number of genes in the pathway in the database/total number of genes in the database). The first 30 enrichment approaches are based on enrichment factors. Pathway in the number of genes ≥ 4, P < .05 is the selection criteria. P-values are expressed in different colors from green to red. The gene counts in the pathway are displayed by circles of various sizes. CPO = cleft palate only, KEGG = Kyoto Encyclopedia of Genes and Genomes, miRNAs = microRNAs.

Figure 5

KEGG pathway analyses show the associated function of the target genes of the miRNAs dysregulated in CLP. Top 30 enriched pathways for the 11232 target genes of the 23 upregulated miRNAs and 6 downregulated miRNAs in CLP. Fisher's exact test is utilized to calculate the enrichment P-value. Based on the first letter of the path name, the terms/paths on the vertical axis are sorted in descending order. Enrichment factors are plotted on the horizontal axis, that is, the number of dysregulated genes/total number of dysregulated genes in the pathway/(number of genes in the pathway in the database/total number of genes in the database). The first 30 enrichment approaches are based on enrichment factors. Pathway in the number of genes ≥ 4, P < .05 is the selection criteria. P-values are expressed in different colors from green to red. The gene counts in the pathway are displayed by circles of various sizes. CLP = cleft lip with palate, KEGG = Kyoto Encyclopedia of Genes and Genomes, miRNAs = microRNAs.

To clarify their relationship among CLO, CPO, and CLP, we illustrated the overlapping top 30 KEGG pathways among the 3 data sets using Venn diagrams (Fig. 6). As shown in Figure 6, 12 common KEGG pathways were enriched in the 3 CLO, CPO, and CLP groups, including Wnt, Notch, and Hedgehog signaling pathway (Fig. 6). Ten specific KEGG pathways were involved in CLO, including Valine, leucine, and adipocytokine signaling pathway (Fig. 6). Ten specific KEGG pathways were involved in CPO, including basal transcription factors, Lysine biosynthesis, and neurotrophin signaling pathway (Fig. 6). Whereas 10 specific KEGG pathways were involved in CLP, including regulation of autophagy, insulin signaling pathway, and mTOR signaling pathway (Fig. 6).

Figure 6

Venn diagrams show the number of distinct and overlapping top 30 KEGG pathways among CLO, CPO, and CLP. The detailed distinct or overlapping KEGG pathways’ information is listed at the right. CLO = cleft lip only, CLP = cleft lip with palate, CPO = cleft palate only, KEGG = Kyoto Encyclopedia of Genes and Genomes.

Discussion

To our knowledge, plasma miRNAs are highly stable.[ Exploring the functions of plasma miRNAs in CLO, CPO, and CLP will help us to understand the physiology and development of these 3 diseases. In this study, we used microarray profiling to deeply assess the interactions and distinctions of plasma miRNAs among CLO, CPO, and CLP patients. A total of 207 miRNAs were upregulated and 98 miRNAs were downregulated in the CLO plasma samples compared to the control. 77 miRNAs were upregulated and 55 miRNAs were downregulated in the CPO plasma samples compared to the control. Around 138 miRNAs were upregulated and 83 miRNAs were downregulated in the CLP plasma samples compared to the control. Further elucidating these plasma miRNAs’ roles in orofacial clefts would provide clues for preventing OFCs.

A lot of evidence has proved that the etiological mechanisms of CLP and CPO are distinct.[ Some studies reveal that nonsyndromic orofacial clefts have identical configurations.[ Recently, CPO was reported to have unrelated genetic factors and should be examined independently.[ In our study, we found that CLO was highly correlated with CLP based on correlation analysis.

Venn diagrams were new method to visualize the complex genetic set relations.[ Here, we found that 57 miRNAs were co-upregulated in CLO, CPO, and CLP plasma samples, whereas 45 miRNAs were co-downregulated in CLO, CPO, and CLP plasma samples.

One limitation of our study was that we performed 4 miRNA microarray chips using 12 samples. Therefore, no P-values were got properly. Thus, these results in Tables 1 and 2 are reported in fold change values instead of P-values. The other limitation was that we used small samples (N = 3 for each group) for microarray chip assay. Whereas 6 selected plasma miRNAs including miR-16-2-3p, miR-215-5p, miR-365a-3p, miR-574-5p, miR-584-5p, and miR-877-5p were validated using qRT-PCR between another 13 CLO and 11 healthy children as we reported previously.[ We also checked another 6 plasma miRNAs including miR-340-5p, miR-877-5p, miR-3648, miR-1260a, miR-494-3p, and miR-1304-3p among 16 CPO, 33 CLP, and 8 healthy children using qRT-PCR.[ The results revealed that qRT-PCR verification data were consistent with the miRNAs microarray data.[ Thus, our miRNAs microarray data are reliable and we could conduct bioinformatics analysis for key genes and pathways participated in CLO, CPO, and/or CLP pathogenesis based on the microarray data.

One of the crucial pathological reasons of orofacial clefts is unusual craniofacial development. MiRNAs played important roles in craniofacial development.[ In addition, SNPs in the miRNA-binding sites or miRNA processing genes were announced to be linked to nonsyndromic orofacial clefts susceptibility.[ In mammals, it is possible for miRNAs to control 30% of the protein-coding genes by means of posttranscriptional silencing; therefore, the dysregulation of miRNAs in CLO, and it is possible for CPO or CLP patients have a far-reaching impact on a wide range of roles in biology. Based on our miRNA microarray data, we found the intersection of differentially expressed miRNAs in CLO, CPO, and CLP. Our results demonstrate the existence of overlap and distinct signaling pathways for which type of OFCs develops and suggest plausible elements responsible for phenotypic heterogeneity, revealing the complex genetic architecture of OFCs. Additionally, it will set up the foundation upon which would be possible for us to discover how the prior diagnostic panels are designed.

Wnt signaling pathway including wnt9a and irf6 has been reported to play key roles in craniofacial morphogenesis.[ IRF6 and Notch ligand Jagged2 signaling were crucial for restricting palatal adhesion and fusion competence.[ Craniofacial abnormalities may be caused by changes in transforming growth factor-beta (TGFβ) signaling.[ TGFβ1 and TGFβ2 are believed to be primarily engaged in controlling the proliferation of the palatal mesenchyme whereas TGFβ3 has been suggested to have a significant impact on controlling the fate of midline epithelial cells during palatal fusion.[ Perturbed hedgehog (HH) signaling particularly mutation in patched1 plays a major role in craniofacial development.[ Based on the KEGG pathway analysis, the predicted target genes of the upregulated and downregulated miRNAs in CLO, CPO, and CLP mainly participated in Wnt, Notch, TGF-beta, and Hedgehog signaling pathways, that is consistent with existing literatures.

Epigenetics is explained as genetic alterations in phenotype or gene expression that result from processes other than alterations in the basic DNA sequence. These processes contain, but are not restricted in, DNA methylation, microRNA impacts, and histone adjustments that change chromatin conformation. Recent studies suggest that different methylation profiles was found in distinct orofacial cleft (OFC) subtypes, representing a prospects with potential in discovering the possible impacts of epigenetic modifications on the aetiology of OFCs and/or as clinically functional biomarkers of OFC subtypes.[

Our study found that distinct plasma miRNAs were differentially expressed in different orofacial cleft (OFC) subtypes. Further analyzing the predictive target genes of nonoverlapping miRNAs in CLO, CPO, and CLP using KEGG pathway analysis, we found that the top 30 enriched enriched pathways were uniquely associated with CLO including nonhomologous end-joining, protein processing in endoplasmic reticulum, glycosphingolipid biosynthesis, valine, leucine and isoleucine biosynthesis, RNA polymerase, endocytosis, sulfur relay system, adipocytokine signaling pathway, protein export, and fatty acid biosynthesis. In contrast, 10 pathways of the top 30 enriched pathways were only linked with CPO, including colorectal cancer, vasopressin-regulated water reabsorption, porphyrin and chlorophyll metabolism, melanogenesis, arrhythmogenic right ventricular cardiomyopathy (ARVC), prostate cancer, renal cell carcinoma, basal transcription factors, lysine biosynthesis, and neurotrophin signaling pathway. Moreover, 10 pathways of the top 30 enriched pathways were only linked with CLP, including regulation of autophagy, caffeine metabolism, insulin signaling pathway, pancreatic cancer, bladder cancer, mTOR signaling pathway, lysine degradation, cyanoamino acid metabolism, hypertrophic cardiomyopathy (HCM), and dilated cardiomyopathy. Further elucidating the complex genetic architecture of OFCs may push us following the fundamental humanitarian aim of stopping orofacial clefts.

Taken together, we provide an integration analysis of these differentially expressed plasma miRNAs in subtypes of nonsyndromic orofacial clefts. Distinct and overlapping signaling pathways were revealed to participate in the pathogenesis of CLO, CPO, and CLP. In the future, urgent need should be taken to examine the mechanism of potential target genes. This may provide new targets for future clinical prevention and treatment.

Author contributions

JZ and WS planned the experiment and conducted the data analyses. NW and JY performed the sample preparation and bioinformatics analyses. NW, JY, and TH performed the statistical analyses. NW and WS wrote and edited the manuscript.

Conceptualization: Weimin Shen.

Data curation: Tao Han, Weimin Shen.

Formal analysis: Tao Han.

Funding acquisition: Jijun Zou.

Investigation: Ni Wu.

Methodology: Ni Wu.

Project administration: Jun Yan.

Software: Jun Yan.

Supervision: Jijun Zou.

Validation: Jijun Zou.

Writing – original draft: Ni Wu.

Writing – review & editing: Weimin Shen.