Microarray analysis reveals key genes and pathways in Tetralogy of Fallot.

The aim of the present study was to identify key genes that may be involved in the pathogenesis of Tetralogy of Fallot (TOF) using bioinformatics methods. The GSE26125 microarray dataset, which includes cardiovascular tissue samples derived from 16 children with TOF and five healthy age‑matched control infants, was downloaded from the Gene Expression Omnibus database. Differential expression analysis was performed between TOF and control samples to identify differentially expressed genes (DEGs) using Student's t‑test, and the R/limma package, with a log2 fold‑change of >2 and a false discovery rate of <0.01 set as thresholds. The biological functions of DEGs were analyzed using the ToppGene database. The ReactomeFIViz application was used to construct functional interaction (FI) networks, and the genes in each module were subjected to pathway enrichment analysis. The iRegulon plugin was used to identify transcription factors predicted to regulate the DEGs in the FI network, and the gene‑transcription factor pairs were then visualized using Cytoscape software. A total of 878 DEGs were identified, including 848 upregulated genes and 30 downregulated genes. The gene FI network contained seven function modules, which were all comprised of upregulated genes. Genes enriched in Module 1 were enriched in the following three neurological disorder‑associated signaling pathways: Parkinson's disease, Alzheimer's disease and Huntington's disease. Genes in Modules 0, 3 and 5 were dominantly enriched in pathways associated with ribosomes and protein translation. The Xbox binding protein 1 transcription factor was demonstrated to be involved in the regulation of genes encoding the subunits of cytoplasmic and mitochondrial ribosomes, as well as genes involved in neurodegenerative disorders. Therefore, dysfunction of genes involved in signaling pathways associated with neurodegenerative disorders, ribosome function and protein translation may contribute to the pathogenesis of TOF.

Introduction

Tetralogy of Fallot (TOF) is a congenital heart disease, with an incidence rate estimated at 5–7/10,000 live births worldwide (1). TOF is characterized by ventricular septal defects, sub-pulmonary and pulmonary stenosis, an over-riding aorta and right ventricular hypertrophy (2). At present, the molecular mechanisms underlying the pathogenesis of TOF remain poorly understood. In past decades, the association between the pathogenesis of this disease and mutations in specific genes, including GATA binding protein 4 (3,4), GATA binding protein 6 (5), zinc finger protein, FOG family member 2 (4,6) and jagged 1 (7) have been reported; although the function of these genes remains controversial among different studies. In addition, a deletion in chromosome 22q11 (8,9) and copy number variations in various chromosomes, such as duplications 1q21.1 and micro deletions in 14q23 (10), have been implicated in the pathogenesis of TOF.

Through the use of genome-wide gene expression microarrays, Bittel et al (11) revealed that the majority of abnormally expressed genes in cardiac tissue samples derived from patients with TOF were involved in compensatory functions, including hypertrophy, cardiac fibrosis and cardiac dilation, and the expression of genes involved in the WNT and Notch signaling pathways were suppressed. Using bioinformatics methods, the present study employed the microarray results submitted by Bittel et al (11) to further identify key genes that may be involved in the pathogenesis of TOF.

Data and methods

Microarray data

The GSE26125 expression profile dataset was downloaded from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) (12). The mRNA used for array hybridization was extracted from the cardiovascular tissue samples of 16 children with TOF and five healthy age-matched control infants (11). The CodeLink Human Whole Genome Bioarray (Applied Microarrays, Inc., Tempe, AZ, USA), which contains >54,000 probes, was employed to analyze the samples (11). In addition, the GPL11329 CodeLink Human Whole Genome Bioarray (Applied Microarrays, Inc.) annotation platform was used.

Data preprocessing and identification of differentially expressed genes (DEGs)

The downloaded data were subjected to background correction, quantile normalization and probe summarization using the robust multi-array average method (13) with the R/Affy package version 3.2.2 in Bioconductor release 3.2 (14). Differential expression analysis of genes between the TOF and control groups was then performed using the Student's t-test with the R/limma package version 3.2.2 in Bioconductor release 3.2 (15). Genes with a log2 fold-change (FC) of ≥2 and a false discovery rate (FDR) of <0.01 were considered to be DEGs. The samples were subsequently clustered based on the identified DEGs using the pvclust R package (16) by calculating the approximate unbiased P-values.

