Gene expression profiling analysis contributes to understanding the association between non-syndromic cleft lip and palate, and cancer.

The present study aimed to investigate the molecular mechanisms underlying non‑syndromic cleft lip, with or without cleft palate (NSCL/P), and the association between this disease and cancer. The GSE42589 data set was downloaded from the Gene Expression Omnibus database, and contained seven dental pulp stem cell samples from children with NSCL/P in the exfoliation period, and six controls. Differentially expressed genes (DEGs) were screened using the RankProd method, and their potential functions were revealed by pathway enrichment analysis and construction of a pathway interaction network. Subsequently, cancer genes were obtained from six cancer databases, and the cancer‑associated protein‑protein interaction network for the DEGs was visualized using Cytoscape. In total, 452 upregulated and 1,288 downregulated DEGs were screened. The upregulated DEGs were significantly enriched in the arachidonic acid metabolism pathway, including PTGDS, CYP4F2 and PLA2G16; and transforming growth factor (TGF)‑β signaling pathway, including SMAD3 and TGFB2. The downregulated DEGs were distinctly involved in the pathways of DNA replication, including MCM2 and POLA1; cell cycle, including CDK1 and STAG1; and viral carcinogenesis, including PIK3CA and HIST1H2BF. Furthermore, the pathways of cell cycle and viral carcinogenesis, with higher degrees of interaction were found to interact with other pathways, including DNA replication, transcriptional misregulation in cancer, and the TGF‑β signaling pathway. Additionally, TP53, CDK1, SMAD3, PIK3R1 and CASP3, with higher degrees, interacted with the cancer genes. In conclusion, the DEGs for NSCL/P were implicated predominantly in the TGF‑β signaling pathway, the cell cycle and in viral carcinogenesis. The TP53, CDK1, SMAD3, PIK3R1 and CASP3 genes were found to be associated, not only with NSCL/P, but also with cancer. These results may contribute to a better understanding of the molecular mechanisms of NSCL/P.

Introduction

Non-syndromic cleft lip, with or without cleft palate (NSCL/P) is one of the most common types of congenital defect and affects 3.4-22.9/10,000 individuals worldwide (1). The interaction between environmental and genetic factors during embryonic development has been identified as the determinant pathogeny of NSCL/P (2).

In previous years, common alleles affecting the susceptibility to this complex disease have been identified using genome-wide association studies (3,4). However, each variant has a low incidence in NSCL/P, which introduces difficulty in determining the expected heritability for the disease (5). There is sufficient evidence that variants in interferon regulatory factor 6 (IRF6) have a substantial impact on the occurrence of NSCL/P (6). For example, a single nucleotide polymorphism (rs642961; G>A) located within an enhancer ~10 kb upstream of the IRF6 transcription initiation site is significantly over-transmitted in NSCL/P, which can disrupt the binding site of transcription factor AP-2α (7), which is a mutation in the autosomal dominant NSCL/P. In addition, mutations of MAFB, ABCA4 (8), VAX1 (9), FGFR2 (10) and SUMO1 (11), as well as the perturbation of the methionine and folate pathways (12), and haplotypes in the Wnt and fibroblast growth factor signaling pathway (13) have all been confirmed to increase the risk of NSCL/P.

In addition, anomalies in cell migration, proliferation, transdifferentiation and apoptosis are considered to be closely associated to the occurrence NSCL/P (14,15). These events are commonly known to be correlated with cancer. Studies have shown that alterations in certain genes, including WNT (16), MSX1 (17), BMP (18) and BCL3 (19), which are considered to be implicated in carcinogenesis, are also involved in NSCL/P (20–23). In 2013, Kobayashi et al (24) showed that, in NSCL/P dental pulp stem cells, BRCA1 and RAD51, targeted by the E2F1 transcription factor, were dysregulated in the developing embryonic orofacial primordial, and are central to lip and palate morphogenesis. In addition, cellular defences against DNA damage may be involved in determining the susceptibility to NSCL/P, which suggests an etiological overlap between this malformation and cancer (24). However, this previous study predominantly investigated differentially expressed genes (DEGs) associated with DNA double-strand break repair and cell cycle control in NSCL/P group samples, which is less convincing for the hypothesis of an etiological overlap between NSCL/P and cancer.

