Identification of miRNA-mRNA-TFs Regulatory Network and Crucial Pathways Involved in Tetralogy of Fallot.

Tetralogy of Fallot (TOF) is the most common cyanotic congenital heart disease. However, its pathogenesis remains unknown. To explore key regulatory connections and crucial pathways underlying the TOF, gene or microRNA expression profile datasets of human TOF were obtained from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database. The differentially expressed mRNAs (DEmRNAs) and microRNAs (DEmiRs) between TOF and healthy groups were identified after data preprocessing, followed by Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Then, we further constructed protein-protein interaction (PPI) network and subnetwork of modules. Ultimately, to investigate the regulatory network underlying TOF, a global triple network including miRNAs, mRNAs, and transcription factors (TFs) was constructed based on the integrated data. In the present study, a total of 529 DEmRNAs, including 115 downregulated and 414 upregulated DEmRNAs, and 7 significantly upregulated DemiRs, including miR-499, miR-23b, miR-222, miR-1275, miR-93, miR-155, and miR-187, were found between TOF and control groups. Furthermore, 22 hub genes ranked by top 5% genes with high connectivity and six TFs, including SRF, CNOT4, SIX6, SRRM3, NELFA, and ONECUT3, were identified and might play crucial roles in the molecular pathogenesis of TOF. Additionally, an miRNA-mRNA-TF co-regulatory network was established and indicated ubiquitin-mediated proteolysis, energy metabolism associated pathways, neurodevelopmental disorder associated pathways, and ribosomes might be involved in the pathogenesis of TOF. The current research provides a comprehensive perspective of regulatory mechanism networks underlying TOF and also identifies potential molecule targets of genetic counseling and prenatal diagnosis for TOF.

Introduction

Cardiac development is an extremely complex process involving the integration of multiple cell lineages into the three-dimensional organ and its seamless connection to the vascular system (Liu and Olson, 2010). Proper embryological heart development and function are precisely controlled by an evolutionarily conserved network of transcription factors (TFs), which are triggered by upstream signaling systems involved in cardiac myogenesis, morphogenesis, and contractility (Olson, 2006; Liu and Olson, 2010). Abnormalities in the networks during the cardiac formation process might lead to the failure of proper cardiac differentiation, migration, and apoptosis, which finally results in congenital heart disease (CHD).

Tetralogy of Fallot (TOF) is the most common cyanotic CHD, occurring in about 3 of every 10,000 live births and accounting for 7–10% of CHD (Bailliard and Anderson, 2009). TOF is caused by the non-uniform separation between bulbos arteriosus and truncus arteriosus during the early stage of embryonic development, resulting in malformations such as obstruction of the right ventricular outflow tract, ventricular septal defect (VSD), override of the ventricular septum by the aortic root, and right ventricular hypertrophy (Bailliard and Anderson, 2009; Wang et al., 2018). The clinical, pathological, phenotypical, and biologic heterogeneity of TOF suggests the complex pathogenesis. The precise molecular pathogenesis of TOF has not been completely identified. However, the latest research is advancing the understanding of the molecular mechanisms of TOF. Today, most investigations into the pathogenesis of TOF have focused on mutations in specific protein-coding genes, such as GATA4, NKX-2.5, JAG1, FOXC2, TBX5, TBX1, etc. (Lalani and Belmont, 2014). In addition, microdeletions of chromosome 22q11.2 and several other copy number variations, including15q11.2, 1q21.1, 2p13.3, 3p25.1, and 16p13.11, are found implicated in the pathological mechanisms of TOF (Lalani and Belmont, 2014).