Functional annotation and pathway enrichment analyses of DEGs

The biological function of the identified DEGs was determined using the following tools: ToppGene (https://toppgene.cchmc.org/) (17), which is a one-stop portal for functional enrichment analysis of gene lists based on the Gene Ontology (GO) database (18); the BioSystems database (19); the BIOCYC database (http://BioCyc.org/) (20); the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (21); the REACTOME pathway database (22,23); WikiPathways September 2015 release (http://www.wikipathways.org/index.php/Wiki Pathways), GenMAPP version 2.1 (Gladstone Institutes, San Francisco, CA, USA), the MSigDB C2 database version 5.0 (http://software.broadinstitute.org/gsea/msigdb), which integrated BioCarta, (Sigma-Aldrich and Signaling Gateway; http://www.biocarta.com/), the PANTHER database version 1.4 (http://pantherdb.org/downloads/index.jsp) (24), Pathway Ontology (http://www.obofoundry.org/ontology/pw.html) (25) and the Small Molecule Pathway Database (SMPDB) version 2.0 (http://www.smpdb.ca) (26,27). FDR values of <0.05 and a gene count of ≥2 were set as the threshold values. GO categories were classified into the following terms: Biological process, molecular function and cellular component.

Gene functional interaction (FI) network and pathway enrichment analyses

ReactomeFIViz (28) is a Cytoscape software version 3.2.0 (29) application that allows researchers to identify network and signaling pathway patterns of interest, search for gene signatures from gene expression data sets, reveal signaling pathways significantly enriched by a list of genes, as well as integrate multiple genomic data types into a pathway context using probabilistic graphical models. In the present study, an FI network was constructed by merging interactions extracted from curated human pathways, with interactions predicted using a machine learning approach. The correlations among genes involved in the same FI were calculated and then used as weights for edges in the FI network. Then, a Monte Carlo Localization graph clustering algorithm was applied to the weighted FI network to generate a sub-network for a list of selected network modules, based on module size and average correlation. Each parameter was set a default value during the analysis as follows: Size of MCL (Markov Cluster Alorithm) clustering result=7; inflation parameter for MCL=5.0; and average correlation=2.5. The gene FI networks were visualized using Cytoscape software (29). Pathway enrichment analysis of each function module was subsequently performed to identify the signaling pathways enriched by genes in each module, with FDR values of <0.05.

Construction of gene-transcription factor regulation networks

The iRegulon plugin (30) in Cytoscape software version 3.2.0, that is associated with and contains information form the integrated databases TRANSFAC, JASPAR, ENCODE (https://genome.ucsc.edu/ENCODE/), SwissRegulon (http://swissregulon.unibas.ch/sr/), and HOMER (http://homer.ucsd.edu/homer/motif/motifDatabase.html), was used to identify transcription factors predicted to regulate the DEGs in the FI network according to the following parameters: A minimum identity between orthologous genes equal to 0.05, and a maximum FDR value of motif similarity equal to 0.001. A larger normalized enrichment score (NES) indicates a higher reliability, and an NES of >3.5 was set as the threshold. The gene-transcription factor pairs were then visualized using Cytoscape software.

Results

Identification of DEGs

Using log2FC values of >2 and FDR values of <0.01 as thresholds, a total of 878 DEGs were identified, including 848 upregulated genes and 30 downregulated genes. Among the 21 samples, 20 were correctly divided into the TOF or control groups by clustering analysis based on the identified DEGs, with an accuracy rate of ~95% (Fig. 1). One sample (TOF_15) was not divided into the TOF group by clustering analysis. Nevertheless, this indicated a relatively good performance of the clustering analysis.

Figure 1.

A cluster dendrogram of AU/BP values demonstrating the division of TOF and control samples based on the identified differentially expressed genes. TOF, Tetralogy of Fallot; AU, approximately unbiased value; BP, bootstrap probability value.

GO functional annotation and pathway enrichment analyses of DEGs

According to the different databases, the upregulated genes were enriched in multiple signaling pathways, including oxidative phosphorylation, the electron transport chain, Huntington's disease, Parkinson's disease, metabolism, the citric acid cycle and respiratory electron transport (Fig. 2). By contrast, the downregulated genes were not significantly enriched in any signaling pathway; however, they were enriched for several GO terms, such as hemoglobin complex, activation of nuclear factor-κB-inducing kinase activity and oxygen transport (Fig. 2).

Figure 2.

GO functional annotation and pathway enrichment analysis of the identified (A) upregulated and (B) downregulated differentially expressed genes between TOF and control group samples. GO, gene ontology; TOF, Tetralogy of Fallot; FDR, false discovery rate; NF-κB, nuclear factor-κB.

Gene FI network and pathway enrichment analyses

The gene FI network contained 7 function modules, each consisting of upregulated genes (Fig. 3). Genes in Modules 0, 1, 2, 3, 5 and 6 were further enriched in ≥ 1 signaling pathways (Table I). Notably, genes enriched in Module 1, including nictotinamide adenine dinucleotide (NADH):ubiquinone oxidoreductase (NDUF) subunit B5, NDUFB8, NDUFA6, NDUFB9, cytochrome c1, NDUFAB1, NDUFC2, NDUFC1, NDUFB1, NDUFV3, NDUFS6, NDUFS4, NADH:ubiquinone oxidoreductase core subunit S (NDUFS)8, NDUFS3, NDUFS2 and NDUFS1 were enriched in the Parkinson's disease, Alzheimer's disease and Huntington's disease pathways. Therefore, according to the KEGG database these genes were associated with neurodegenerative disorders. Genes enriched in Module 0, included the mitochondrial ribosomal protein-large (MRPL) and MRP-small (MRPS) gene families. Module 3 included ribosomal protein (RP)S3A, RPL8, RPL15, RPL37A, RPL23A, RPS27L, RPL36AL and signal recognition particle 14. Module 5 consisted of RPS26, RPL3L, signal sequence receptor subunit (SSR)4, Sec6a translocon γ subunit, RPS27A, SSR3 and RPS24. Genes in Modules 0, 3 and 5 were dominantly enriched in pathways associated with ribosomes and/or protein translation (Table I).

Figure 3.

Gene function interaction network of the identified differentially expressed genes between TOF and control groups. TOF, Tetralogy of Fallot.

Table I.

Pathway enrichment analysis of genes in the gene functional modules.

A, Module 0
Pathway enrichment	Protein from module	False discovery rate
Mitochondrial translation (R)	34	<5.000×10−4
Ribosome (K)	20	<5.000×10−4
B, Module 1
Pathway enrichment	Protein from module	False discovery rate
The citric acid cycle and respiratory electron transport (R)	16	<1.667×10−4
Parkinson's disease (K)	16	<1.667×10−4
Alzheimer's disease (K)	16	<1.667×10−4
Huntington's disease (K)	16	<1.667×10−4
Non-alcoholic fatty liver disease (K)	16	<1.667×10−4
C, Module 2
Pathway enrichment	Protein from module	False discovery rate
Hedgehog ligand biogenesis (R)	8	<3.333×10−4
Regulation of apoptosis (R)	8	<3.333×10−4
Proteasome (K)	8	<3.333×10−4
Degradation of beta-catenin by the destruction complex (R)	8	<2.500×10−4
Metabolism of amino acids and derivatives (R)	9	<1.667×10−4
D, Module 3
Pathway enrichment	Protein from module	False discovery rate
Ribosome (K)	7	<1.000×10−3
SRP-dependent co-translational protein targeting to membrane (R)	6	<5.000×10−4
Eukaryotic translation termination (R)	5	<3.333×10−4
Eukaryotic translation elongation (R)	5	<2.500×10−4
Nonsense-mediated decay (R)	5	<2.000×10−4
E, Module 5
Pathway enrichment	Protein from module	False discovery rate
SRP-dependent co-translational protein targeting to membrane (R)	7	<1.000×10−3
Eukaryotic translation termination (R)	4	<5.000×10−4
Eukaryotic translation elongation (R)	4	<3.333×10−4
Nonsense-mediated decay (R)	4	<2.500×10−4
Eukaryotic translation initiation (R)	4	<2.000×10−4
F, Module 6
Pathway enrichment	Protein from module	False discovery rate
Processing of capped intron-containing Pre-mRNA (R)	5	<1.000×10−3
RNA polymerase II transcription (R)	3	4.000×10−3
Processing of capped intron-less pre-mRNA (R)	2	8.330×10−3
Transcriptional regulation of pluripotent stem cells (R)	2	1.230×10−2
Regulatory RNA pathways (R)	2	4.200×10−2

K indicates that the corresponding pathways were searched from the Kyoto Encyclopedia of Genes and Genomes database; R indicates that the corresponding pathways were searched from the Reactome database; SRP, signal recognition particle.

Using the iRegulon plugin, eight transcription factors, including Xbox binding protein 1 (XBP1), ETS transcription factor, forkhead box O1, interferon regulatory factor 8, TATA-box binding protein associated factor 1, zinc finger protein 143, nuclear transcription factor Y subunit α (NFYA) and YY1 transcription factor (YY1) were predicted to regulate DEGs in the FI network (Fig. 4). Among them, XBP1 was an upregulated gene, and its target genes included those encoding ribosomal proteins (including RPS24, RPL37A, RPS27A, RPS27L, RPL23A and RPL36AL), mitochondrial ribosomal proteins (MRPS35, MRPL14, MRPL23 and MRPS18B), as well as NDUFB8, NDUFB5, NDUFS4, NDUFS1, proteasome 26S subunit and non-ATPase 13. In addition, a further two transcription factors (NFYA and YY1) demonstrated differential expression patterns (Log2FC >1) in the TOF samples when compared with the normal controls, although their Log2FC values were lower than the threshold (Table II).

Figure 4.

Gene-transcription factor regulation networks. White triangles represent transcription factors and colored circles represent target genes.

Table II.

Transcription factors predicted to regulate genes in the gene functional interaction networks.

Transcription factor	Log2FC	P-value	Adjusted P-value
NFYA	−1.70543	0.017312	0.066934
XBP1	2.513494	0.000124	0.001634
TAF1	0.065708	0.873505	0.922746
FOXO1A	0.315049	0.449193	0.607251
ZNF143	0.675544	0.103911	0.227741
IRF8	−0.48724	0.423961	0.585214
YY1	1.462285	0.001838	0.012075
ELK4	0.281777	0.31807	0.482523

FC, fold-change; NFYA, nuclear transcription factor Y subunit α; XBP1, Xbox binding protein 1; TAF1, TATA-box binding protein associated factor 1; FOXO1A, forkhead box O1A; ZNF143, zinc finger protein 143; IRF8, interferon regulatory factor 8; YY1, YY1 transcription factor; ELK4, ETS transcription factor.

Discussion

Previously, Bittel et al (11) identified 1,062 DEGs using threshold log2FC and FDR values of ≥2 and <0.05, respectively. In the present study, 878 DEGs were identified using an FDR threshold value of <0.01 and log2FC threshold values >2. Gu et al (31) identified 499 DEGs using the same microarray dataset (GSE26125) with the narrower log2FC threshold values of >2.5 and P-values of <0.01. Therefore, the number of DEGs identified depended on the threshold values used. In addition, the methods and software used for the bioinformatics analyses among these two previous studies and the present study varied. Bittel et al (11) used the Ingenuity Pathway Analysis tool for ontological assessments, which is a curated database and analytical bioinformatics system for identifying interactions, functions and interconnections (networks) between biological molecules. In addition, the differential expression patterns of several genes involved in the WNT or Notch signaling pathways were validated using reverse transcription-quantitative polymerase chain reaction analysis. In the present study, as well as the functional annotation of individual genes, an FI network was constructed and pathway enrichment analysis for each function module was performed. This was used to identify transcription factors predicted to regulate genes in the FI network. Gu et al (31) searched for potential small-molecule drugs by mapping the identified DEGs to the Connectivity Map database.

Bittel et al (11) demonstrated that the majority of the DEGs with abnormal expression were involved in compensatory functions, including hypertrophy, cardiac fibrosis and cardiac dilation, while the WNT and Notch signaling pathways, which are involved in spatial and temporal cell differentiation, appeared to be suppressed (11). In the present study, an FI network based on the identified DEGs was constructed and the biological functions of these genes were investigated using different databases. The results demonstrated that these genes were dominantly enriched in signaling pathways associated with ribosomes and protein translation, as well as neurodegenerative disorders.

The MRPL and MRPS family genes in Module 0 and S-ribosomal protein (RPS) and L-ribosomal protein (RPL) family genes in Modules 3 and 5 were dominantly enriched in signaling pathways associated with ribosomes and protein translation. MRPL family genes, including MRPL40, MRPL10, MRPL13, MRPL12, MRPL14, MRPL17, MRPL16, MRPL54, MRPL18, MRPL38, MRPL39, MRPL33 and MRPL35, encode RPs that form the large subunit of mitochondrial ribosomes. RPS and MRPS family genes, including MRPS35, MRPS16, MRPS33, MRPS15, MRPS14, MRPS27, MRPS23, MRPS7 and MRPS18B, encode RPs that form the small subunit of mitochondrial ribosomes. Meanwhile, RPS26, RPS27A, RPS3A and RPS24 encode RPs of the small ribosomal subunit, and RPL8, RPL15, RPL37A, RPL23A and RPL3L, encode RPs of the large subunit. This suggests that the abnormal expression of cytoplasmic and mitochondrial ribosomes may contribute to the pathogenesis of TOF, which is consistent with the results of a previous study that reported an increase in ribosome number in patients with congenital heart disease (32). In addition, O'Brien et al (33) observed a marked increase in the expression of small nucleolar RNAs (snoRNAs) in the right ventricular myocardium of 16 infants with nonsyndromic TOF, and demonstrated that the target nucleotides of the differentially expressed snoRNAs were primarily 28S and 18S ribosomal RNAs. These results, together those of the present study suggests that the differential expression of genes encoding ribosome subunits may be associated with the dysregulation of snoRNAs. In a previous study that investigated mutations in RPL5 and RPL11 genes in Czech patients with Diamond-Blackfan anemia (34), a mutation in RPL5 in a patient with TOF was reported. Therefore, this study may support the involvement of RPL genes in the pathogenesis of TOF.

In the present study, multiple upregulated genes encoding the subunits of the NADH dehydrogenase complex, including NDUFB8, NDUFB5, NDUFS4 and NDUFS1, were enriched in three signaling pathways associated with Parkinson's disease, Alzheimer's disease and Huntington's disease. This indicates that TOF may share common mechanisms with neurodegenerative disorders. However, the association between NDUFB and NDUFS family genes and TOF has seldom been reported, except for NDUFB5, which has been confirmed to be expressed in mouse heart tissues (35). Despite the lack of direct evidence, previous studies have supported a connection between TOF and neurodegenerative disorders. For instance, Brown et al (36) reported a case of a male infant with infantile neurodegeneration and TOF, and Jinnou et al (37) reported a case of a male infant with pontocerebellar hypoplasia and TOF. These two cases suggest that the pathological mechanisms underlying TOF and neurodegenerative disorders may share common features.

Among the eight transcription factors predicted by the iRegulon plugin in the present study, XBP1 expression was observed to be upregulated in patients in the TOF group. Notably, this gene was predicted to regulate genes encoding the subunits of cytoplasmic and mitochondrial ribosomes, as well as genes involved in neurodegenerative disorders. The XBP1 protein is characterized by its ability to bind the conserved transcriptional Xbox element, which is present in the promoter of the human leukocyte antigen DRα (38). XBP1 is known to be a marker of endoplasmic reticulum stress, a phenomenon that manifests with the accumulation of unfolded proteins in the endoplasmic reticulum (39), which frequently occurs during ischemia/reperfusion following myocardial ischemia (40). However, XBP1 is not a known regulator of the aforementioned genes; therefore, further studies that explore the role of XBP1 in the pathogenesis of TOF in more detail are required.

In conclusion, the results of the present study suggest that the dysregulation of genes encoding the mitochondrial and cytoplasmic ribosomal subunits may contribute to the pathogenesis of TOF via signaling pathways associated with ribosomes and protein translation. In addition, genes encoding the NADH dehydrogenase complex may contribute to the pathogenesis of this disease via neurodegenerative disorder-associated signaling pathways. As the transcription factor XBP1 was predicted to be implicated in the regulation of genes involved in these signaling pathways, it may therefore be involved in the pathogenesis of TOF. However, this is yet to be validated in future studies. The results of the present study provide an in-depth insight into the molecular mechanisms underlying the pathogenesis of TOF.