In the present study, the microarray data deposited by Kobayashi et al were downloaded to further reveal the interplay between the DEGs in NSCL/P samples and cancer genes, and to identify the precise nosogenesis of NSCL/P. Subsequently, pathway enrichment analysis and pathway interaction analysis of the DEGs were performed, and a cancer-associated protein-protein interaction (PPI) network for the DEGs was constructed. The results of these investigations may assist in elucidating the etiology of NSCL/P, and provide more information on the correlation between the mechanisms of NSCL/P and cancer.

Materials and methods

Affymetrix microarray data

The gene expression profile data of GSE42589 (24) were obtained from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/), based on the GPL6244 [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array platform (Affymetrix, Santa Clara, CA, USA; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6244). In total, 13 dental pulp stem cell samples were available for further analysis, including seven dental pulp stem cell samples collected from children with NSCL/P in the exfoliation period, and six controls obtained from healthy children in the exfoliation period. This study was approved by the Biosciences Institute Research Ethics Committee (Protocol 037/2005) at the University of São Paulo (São Paulo, Brazil), and all the patients or legal guardians signed informed-consent documents (22). All samples were cultured in Dulbecco's modified Eagle's medium/Ham's F-12 medium (Life Technologies; Thermo Fisher Scientific, Inc., Waltham, MA, USA) supplemented with 15% fetal bovine serum (HyClone, Logan, UT, USA), 1% non-essential amino acids solution (Life Technologies; Thermo Fisher Scientific, Inc.) and 1% penicillin-streptomycin solution (Life Technologies; Thermo Fisher Scientific, Inc.), in a humidified incubator at 37°C and 5% CO2.

CEL files and the probe annotation files were downloaded, and the gene expression data of all samples were normalized using the Robust Multi-array Average (25) algorithm of the Bioconductor Affy package in R (http://www.bioconductor.org/packages/release/bioc/html/affy.html) (26).

Screening of DEGs

The RankProd method (27) in the Bioconductor package was used to identify genes, which were significantly differentially expressed in the NSCL/P dental pulp stem cells. The raw P-value was adjusted into the false discovery rate (FDR) using the Benjamin and Hochberg method (28), and only the genes within the cut-off criteria of |log2 fold change| >1 and FDR<0.05 were selected as DEGs.

Pathway enrichment analysis

To identify the significant metabolic pathways for the DEGs, the screened DEGs were submitted to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/kegg1.html) for pathway enrichment analysis (29). An FDR <0.1 was used as the cut-off criterion.

Pathway interaction analysis

Pathway interactions were analyzed based on the association between the DEGs in the pathways, which was determined from the protein-protein interaction (PPI) network for the DEGs, obtained from the human protein reference database (http://www.hprd.org/) (30). The pathway interaction network was visualized using Cytoscape (http://cytoscape.org/) (31).

Construction of the cancer-associated PPI network

Cancer genes were obtained from the a database of Functional Census of Human Cancer Genes (http://210.46.85.180:8080/fcensus/) (32), which provides multiple dimension information for cancer genes, including cancer type, cancer gene type, mutation type and mutation frequency, calculated from high-throughput mutational screens of cancer genomes. The R package was used to obtain the interactions between the DEGs and cancer genes. The five DEGs exhibiting the highest degree were selected to construct the PPI network, and the network was visualized using Cytoscape.

Results

Identification of DEGs

Based on the cut-off criteria used for determination of the DEGs, a total of 1,740 DEGs were identified in the NSCL/P samples, including 452 upregulated DEGs and 1,288 downregulated DEGs.

Pathway enrichment analysis of the upregulated and down- regulated DEGs

The upregulated DEGs were significantly enriched in three pathways: Seven DEGs, including TGFB2, TGFB3 and VCAM1, were enriched in the hsa05144 malaria pathway (FDR=6.27E-02); seven DEGs, including PTGDS, PTGIS, CYP4F2, PTGES and PLA2G16 were enriched in the hsa00590 arachidonic acid metabolism pathway (FDR=9.75E-02); and eight DEGs, including ID2, ID4, SMAD3 and TGFB2, were involved in the hsa04350 TGF-β signaling pathway (FDR=9.75E-02; Table I).

Table I

Pathway enrichment analysis for the upregulated and downregulated differentially expressed genes.

Category	Term	Description	Count	FDR	Genes
Upregulated	hsa05144	Malaria	7	6.27E-02	HBB, HGF, ICAM1, CCL2, TGFB2, TGFB3, VCAM1
	hsa00590	Arachidonic acid metabolism	7	9.75E-02	GPX3, PTGDS, PTGIS, CYP4F2, PTGES, PLA2G16, PLA2G2C
	hsa04350	TGF-β signaling pathway	8	9.75E-02	BMP4, ID2, ID4, SMAD3, SMAD6, TGFB2, TGFB3, BAMBI
Downregulated	hsa03030	DNA replication	18	7.58E-10	MCM2, MCM4, POLA1, POLA2, POLD3, POLE, PRIM1, PRIM2, RFC2, RFC3
	hsa03460	Fanconi anemia pathway	20	7.83E-09	BRCA1, FANCD2, RAD51, TOP3A, USP1, UBE2T, POLK, FANCL, RMI1, BRIP1
	hsa05322	Systemic lupus erythematosus	30	2.47E-08	HIST1H2AI, HIST1H2AK, HIST1H2AC, HIST1H2AB, HIST1H2BE, HIST1H3I, HIST1H3G, HIST1H3J, HIST1H2AH, HIST2H4B
	hsa05034	Alcoholism	32	9.07E-06	CALM2, SLC29A1, HIST1H2BB, PKIA, HIST1H2AI, HIST1H2AK, HIST1H2AC, HIST1H2AB, HIST1H3I, HIST1H3G
	hsa03008	Ribosome biogenesis in eukaryotes	20	9.07E-06	DKC1, UTP14A, WDR3, GTPBP4, GNL3, POP5, NOP58, GNL3L, UTP6, EFTUD1
	hsa03440	Homologous recombination	12	1.78E-05	BLM, BRCA2, MRE11A, NBN, RAD51, TOP3A, XRCC2, POLD3, RAD54B, EME1
	hsa04110	Cell cycle	24	1.75E-04	CDK1, CDC6, ANAPC10, MCM2, MCM6, ORC1, TGFB1, SMC3, CCNE2, STAG1
	hsa03430	Mismatch repair	9	4.06E-04	MSH6, MLH1, MSH2, RFC2, RFC3, RFC4, RFC5, EXO1, POLD3
	hsa05203	Viral carcinogenesis	30	2.12E-03	CASP3, PIK3CA, PIK3R1, PMAIP1, RASA2, RBL1, HIST1H2BL, HIST1H2BF, HIST1H2BE, HIST1H2BI
	hsa03013	RNA transport	24	3.37E-03	EIF1AX, EIF4EBP1, EIF3E, NUP88, TPR, EIF3J, EIF4E2, NUPL1, PHAX, NDC1
	hsa00240	Pyrimidine metabolism	18	4.79E-03	CTPS1, DCK, PNP, POLA1, POLE, POLE2, PRIM1, PNP, RRM2, TYMS
	hsa03420	Nucleotide excision repair	10	1.65E-02	ERCC4, GTF2H1, MNAT1, POLE, POLE2, RFC2, RFC3, RFC4, RFC5, POLD3
	hsa03040	Spliceosome	19	2.25E-02	SRSF7, SNRPD1, SNRPG, PRPF3, SNRNP40, RBM8A, THOC1, PPIE, PPIH, SRSF10
	hsa00900	Terpenoid backbone biosynthesis	6	3.47E-02	FNTB, HMGCR, HMGCS1, IDI1, ZMPSTE24, PDSS1
	hsa03018	RNA degradation	13	3.89E-02	HSPA9, C1D, EXOSC8, EXOSC2, LSM5, CNOT10, LSM3, CNOT7, EXOSC3, PAPD5, PNPT1
	hsa05323	Rheumatoid arthritis	14	4.55E-02	ATP6V1E1, IL1B, IL6, IL8, IL11, MMP1, MMP3, TGFB1, VEGFA, ATP6V1D
	hsa05202	Transcriptional misregulation in cancer	23	6.29E-02	IGFBP3, IL6, IL8, PAX3, CDK14, PLAT, PLAU, HIST1H3I, HIST1H3G, HIST1H3J

FDR, false discovery rate; TGF-β, transforming growth factor-β.

The downregulated DEGs were significantly enriched in 17 pathways. DEGs, including MCM2, MCM4, PRIM1, POLA1 and POLA2 were enriched in the pathway of hsa03030 DNA replication (FDR=7.58E-10); DEGs, including UTP6, GTPBP4 and GNL3 were correlated with hsa03008 ribosome biogenesis in eukaryotes (FDR=9.07E-06); DEGs including RAD51 and TOP3A, were associated with hsa03440 homologous recombination (FDR=1.78E-05); DEGs, including CDC6, MCM2, CDK1, STAG1 and ANAPC10, were implicated in hsa04110 cell cycle (FDR=1.75E-04); DEGs, including PIK3R1, CASP3, HIST1H2BL and HIST1H2BF, were enriched in the hsa05203 viral carcinogenesis pathway (FDR=2.12E-03); and DEGs, including POLA1, PNP and PRIM1, were implicated in hsa00240 pyrimidine metabolism (FDR=4.79E-03; Table I).

In the pathway interaction network, the pathways of the hsa04110 cell cycle and hsa05203 viral carcinogenesis had the highest degrees of interaction, and interacted with other pathways, including hsa03030 DNA replication, hsa03008 ribosome biogenesis in eukaryotes, hsa00590 arachidonic acid metabolism, hsa00240 pyrimidine metabolism, and hsa03018 RNA degradation. In addition, the hsa04110 cell cycle pathway was found to interacted with the pathways of hsa05203 viral carcinogenesis and hsa04350 TGF-β signaling pathway (Fig. 1).

Figure 1

Pathway interaction network for the upregulated and downregulated DEGs. The size of the nodes indicates the number of pathway interactions. DEGs, differentially expressed genes; TGF-β, transforming growth factor-β.

Analysis of the cancer-associated PPI network

In total, 2,617 interaction pairs of DEGs and cancer genes were screened in the present study. The five DEGs with the highest degree were TP53, SMAD3, PIK3R1, CASP3 and CDK1, and their degrees were 104, 75, 67, 60 and 54, respectively.

With the exception of CASP3, the four residual genes (TP53, SMAD3, PIK3R1 and CDK1) were not only DEGs of NSCL/P, but were also identified as cancer genes. In the PPI network, TP53 was directly associated with E2F1, GNL3, PRIM1, PNP and POLA1; CDK1 was associated with E2F1, MCM4 and POLA1; CASP3 was directly associated with MLH1, BLM, BRCA1 and RAD51, as well as TP53 and PIK3R1. SMAD3, and TP53, were associated to certain cancer genes, including TP73, IL16, MAPK1, MAPK9 and CDK2; and SMAD3 was found to interact with ANAPC10 (Fig. 2).

Figure 2

Cancer-associated PPI network for the five DEGs with the highest degree. Yellow nodes represent the DEGs in both NSCL/P and cancer; red nodes represent the DEGs only in NSCL/P; blue nodes represent the cancer genes. The size of the nodes indicate the degree of interaction of the DEGs; PPI, protein-protein interaction; DEGs, differentially expressed genes; NSCL/P, non-syndromic cleft lip, with or without cleft palate.

Discussion

In the present study, 452 DEGs were identified to be significantly upregulated and 1,288 were found to be downregulated in the NSCL/P samples. According to the analysis of the cancer-associated PPI network, the five DEGs with the highest degrees were TP53, SMAD3, PIK3R1, CASP3 and CDK1. Among these, TP53, SMAD3, PIK3R1 and CDK1 were not only DEGs for NSCL/P, but were also associated with cancer.

TP53, encoding the p53 protein, acts as a tumor suppressor, and its loss of function is a precondition for almost all types of cancer (33). The effector functions of p53 range from arresting the cell cycle to inducing more substantial events, including senescence or apoptosis (34). A previous study demonstrated that the transcription factor p63, a homologue of p53, can transactivate IRF6 by binding to an upstream enhancer element, whose genetic variation is associated with increased susceptibility to cleft lip (35). It is also possible that p63 may be an important upstream regulator of desmosomal cell adhesion, which may contribute to the skin fragility observed in patients with cleft lip and palate (36). In addition, the L514F mutation in the sterile α-motif region of p63 can interrupt the binding of p63 to the RNA-processing protein, ABBP1, which leads to aberrant splicing of the keratinocyte growth factor receptor and inhibition of epithelial differentiation (37). The present study observed that TP53 was directly associated with GNL3, PRIM1, PNP and POLA1. GNL3, encoding guanine nucleotide binding protein-like 3, was enriched in the pathway of ribosome biogenesis in eukaryotes; PRIM1, encoding polypeptide 1 of DNA primase; PNP, encoding purine nucleoside phosphorylase; and POLA1, encoding the catalytic subunit of DNA polymerase, were enriched in the pathway of pyrimidine metabolism. PRIM1 and POLA1 were also associated with DNA replication. These pathways were all involved in the process of cell proliferation Normal palate and orofacial morphogenesis requires mesenchymal cell proliferation and differentiation, and inhibiting the progression of cell cycle between the G1 and S phases in human embryonic palatal mesenchymal cells may induce cleft palate (38). Thus, TP53 may be key in NSCL/P by modulating ribosome biogenesis, pyrimidine metabolism and DNA replication via interactions with GNL3, PRIM1, PNP and POLA1.

In addition, CDK1 was found to interact with POLA1, as well as MCM4, which were enriched in DNA replication. CDK1 encodes cyclin-dependent kinase 1, a catalytic subunit of M-phase promoting factor, which is crucial for G1/S and G2/M phase transitions in eukaryotic cell cycle (39). In addition, CDK1 and TP53 were observed to be associated with E2F1. E2F1, a master regulator of cell cycle, can promote the G1/S transition, transactivating a variety of genes involved in chromosomal DNA replication, including its own promoter (40). Increased E2F1 activity can promote tumorigenesis (41). A previous study reported that E2F1 may be involved during murine palatogenesis (42). The present study also observed that the pathway of the cell cycle interacted with the pathway of viral carcinogenesis, and PIK3R1 and CASP3 were enriched in viral carcinogenesis. PIK3R1, encoding the p85α regulatory subunit of phosphoinositide-3-kinase is known to be associated with a series of cellular processes associated with malignant behavior, including proliferation, adherence, transformation and survival (43). A PIK3R1 mutant has been identified in glioblastoma, ovarian cancer and colon cancer (44,45). In addition, CASP3, encoding a member of the cysteine-aspartic acid protease (caspase) family, is important in the extrinsic and intrinsic apoptotic pathways (46). A previous study suggested that increased expression of CASP3 is associated with tumors of the mouth (47). Thereby, CDK1, together with POLA1, MCM4, E2F1, PIK3R1 and CASP3 may not only be critical in the development of NSCL/P, but also in cancer.

SMAD3, a member of the SMAD family, was enriched in the TGF-β signaling pathway. SMAD family members are essential intracellular signaling components of the TGF-β superfamily (48). TGF-β is a cytokine, which controls the proliferation, differentiation, migration and apoptosis of several different cell types, and is important in mediating epithelial-mesenchymal transformation during the normal fusion of the lip and palate (49,50). It has been confirmed that TGF-β3 can promote fetal cleft lip repair and fusion in mouse fetuses by increasing the availability of mesenchymal cells or inducing expression of cyclin D1 (49). Also, TGF-β can promote tumor cell proliferation by stimulating the production of autocrine mitogenic factors, such as platelet-derived growth factor B (51). TGF-β can contribute to tumor invasion by inducing eithelial-mesenchymal transition (52). TGF-β can also enhance cell motility by cooperating with ERBB2, which is observed to be overexpressed in breast cancer cells (53). In addition, TGF-β can suppress immunity in patients with human glioma via decreasing the expression of the activating immunoreceptor, NKG2D, in CD8+ T cells and natural killer cells, and repressing the expression of the NKG2D ligand, MICA (54). In the present study, SMAD3 was found to interact with ANAPC10, which was enriched in the pathway of cell cycle. Thereby, SMAD3 may be an important gene in the development of NSCL/P and cancer.

In conclusion, 452 upregulated and 1,288 downregulated DEGs were identified in the present study. Five important DEGs, including TP53, CDK1 and SMAD3, may be associated with both NSCL/P and cancer. These results suggested correlation between the pathogenesis of NSCL/P and cancer, which may provide novel information for the clinical diagnosis of NSCL/P.