Recently, several studies have revealed that microRNAs (miRNAs) play crucial roles in cardiac signaling and transcriptional pathway modulating multiple cardiac development, function, and disease (Liu and Olson, 2010). Cardiac development depends on the proper spatiotemporal expression of particular miRNAs. miR-1 was the first reported miRNA to be implicated in cardiac development, which targets DLL1 and HAND2 and further promotes the differentiation of embryonic stem cells into cardiac lineage (O’Brien et al., 2012). MiR-133a-1 and miR-133a-2 double-knockout mice have abnormal expression of SRF and CCND2, leading to late embryonic and neonatal lethality due to chamber dilatation and VSD (Liu et al., 2008). Deletion of either miR-106a-363 or miR-106b-25 combined with the miR-17∼92 null allele can result in embryonic lethality, accompanied by severe VSD, atrial septal defects, and thin-walled myocardium (Ventura et al., 2008). Taken together, these studies revealed that miRNAs play vital roles in vertebrate heart development. However, a comprehensive profiling of miRNA regulatory networks in TOF has not been reported to date.

The molecular mechanism and critical regulators involved in TOF are poorly deciphered. In the present study, to investigate the regulatory networks and key pathways underlying TOF, we constructed the protein–protein interaction (PPI) and an integrative miRNA–mRNA–TF regulatory network to further elucidate the key mechanisms underlying TOF.

Materials and Methods

Collection of Datasets

All of the TOF datasets were obtained from the Gene Expression Omnibus (GEO) (Clough and Barrett, 2016) with the following criteria: (1) mRNA or miRNA expression profiles by microarray or RNA-seq could be available; and (2) TOF patients and healthy controls were estimated. Consequently, a total of four microarray datasets [GSE35776 (O’Brien et al., 2012), GSE26125 (Bittel et al., 2011), GSE40128 (Zhang et al., 2013), and GSE35490 (Bittel et al., 2014)] and one RNA-seq dataset [GSE36761 (Grunert et al., 2014)] were included in this study, in which GSE40128 and GSE35490 were used for its miRNA expression profile. A total of 107 samples (32 controls and 75 cases) were analyzed in this study; details are provided in Table 1.

TABLE 1

Microarray and RNA-seq datasets used in this study and their experimental design.

	GEO accession number	Sample size		Platform
		Controls	Patients
1	GSE35776	8	16	(HuEx-1_0-st) Affymetrix Human Exon 1.0 ST Array
2	GSE36761	7	22	Illumina Genome Analyzer
3	GSE26125	5	16	CodeLink Human Whole Genome Bioarray
4	GSE40128	3	5	(miRNA-1) Affymetrix Multispecies miRNA-1 Array
5	GSE35490	8	16	(miRNA-1) Affymetrix Multispecies miRNA-1 Array

Identification of Differentially Expressed mRNAs and miRNAs

For microarray data, background correction, quantile normalization, and determination of expression levels were conducted using R/Affy (Gautier et al., 2004) or Oligo (Carvalho and Irizarry, 2010) packages. Differential expression analysis of mRNAs or miRNAs between the TOF and control groups was conducted using the R/Limma package (Diboun et al., 2006).

For RNA sequencing dataset (GSE36761), adapter sequences of raw reads were trimmed by Trimmomatic (Bolger et al., 2014). STAR was utilized to align the short reads to the reference human genome (hg38) (Dobin et al., 2013). Accepted mapped reads were then summarized at gene level counts using featureCounts software (Liao et al., 2014). Differential gene expression analysis between the TOF and control groups was conducted using R package DESeq2 (Love et al., 2014).

The adjusted P < 0.05 and |log Fold Change (FC)| > 2 were set as the cut-off criteria for differentially expressed mRNAs (DEmRNAs), as well as differentially expressed miRNAs (DEmiRs) with thresholds of |log FC| > 1.5 and adjusted P < 0.05.

miRNA–Target Regulatory Network Analysis

miRWalk software is a comprehensive database that provides predicted and validated miRNA binding sites of known genes of human, rat, cow, mouse, and dog (Dweep and Gretz, 2015). Experimentally verified miRNA–gene regulatory pairs were considered as targets of DEmiRs. Then, we extracted the genes as significant DEmRNAs by intersecting the target genes of the DEmiRs and the overlapping genes of the DEmRNAs. The miRNA–target gene regulatory network was constructed using the Cytoscape software (Shannon et al., 2003).

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analysis

To further understand the gene function in TOF, Gene Ontology (GO) analysis (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (Kanehisa and Goto, 2000) of the DEmRNAs were performed using Metascape (Zhou et al., 2019), with P < 0.05 as the cut-off criterion.

PPI Network Construction and Identification of Hub Genes

The STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database provides a critical assessment and integration of PPI, including direct (physical) as well as indirect (functional) associations (Szklarczyk et al., 2015). The PPI networks of DEmRNAs were identified using the STRING database, and PPI score threshold was set as 0.4 to obtain interactions with higher confidence (>0.4). As nodes with a higher degree of connectivity contributing more to the stability of the network, the top 5% DEmRNAs ranked by the degree connectivity were defined as hub genes by using the plugin CentiScaPe (Scardoni et al., 2009). The Molecular Complex Detection (MCODE) (Bader and Hogue, 2003) plugin was applied to identify significant modules with node score cutoff = 0.2, K-Core = 4, max depth = 100 and degree cut-off = 4 set as the cut-off criteria. The PPI network was visualized by Cytoscape software (Shannon et al., 2003).

miRNAs–mRNAs–TFs Regulatory Network Analysis

MicroRNAs can regulate gene expression at the posttranscriptional stage (Slezak-Prochazka et al., 2010), whereas TFs can activate or repress transcription at a pretranscriptional stage (Zaret and Carroll, 2011). The miRNA–target regulatory network was visualized to depict interactions between miRNAs and their potential targets using Cytoscape 3.7.1.

The plugin iRegulon (Janky et al., 2014) in Cytoscape software, which included multiple human databases of the TF–target pairs such as Transfac and Encode, was used to identify TFs regulating DEmRNAs in the miRNA–target regulatory network. In this study, the parameters for TF motif enrichment analysis were set as the minimum identity between orthologous genes ≥0.05, maximum false discovery rate (FDR) on motif similarity ≤0.001, and TF motifs with normalized enrichment score (NES) >3 were considered as the threshold value for the selection of potential associations. TransmiR is a TF–miRNA regulation database, through which regulatory relations between TFs and miRNAs can be identified (Tong et al., 2019). Ultimately, the miRNA–mRNA—-TF regulatory network was constructed utilizing Cytoscape 3.7.1.

Results

Identification of DEmRNAs and DE miRNAs (DEmiRs) for TOF

Compared with the normal samples, following quality control and data normalization, a total of 1735 DEmRNAs including 1315 upregulated and 420 downregulated DEmRNAs were identified based on the criteria. miRNA expression analysis showed that a total of seven upregulated miRNAs were identified as significant DEmiRs: miR-155, miR-222, miR-23b, miR-499-3p, miR-1275, miR-93, and miR-187. The volcano plots of the DEmRNAs or DEmiRs are presented in Supplementary Figure S1.

We obtained 6261 targeted genes of the 7 significant DEmiRs from the miRWalk database. Four hundred and fourteen genes were found as the intersection genes between these target genes and the overall upregulated DEmRNAs. Therefore, a total of 529 DEmRNAs, including the 414 upregulated and 115 downregulated genes, were identified as the final sets of significant DEmRNAs (Supplementary Table S1).

GO and KEGG Pathway Enrichment Analysis of DEmRNAs

Functional and pathway enrichment analyses were applied to explore the pathogenesis in TOF. The DEmRNAs enriched according to the P-value are visualized in Figure 1A (Supplementary Table S2). The significantly top five GO terms in DEmRNAs were striated muscle cell differentiation, mitochondrial matrix, cofactor metabolic process, myofibril assembly, and generation of precursor metabolites and energy. The significant top five PPI pathways of DEmRNAs were ubiquitin-mediated proteolysis, non-alcoholic fatty liver disease, Parkinson’s disease, Huntington’s disease, and citrate cycle (Figure 1B and Supplementary Table S3).

FIGURE 1

Functional enrichment analyses of differentially expressed genes (DEmRNAs). (A) Enriched ontology clusters. Each term is represented by a circle node, where its size is proportional to the number of input genes falling into that term, and its color represents its cluster identity (i.e., nodes of the same color belong to the same cluster). Terms with a similarity score >0.3 are linked by an edge (the thickness of the edge represents the similarity score). (B) The top 20 KEGG enrichment analysis for DEmRNAs. The numbers on the x-axis are the names of pathways. The numbers on the y-axis are gene counts.

PPI Network Construction and Module Analysis

To investigate the relationship of the DEmRNAs, a PPI network was constructed utilizing Cytoscape software based on the STRING database. In total, 448 nodes and 1693 PPI relationships were obtained (Figure 2). Twenty-two genes ranked by top 5% genes with high connectivity were defined as hub genes for TOF, all of which were upregulated (Table 2 and Supplementary Table S4). The 22 hub genes interact directly with 337 edges. GO analysis results showed that hub genes were significantly enriched in ubiquitin–protein transferase activity and ubiquitin protein ligase activity in the molecular function ontology, cytosolic part and ubiquitin ligase complex in cellular components ontology, and protein polyubiquitination and cellular protein catabolic process in biological processes ontology. The results of KEGG analysis showed that hub genes were significantly enriched in ubiquitin-mediated proteolysis and protein processing in endoplasmic reticulum (ER) and ribosomes (Figure 3 and Supplementary Table S5). These hub genes might play a crucial role in TOF and should be studied further.

FIGURE 2

The protein–protein interaction (PPI) network of differentially expressed genes (DEmRNAs). Larger circles with red colors represent hub genes, and smaller hexagons with green colors represent non-hub genes.

TABLE 2

Hub genes identified in the protein–protein interaction network.

Gene	Degree	Log FC	adj.P.Val	Regulation
RPS27A	70	2.2035714	0.0008567	Up
HSPA8	46	2.6819216	0.0001913	Up
SKP1	38	3.0387064	0.0004551	Up
HSPA9	36	3.0767256	0.0001602	Up
RPS2	36	2.7716774	0.0008141	Up
RPL8	33	2.5430817	0.0005055	Up
SOD2	32	3.2653954	7.139E−05	Up
FBXO32	31	4.3608082	3.498E−06	Up
TCP1	31	3.2699954	0.0001236	Up
UBE2N	29	2.2827098	0.0012278	Up
UBE3A	28	2.0175226	0.0059344	Up
PSMC2	28	2.1476242	1.475E−13	Up
FBXW11	27	2.1880644	0.000518	Up
DLD	27	4.1483123	4.553E−05	Up
CCT5	27	2.5370472	0.0029727	Up
RPS24	27	3.2131167	0.000303	Up
SIAH2	26	2.4161537	0.0001602	Up
UBE2D3	26	2.3856558	0.0007899	Up
VDAC1	26	3.971706	4.006E−06	Up
RPL7A	26	2.668917	0.00179	Up
ZBTB16	25	2.4423304	0.0034434	Up
UBE2H	24	2.427039	0.000537	Up
UBR1	24	2.065815	0.0006099	Up
RCHY1	24	2.8119357	6.861E−05	Up

FIGURE 3

Functional enrichment analyses of the hub genes. The numbers on the x-axis are the names of pathways or GO terms. The numbers on the y-axis are gene counts. BP, biological progress; CC, cellular component; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

A total of three modules from the PPI network satisfied the criteria of an MCODE computed K-score = 4 and degree cut-off = 4 (Figure 4). Fifty-nine MCODE-related genes were identified and which contained the 17 hub genes. Module 1 (Figure 4A) included 21 nodes and 210 PPI pairs; KEGG analysis showed that these genes were significantly enriched in ubiquitin-mediated proteolysis; and GO analysis showed that these genes were mainly associated with protein polyubiquitination, proteasomal protein catabolic process, contractile fiber part, sarcomere, and ubiquitin-protein transferase activity (Figure 4D). Module 2 included 27 nodes and 84 PPI pairs (Figure 4B), and KEGG analysis showed that these genes were enriched in energy metabolism associated pathways (oxidative phosphorylation, citrate cycle, and carbon metabolism), and neurodevelopmental disorders (Parkinson’s disease, Huntington’s disease, and Alzheimer’s disease), and GO terms were primarily in cellular respiration, organelle localization, oxidoreductase complex, mitochondrial protein complex, oxidoreductase activity, and ATPase activity (Figure 4E). Module 3 included 11 nodes and 32 PPI pairs (Figure 4C); KEGG analysis showed that these genes were significantly enriched in ribosomes, and GO analysis revealed that module 3 related genes were mainly associated with establishment of protein localization to organelles, translation initiation, cytosolic part chaperone complex, structural constituents of ribosomes, and unfolded protein binding (Figure 4F). The detailed GO and KEGG pathway analysis information of three modules is shown in Supplementary Table S6.

FIGURE 4

Module analysis of differentially expressed genes (DEmRNAs). (A) Module 1. (B) Module 2. (C) Module 3. (D) Functional enrichment analyses of Module 1. (E) Functional enrichment analyses of Module 2. (F) Functional enrichment analyses of Module 3. Larger circles with red colors represent hub genes, smaller hexagons with green colors represent non-hub genes, and black arrows represent seed genes. BP, biological progress; CC, cellular component; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

miRNA–TF–Target Regulatory Network Analysis

To further understand the regulatory relationship between DEmiRs and DEmRNAs, the miRNA–target regulation network was constructed as shown in Figure 5. miR-222 had the highest connectivity with 529 target DEmRNAs. Moreover, several miRNAs were predicted to have common targets, for example, miR222 and miR93 have 139 common target genes, and miR222 and miR23b have 90 common target genes (Figure 5).

FIGURE 5

miRNA–target regulatory networks of differentially expressed genes (DEmRNAs). The pink V-shapes represent miRNAs, red circles represent upregulated DEmRNAs, and green circles represent downregulated DEmRNAs.

Six TFs were predicted to regulate DEmRNAs by iRegulon analysis, including SRF, CNOT4, SIX6, SRRM3, NELFA, and ONECUT3. To determine the regulatory connections between TFs and DEmRNAs, the TF–target regulation network of 302 nodes and 556 interaction pairs was constructed (Figure 6), SRF targeted 254 genes, CNOT4 targeted 60 genes, SIX6 targeted 84 genes, SRRM3 targeted 44 genes, NELFA targeted 46 genes, ONECUT3 targeted 68 genes. The top 20 nodes with the highest degree are listed in Table 3, including HSPA8, SKP1, HSPA9, FBXO32, UBE3A, FBXW11, VDAC1, UBE2D3, SIAH2, ZBTB16, UBE2H, HIF1A, UQCRC2, MDH2, SQSTM1, FBXO22, ASB11, HUWE1, HECW2, and RAB1A, which were coregulated by the three TFs.

FIGURE 6

The TF–target regulatory network of differentially expressed genes (DEmRNAs). Orange triangles represent TFs, red circles represent upregulated genes, and green circles represent downregulated genes. TFs, transcription factors.

TABLE 3

Top 20 nodes with higher degree in transcription factor–target regulatory network.

Node name	Regulation	Degree	TF
HSPA8	Up	46	SIX6
SKP1	Up	38	SIX6
HSPA9	Up	36	SRF
FBXO32	Up	31	SRF
UBE3A	Up	28	SRF
FBXW11	Up	27	SRF
VDAC1	Up	26	SRF
UBE2D3	Up	26	SRF
SIAH2	Up	26	SRF
ZBTB16	Up	25	SIX6
UBE2H	Up	24	SRF
HIF1A	Down	23	SRF
UQCRC2	Up	23	SIX6
MDH2	Up	23	SRF
SQSTM1	Up	23	SIX6
FBXO22	Up	23	SRF
ASB11	Up	23	NELFA
HUWE1	Up	22	SIX6
HECW2	Up	22	SRF
RAB1A	Up	21	SIX6

TransmiR analysis revealed that SRF regulated miR-23b and miR-93; miR222 and miR-23b are regulated by MITF, but also target MITF; BCL6 is the target of miR-93 and miR-155, which also regulates miR-93 and miR-155; ZBTB16 could repress miR-222, while FOSL1 could activate miR-222. Finally, the miRNA–mRNA–TF regulatory network was established with 7 miRNAs, 529 DEmRNAs, and 6 TFs through Cytoscape 3.7.1 (Figure 7).

FIGURE 7

miRNA–mRNA–TF regulatory network of differentially expressed genes (DEmRNAs). The pink diamonds represent miRNAs, orange triangles represent TFs, red circles represent upregulated genes, and green circles represent downregulated genes. TFs, transcription factors.

Discussion

Heart development is an extremely complex process involving networks of various signaling pathways, the downstream TFs, protein-coding genes, and posttranscriptional regulation for cardiac myogenesis, morphogenesis, and contractility. Owing to the lack of specific molecule mechanisms, genetic counseling and prenatal diagnosis of TOF have been a major challenge. Therefore, it is imperative to explore the specific molecular mechanism for TOF, which will be valuable for risk assessment and even prognosis prediction. Recently, certain dysregulated mRNAs, miRNAs, and TFs have been reported in TOF individually. Notably, transcriptional and posttranscriptional regulation plays important roles in TOF development.

In our study, an integrative analysis of microarray and RNA-sequencing datasets was conducted to identify the regulatory network underlying TOF. A total of 529 DEmRNAs, 7 DEmiRs, and 6 TFs were identified as significant elements for the TOF. Based on 529 DEmRNAs, we obtained 22 hub genes and three significant modules from PPI networks using MCODE analysis. Eleven hub genes were significantly enriched in protein polyubiquitination, including SKP1, FBXO32, UBE3A, PSMC2, RPS27A, SIAH2, UBE2D3, UBE2N, FBXW11, DLD, and ZBTB16. The ubiquitin–proteasome system (UPS) is important for the regulation of cellular protein degradation and synthesis (Al-Hassnan et al., 2016). Impaired UPS has been involved in the etiology of human cardiovascular diseases, such as cardiomyopathy (CMP), heart failure, and cardiac ischemia. Skp1/Cul1/F-box (SCF) ubiquitin E3 ligase complex plays a crucial role in ubiquitination of cardiac proteins (Jang et al., 2011). FBXO32, which encodes a muscle-specific ubiquitin ligase, has been implicated in the pathogenesis of CMP through the UPS (Al-Hassnan et al., 2016). UBE3A is one of the important members of UPS; Ube3a mutant mice had deficits in Ca2+/calmodulin-dependent kinase II (CaMKII), while CaMKII played a vital role in cardiac hypertrophy (Cheng et al., 2019). ZBTB16, a Kruppel-like zinc finger protein, highly expressed in the heart, is an important TF that regulates cardiac hypertrophy through the Ang II receptor 2 in response to angiotensin II (Wang et al., 2012). GO function annotation also indicates that many hub genes (HSPA8, PSMC2, RPL7A, RPL8, RPS2, RPS24, RPS27A, HSPA9, TCP1, CCT5) are significantly enriched in regulation of transcription, posttranscriptional regulation and translational regulation, including mRNA catabolic process, RNA catabolic process, ribosomes, translational initiation, nuclear-transcribed mRNA catabolic process, regulation of mRNA stability, DNA biosynthetic process, and rRNA processing. Elevated rates of protein synthesis have been demonstrated to contribute to cardiomyocyte hypertrophy. Accelerated rates of protein synthesis can be achieved by altering the translational efficiency of existing ribosomes and/or by accelerating the rate of synthesis of new ribosomes, indicating a crucial role for ribosome biogenesis in the regulation of cardiac growth (Brandenburger et al., 2003). Abnormal expression of ribosomes can affect the development and function of the cardiomyocyte and may lead to right ventricular hypertrophy of TOF. Members of the heat shock protein 70 family (HSPA8 and HSPA9) have also been proven to be essential to maintain cardiomyocyte protein homeostasis (Fang et al., 2017). Underexpression of HSPA8 and HSPA9 will disrupt cardiomyocyte homeostasis and interfere with critical cellular functions by small aggregates of proteins, which finally contributes to heart failure and hereditary heart diseases. However, the biochemical property and the molecular function of other hub genes have not been fully elucidated, and further functional studies are crucial to better understand the mechanism underlying the pathogenesis of TOF.

Gene Ontology analysis of DEmRNAs showed that DEmRNAs are mainly involved in striated muscle cell differentiation, mitochondrial matrix, cofactor metabolic process, myofibril assembly, and generation of precursor metabolites and energy, and these biological processes mainly contribute to cardiac cell division and differentiation. KEGG pathway analysis demonstrated abnormal molecular expression of several pathways may contribute to the pathogenesis of TOF, such as ubiquitin-mediated proteolysis pathways, energy metabolism associated pathways (citrate cycle, cysteine and methionine metabolism, glycolysis/gluconeogenesis, and PPAR signaling pathway), and pathways associated with neurodevelopmental disorders (Parkinson’s disease, Huntington’s disease, and Alzheimer’s disease). Ubiquitin-mediated proteolysis is essential to maintain the integrity of cellular proteins that make up sarcomeres, mitochondria, ER, and membranes and therefore is critical to cardiomyocyte survival and functioning (Pagan et al., 2013; Li et al., 2018). Insufficient ubiquitin-mediated proteolysis might be a pathogenic factor leading to cardiac pathogenesis (Li et al., 2018). Myocardial energy metabolism is an important determinant of cardiac structure and function (Finck, 2007). Energy metabolism associated pathways affect the normal cellular biological processes, including cell survival/death, differentiation/proliferation, and DNA repair (Dhar-Chowdhury et al., 2007; Hanahan and Weinberg, 2011; Tran and Wang, 2019). Pathological alterations of energy metabolism pathways could result in impaired signaling transduction, perturbed ion and redox homeostasis, and contractile dysfunction. Future studies in the specific role of energy metabolism will shed new light on the pathogenesis of TOF. Fifteen genes, including SLC25A4, ATP5F1A, ATP5F1B, NDUFA10, NDUFB5, POLR2E, REST, SDHD, SOD2, UQCRC2, LPL, VDAC1, VDAC2, UQCR11, and GRIN2A, were also enriched in neurodevelopmental disorders associated pathways. New research shows neurodevelopmental disorders occur in 10% of all children with CHD and in 50% of children with severe CHD (Marino et al., 2012), and de novo mutations account for at least 20% of CHD patients with associated extracardiac or neurodevelopmental abnormalities (Homsy et al., 2015), indicating a striking shared genetic etiology between TOF and neurodevelopmental disorders.

Among the six TFs predicted in the present study, SRF has been previously reported to contribute to differentiation, maturation, and homeostasis of cardiomyocytes, Srf knockout mice show severe defects vascular smooth muscle cells lacking organizing actin/intermediate filament bundles and die at the gastrulation stage (Schlesinger et al., 2011). CNOT4 functions as a ubiquitin-protein ligase (E3) (Deshaies and Joazeiro, 2009), interacting with a subset of ubiquitin−conjugating enzymes (E2s). SIX6 can mediate the repression of p27/Kip1, an inhibitor of the cell cycle G1/S transition, thereby regulating cellular organ−specific proliferation (Conte et al., 2005). SRRM3 is a alternative splicing regulator required for motor coordination (Nakano et al., 2019), potentially related to myocardial systolic motion. NELFA, a member of the NELF, cause promoter proximal pausing on the HSP70 gene (Wu et al., 2003), while HSPA8 and HSPA9 are identified hub genes of higher degrees. ONECUT3 is essential in multiple tissues for proper specification and development of unique cell types (Kropp and Gannon, 2016). However, whether these TFs are implicated in the pathogenesis of TOF can be investigated in future studies.

MicroRNAs, one of the most important endogenous epigenetic factors, inhibit posttranscriptional gene expression of target genes. Recently, several studies have revealed miRNAs as governing gene expression during cardiac development (Small and Olson, 2011). Seven upregulated DEmiRs were identified in the present study, including miR-499, miR-155, miR-23b, miR-222, miR-1275, miR-93, and miR-187. Among of them, miR-499, miR-155, miR-23b, and miR-222 have been demonstrated as participating in cardiac development and function. miR-499 regulates the mitochondrial fission machinery in cardiomyocytes by suppressing calcineurin-mediated dephosphorylation of DRP1 (Wang et al., 2011). Similarly, miR-1275 was found to correlate with decreased mitochondrial membrane potential (Nymark et al., 2015). miR-155 is a critical microRNA that represses MEF2A expression (Seok et al., 2011), and Mef2a knockout mice exhibit dilated right ventricle, myofibrillar fragmentation, mitochondrial disorganization, and cardiac sudden death (Naya et al., 2002). miR-222 is upregulated in ventricular outflow tract tissues from infants with non-syndromic TOF, and miR-222 is involved in proliferation and differentiation of cardiomyocytes in in vitro experiments (Zhang et al., 2013). miR-23b was found upregulated in cyanotic CHD patients, and the upregulated miR-23b increases cardiomyocyte apoptosis and reduces cell growth under hypoxic conditions (Li et al., 2017). miR-93 was identified as highly expressed in the fetal mouse heart at four key time points (embryonic day E12.5, E14.5, E16.5, and E18.5) (Cao et al., 2012). miR-187 and its target genes Itpkc, Tbl1xr1, and Lrrc59 were reported to play regulatory roles in cardiomyocyte apoptosis and cardiac inflammation (Ektesabi, 2019). Therefore, we speculate that these miRNAs mainly contribute to abnormal cardiac cell division, differentiation, and apoptosis for TOF pathogenesis.

Conclusion

In summary, we performed an integrated analysis, based on four microarray datasets and one RNA-seq dataset, to identify specific TFs, miRNAs, and mRNAs that may play a pivotal role in TOF. Moreover, we established a TF–miRNA–mRNA network in TOF. However, there are certain limitations in our present study. The sample number of TOF specimens from patients was relatively small, which may reduce the credibility of the miRNA–mRNA–TF coregulatory network analysis. In addition, the present study lacked further experiments to verify the analysis results as a solid foundation. All in all, the current research provides a comprehensive perspective of regulatory mechanism networks underlying TOF and also identifies potential molecule targets of genetic counseling and prenatal diagnosis.

Data Availability Statement

The datasets generated for this study can be found in the Gene Expression Omnibus (GEO): GSE35776, GSE36761, GSE26125, GSE40128, GSE35490.

Author Contributions

FL conceived and designed the experiments. GY, BZ, BW, and QF analyzed and interpreted the data. GY and BZ wrote the manuscript. All the authors reviewed the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.