Abstract

Atrial septal defects (ASDs) are a common human congenital heart disease (CHD) that can be induced by genetic abnormalities. Our previous studies have demonstrated a genetic interaction between Tbx5 and Osr1 in the second heart field (SHF) for atrial septation. We hypothesized that Osr1 and Tbx5 share a common signaling networking and downstream targets for atrial septation. To identify this molecular networks, we acquired the RNA-Seq transcriptome data from the posterior SHF of wild-type, Tbx5+/, Osr1+/−, Osr1−/− and Tbx5+/−/Osr1+/− mutant embryos. Gene set analysis was used to identify the Kyoto Encyclopedia of Genes and Genomes pathways that were affected by the doses of Tbx5 and Osr1. A gene network module involving Tbx5 and Osr1 was identified using a non-parametric distance metric, distance correlation. A subset of 10 core genes and gene–gene interactions in the network module were validated by gene expression alterations in posterior second heart field (pSHF) of Tbx5 and Osr1 transgenic mouse embryos, a time-course gene expression change during P19CL6 cell differentiation. Pcsk6 was one of the network module genes that were linked to Tbx5. We validated the direct regulation of Tbx5 on Pcsk6 using immunohistochemical staining of pSHF, ChIP-quantitative polymerase chain reaction and luciferase reporter assay. Importantly, we identified Pcsk6 as a novel gene associated with ASD via a human genotyping study of an ASD family. In summary, our study implicated a gene network involving Tbx5, Osr1 and Pcsk6 interaction in SHF for atrial septation, providing a molecular framework for understanding the role of Tbx5 in CHD ontogeny.

Introduction

Atrial septation is the division of a common atrium into the left and right atria that are observed in a structurally mature, four-chambered tetrapod heart. The atrial septum is formed by three distinct structures that undergo concurrent morphogenesis: the atrioventricular endocardial cushions, the primary atrial septum (PAS) and the spina vestibuli or dorsal mesenchymal protrusion (SV/DMP) (1–3). PAS and the SV/DMP eventually contact and fuse with the atrioventricular endocardial cushions after the atrioventricular canal has expanded rightward. With this fusion, atrial septation completes and ensures that blood flows only through the right atrioventricular canal into the right ventricle.

The failure of atrial septation results in an atrial septal defect (ASD), which is one of the most common forms of human congenital heart disease (CHD). This disease accounts for up to 40% of clinically relevant acyanotic shunts in adults (4). Missing one or more of the three atrial septal components, ASDs are classified by distinct structural defects. Common ASDs are characterized by a common atrium that lacks the PAS and the SV/DMP. Secundum ASDs result from the absence or hypoplasia of the PAS. Primum ASDs, a form of atrioventricular septal defect (AVSD), occur in the DMP region adjacent to the mitral and tricuspid valves. Historically, AVSDs have been attributed to developmental abnormalities of the atrioventricular endocardial cushions or their derivatives (3). However, recent work has suggested that deficiencies of the DMP could explain some forms of AVSDs (5–7).

Using mouse models, several studies have demonstrated sources of progenitor cells for atrial septal developemnt. A role for the second heart field (SHF) cardiac progenitor cells in atrial septation has been identified (8–13). The developing DMP is marked by SHF markers Mef2c-AHF:Cre (13–16) and Isl1 (16), and is generated by the cell lineage receiving Hedgehog (Hh) signaling (12,17,18). Progenitors of much if not all of the atrial septum, including the PAS and SV/DMP, are SHF derivatives from E8 to E10 (6,18).

Studies in vertebrate model organisms have identified important transcription factors involved in atrial septation, including Gata4, Nkx2–5 and Tbx5, in which mutations were identified and shown to cause ASDs in humans (19–22). These transcription factors are thought to control a conserved program that triggers activation and expression of specific signaling molecules in regulating the genesis of the cardiac septum from mesodermal stem cells. Tbx5 is a member of T-box transcription factor family and acts as a key regulator during early cardiac morphogenesis. Haploinsufficiency of Tbx5 causes Holt–Oram syndrome (HOS) in humans. HOS is an autosomal-dominant disease affecting one of every 100 000 live births. HOS is also characterized by forelimb deformities combined frequently with congenital heart defects. It is known that Tbx5 binds to the T-box binding elements and activates the transcription of its downstream target genes by acting as a co-activator. Furthermore, Tbx5 interacts with Nkx2–5 and Gata4, and positively regulates transcription within the developing heart. Previous work has shown that Tbx5 is required in the posterior SHF (pSHF) for atrial septation (23). Our recent study has elucidated an interaction between Tbx5 and Osr1, required in E9.5 pSHF for cardiac progenitor proliferation and for atrial septation (24). Osr1 is a known transcription factor whose mutation causes common atrium in mouse embryos (24,25). However, the genes that are closely involved in the interaction network of Tbx5 and Osr1 in atrial septal progenitors during atrial septation remain almost unknown.

In this study, we investigated the gene regulatory module involving Tbx5 and Osr1 in the pSHF development using a systems biology approach. The global transcriptional fluctuation induced by the various dose levels of Tbx5 and Osr1 in the pSHF allowed us to construct a gene network model involving Tbx5 and Osr1 interaction.

Results

Tbx5+/−, Tbx5+/−/Osr1+/− and Osr1−/− mice show global gene-expression variations in the SHF at E9.5

Previous work has shown that Tbx5 directly regulates the transactivation of Osr1 in the pSHF for atrial septation (8,24). We have evaluated the incidence of ASDs in mutant embryos of Tbx5 and Osr1 in order to understand the interactions between Tbx5 and Osr1 (Fig. 1). The Tbx5 haploinsufficiency had 40% of the embryos developing AVSDs (Fig. 1C and H) as reported earlier (8,24,25). The Osr1tm1RJ/+ mice, referred to as Osr1+/− mice, have a lacZ gene inserted into the first coding exon of the Osr1 gene, which results in the loss of Osr1-gene function (25). Although only 50% of Osr1−/− embryos survive through E13.5 (18/37), 75% of them displayed AVSDs (Fig. 1E and J), whereas no AVSDs were seen in either Osr1+/− (Fig. 1B and G) or wild-type embryos (Fig. 1A and F). There was significantly higher incidence of AVSDs of Tbx5+/−/Osr1+/− mouse embryos (11/11, 100%, Fig. 1D and I) than in either Osr1+/− embryos (0/18, P = 0.000) or Tbx5+/− embryos (5/15, P = 0.0001).

Figure 1.

Tbx5 and Osr1 interact in the pSHF. (AJ) Histology of Tbx5 and Osr1 transgenic mouse embryo hearts at E13.5 An asterisk (*) indicates missing structures of atrial septation. Magnification in (A–E) is ×40; magnification in (F–J) is ×100. PAS, primary atrial septum; DMP, dorsal mesenchymal protrusion. (K) Expression levels of Gapdh, Nkx2-5, Isl1, Tbx1, Tbx5, Osr1 and Myl2 in the micro-dissected heart. The pSHF and aSHF of E9.5 wild-type mouse embryos were evaluated by RNA-Seq analysis. Different letters indicate differential gene expression (P < 0.05) between the tissues (n = 6). (L) Expression levels of Tbx5 and Osr1 in the micro-dissected pSHF of E9.5 mouse embryos were evaluated by real-time PCR analysis. Expression levels of Tbx5+/−, Tbx5+/−/Osr1+/−, Osr1+/− or Osr1−/− versus wild-type embryos were calculated by the ΔΔCt method (26). Different letters indicate differential gene expression (P < 0.05) between the genotypes (n = 3).

Figure 1.

Tbx5 and Osr1 interact in the pSHF. (AJ) Histology of Tbx5 and Osr1 transgenic mouse embryo hearts at E13.5 An asterisk (*) indicates missing structures of atrial septation. Magnification in (A–E) is ×40; magnification in (F–J) is ×100. PAS, primary atrial septum; DMP, dorsal mesenchymal protrusion. (K) Expression levels of Gapdh, Nkx2-5, Isl1, Tbx1, Tbx5, Osr1 and Myl2 in the micro-dissected heart. The pSHF and aSHF of E9.5 wild-type mouse embryos were evaluated by RNA-Seq analysis. Different letters indicate differential gene expression (P < 0.05) between the tissues (n = 6). (L) Expression levels of Tbx5 and Osr1 in the micro-dissected pSHF of E9.5 mouse embryos were evaluated by real-time PCR analysis. Expression levels of Tbx5+/−, Tbx5+/−/Osr1+/−, Osr1+/− or Osr1−/− versus wild-type embryos were calculated by the ΔΔCt method (26). Different letters indicate differential gene expression (P < 0.05) between the genotypes (n = 3).

Since Tbx5 and Osr1 interact in pSHF at E9.5 for atrial septation (24), we collected pSHF at E9.5 by microdissection. To validate that our microdissection correctly isolates the pSHF, aSHF and heart, we analyzed gene expression data of these three isolated tissues at E9.5 from six isolated biological replicates via RNA-Seq analysis (Fig. 1K). Nkx2-5 and Myl2 mark the heart at E9.5, while Isl1 and Tbx1 are widely accepted as markers for SHF (27,28). It is also reported that Tbx5 and Osr1 are strongly expressed in pSHF, but rarely expressed in the aSHF (8,25), whereas Tbx1 is more highly expressed in aSHF (27–29). Our RNA-Seq analysis results confirmed this tissue-specific expression pattern (Fig. 1K).

We tested the hypothesis that the different incidence of ASDs in our mutant mouse embryos was associated with the Tbx5- and Osr1-expression variations via evaluating the expression levels of Tbx5 and Osr1 in pSHF at E9.5. The Tbx5 expression level had no changes between Tbx5+/− and Tbx5+/−/Osr1+/− or among Osr1+/−, Osr1+/− and wild-type (Fig. 1L). Osr1 showed a dose gradient across wild-type, Osr1+/−, Tbx5+/−, Tbx5+/−/Osr1+/− and Osr1−/−. As verified using real-time polymerase chain reaction (PCR) analysis, when the wild-type dose level was set to 1 (100%), the expression levels of Osr1 for Osr1+/−, Tbx5+/−, Tbx5+/−/Osr1+/− and Osr1−/− mice were 0.73 ± 0.15, 0.54 ± 0.17, 0.47 ± 0.19 and 0.25 ± 0.08, respectively (Fig. 1L). Therefore, we obtained five different levels, although some of them being close, of Osr1 expression gradients using different transgenic mouse lines.

In order to identify important genes involved in atrial septation, we performed RNA-Seq experiments for E9.5 pSHF of wild-type, Osr1+/−, Tbx5+/−, Tbx5+/−/Osr1+/− and Osr1−/− mouse embryos. For each sample, 20–25 million short reads were generated using Illumina HiSeq 2000, and >90% of short reads were mapped to the mouse reference genome GRCm38. Principal components analysis (PCA) showed that the five genotypes have distinct expression patterns (Fig. 2A). The first principal component (PC1) accounted for 57.95% of variation (standard deviation) in all transcriptome data, and the second principal component (PC2) accounted for 16.68% of variation. The heat map for the top 6989 variable genes (standard deviation > 1) also showed that dose changes in Tbx5 and Osr1 affected the global gene expression (Fig. 2B). The Tbx5 heterozygotes and the Osr1 heterozygotes were clustered in a sub-group, which further clustered together with Tbx5 and Osr1 compound heterozygotes. This gene heterozygotes group was more closely related to wild-type than Osr1−/− samples.

Figure 2.

Global gene-expression variations in the SHF of Tbx5+/−, Tbx5+/−/Osr1+/−, Osr1+/− and Osr1−/− mice embryos at E9.5. (A) The six samples were plotted on the first two principal components, PC1 and PC2. (B) The heat map of global gene expression levels for the 6989 genes that have the largest variation across all samples (standard deviation > 1).

Figure 2.

Global gene-expression variations in the SHF of Tbx5+/−, Tbx5+/−/Osr1+/−, Osr1+/− and Osr1−/− mice embryos at E9.5. (A) The six samples were plotted on the first two principal components, PC1 and PC2. (B) The heat map of global gene expression levels for the 6989 genes that have the largest variation across all samples (standard deviation > 1).

Osr1 dosage affects gene-signaling pathways and gene regulatory networks in pSHF

The differentially expressed genes (DEGs) related with Osr1 dosage were identified using regression analysis across all RNA-Seq samples against Osr1 expression levels. Regression analysis and identification of DEGs was also performed using R programing language. We found 431 genes up-regulated and 407 genes down-regulated by Osr1 using a P-value cutoff of 0.05. We performed real-time PCR analysis for 55 selected transcripts. Thirty-nine genes were confirmed to have significant expression changes (70%). Among those genes, Osr1 and Cdk6 have been previously suggested as direct transcription targets of Tbx5 in the pSHF (23).

In order to identify the pathways altered by the Osr1 dosage, we conducted pathway analysis using a hypergeometric distribution-based gene set enrichment analysis method (30). The pathways were compiled from Kyoto Encyclopedia of Genes and Genomes (KEGGs, www.genome.jp/kegg). A set of 16 pathways including cell cycle, adherents junction, Wnt-signaling pathway and Hedgehog (Hh)-signaling pathway were found to be positively associated with the expression level of Osr1 (Table 1A), and 6 pathways were found to be negatively associated (Table 1B).

Table 1.

Gene-signaling pathways associated with the level of Osr1 in pSHF

Pathway P-value 
A. Pathways whose integrity are positively associated with Osr1 expression 
Melanogenesis 0.00741 
Cell cycle 0.003475 
Endocytosis 0.004134 
Wnt-signaling pathway 0.001436 
Hh-signaling pathway 0.008874 
Pathways in cancer 0.002459 
Basal-cell carcinoma 0.001791 
Adherens junction 0.000148 
Thyroid cancer 0.000744 
Fc gamma R-mediated phagocytosis 0.005358 
GnRH-signaling pathway 0.00741 
Vasopressin-regulated water reabsorption 0.001092 
Asthma 0.0067 
Oocyte meiosis 0.004699 
Antigen processing and presentation 0.008874 
Bile secretion 0.000111 
B. Pathways whose integrity are negatively associated with Osr1 expression 
Cytokine–cytokine receptor interaction 0.005898 
Purine metabolism 0.00765 
Pyrimidine metabolism 5.50E−05 
RNA polymerase 6.20E−05 
Basal transcription factors 0.000438 
Ubiquinone and other terpenoid-quinone biosynthesis 1.02E−08 
Pathway P-value 
A. Pathways whose integrity are positively associated with Osr1 expression 
Melanogenesis 0.00741 
Cell cycle 0.003475 
Endocytosis 0.004134 
Wnt-signaling pathway 0.001436 
Hh-signaling pathway 0.008874 
Pathways in cancer 0.002459 
Basal-cell carcinoma 0.001791 
Adherens junction 0.000148 
Thyroid cancer 0.000744 
Fc gamma R-mediated phagocytosis 0.005358 
GnRH-signaling pathway 0.00741 
Vasopressin-regulated water reabsorption 0.001092 
Asthma 0.0067 
Oocyte meiosis 0.004699 
Antigen processing and presentation 0.008874 
Bile secretion 0.000111 
B. Pathways whose integrity are negatively associated with Osr1 expression 
Cytokine–cytokine receptor interaction 0.005898 
Purine metabolism 0.00765 
Pyrimidine metabolism 5.50E−05 
RNA polymerase 6.20E−05 
Basal transcription factors 0.000438 
Ubiquinone and other terpenoid-quinone biosynthesis 1.02E−08 

Our previous studies have shown that integrity of Hh-signaling pathway is dependent on the level of Tbx5 and Osr1 (8,24). Here, we performed real-time PCR analysis to quantify expression of several key genes of the cell-cycle pathway and Wnt-signaling pathway, using the pSHF of wild-type, Osr1+/− and Tbx5+/−/Osr1+/− mouse embryos at E9.5. Tested genes were selected by combining the gene-expression profile of our RNA-sequencing data of E9.5 pSHF Osr1 mutant embryos, the microarray result of E9.5 pSHF of Tbx5+/− mouse embryos (8) and the ChIP-Seq data of Tbx5 in HL-1 cells (31). Results of this experiment could enable us to verify our pathway analysis. Cdk6 and Ccnd2, critical regulators for G1-S phase transition, were significantly reduced in all three transgenic mouse embryos. The Ccnd3 level was reduced in pSHF of both Osr1+/−and Tbx5+/−/Osr1+/− embryos, and the Ccnd1 level was reduced in pSHF of Osr1+/− embryos only (Fig. 3A). Interestingly, Cdk6 expression was further down-regulated in Tbx5+/−/Osr1+/− mouse embryos, compared with littermate controls, Osr1+/− or Tbx5+/− mouse embryos (Fig. 3A). For those examined genes involved in or regulated by the Wnt-signaling pathway, we observed significantly decreased expression of Wnt2, Wnt3a, Dkk1, Vegf-α (32–34) and Fn1 (35,36) in Osr1+/− embryos. Wnt2 expression was significantly reduced in all three mutant lines (Fig. 3B).

Figure 3.

Expression levels of key genes involved in the cell-cycle pathway and Wnt-signaling pathway. Expression of cell-cycle genes (A) and Wnt pathway genes (B) in the micro-dissected pSHF of E9.5 Tbx5+/−, Tbx5+/−/Osr1+/− or Osr1+/− versus wild-type embryos were evaluated by real-time PCR analysis and calculated by the ΔΔCt method (26). Asterisk (*) indicates significant difference (P < 0.05) (n = 3–5).

Figure 3.

Expression levels of key genes involved in the cell-cycle pathway and Wnt-signaling pathway. Expression of cell-cycle genes (A) and Wnt pathway genes (B) in the micro-dissected pSHF of E9.5 Tbx5+/−, Tbx5+/−/Osr1+/− or Osr1+/− versus wild-type embryos were evaluated by real-time PCR analysis and calculated by the ΔΔCt method (26). Asterisk (*) indicates significant difference (P < 0.05) (n = 3–5).

A Tbx5–Osr1 network module is predicted in pSHF and validated by time-course microarray for P19CL6 cells

Using the pSHF RNA-Seq data of wild-type, Osr1+/−, Tbx5+/−, Osr1−/− and Tbx5+/−/Osr1+/− mouse embryos at E9.5, we constructed a co-expression network by applying the distance correlation (DC) matrix between all gene pairs. DC takes into account not only linear relationships, but also non-linear relationships between genes (37,38). The DC matrix was then transformed into an adjacency matrix and hierarchical clustering was used to identify the network hub genes. The network hub containing Osr1 was compared with the gene set of 255 CHD related genes curated by Lage's research group (39), and the resulting gene networks module includes 10 core genes, Tbx5, Osr1, Tll1, Pcsk6, Dtna, Angpt1, Hand2, Mapk7, Trex1 and Fgf8 (Fig. 4A). When the DCs were calculated between the network module genes, Pcsk6 showed a short distance to Tbx5 and Angpt1, with high DCs equal to 0.708 and 0.714, respectively (Fig. 4A).

Figure 4.

A Tbx5–Osr1 network module is involved in the pSHF development. (A) 10 genes were found to form a network module that had intensive interactions between them. Each network edge (black line) represents the short distance (high-DC value) between the two genes. (B) The heat map of the expression levels of the network module genes for time-course P19CL6 microarray experiments. Red represents high-gene expression and green represents low expression.

Figure 4.

A Tbx5–Osr1 network module is involved in the pSHF development. (A) 10 genes were found to form a network module that had intensive interactions between them. Each network edge (black line) represents the short distance (high-DC value) between the two genes. (B) The heat map of the expression levels of the network module genes for time-course P19CL6 microarray experiments. Red represents high-gene expression and green represents low expression.

In order to validate the network module genes, we re-analyzed the previous published time-course microarray data of P19CL6 cells (40). P19CL6 cell constitutes a strain of pluripotent embryonic stem cells that can differentiate into spontaneously beating cardiac myocytes in the presence of dimethyl sulfoxide (DMSO). Total RNA samples were collected from P19CL6 cells growing in the presence of DMSO for 0, 4, 6, 8, 10 or 12 days and microarray data were collected using Affymetrix Mouse Genome 430 2.0 array. We performed analysis of variance (ANOVA) and identified 4,816 genes that were differentially expressed across the time points by controlling false discovery rate at 10%. Most of the network module genes have been found to be differentially expressed in the time-course, including Dtna, Mapk7, Tbx5, Angpt1, and Hand2. As shown in the heatmap in Figure 4B, the expression levels of Angpt1, Tbx5 and Hand2 increased in the later stages of the differentiation, whereas Dtna and Mapk7 had higher expression levels in early stages. Trex1 data are missing because no probe is designed for the Trex1 gene in the Affymetrix 430 2.0 array.

The Tbx5–Osr1 network module is validated by gene-expression variations in pSHF of Tbx5 and Osr1 mutant mouse embryos at E9.5

As this network module is based on the gene-expression variations between the mouse lines, we performed real-time PCR using our pSHF of wild-type, Tbx5+/−, Osr1−/− and Tbx5+/−/Osr1+/− mouse embryos at E9.5 to validate our network module at the molecular level. According to our network model, the genes that had the shortest distance to Tbx5/Osr1 would have the highest possibility of being their direct downstream targets. We showed that Osr1, Pcsk6 and Tll1, which directly link to Tbx5, were expressed significantly lower within the pSHF, but not within the heart of Tbx5+/ embryos compared with wild-type embryos (Fig. 5A). The other module genes including Dtna, Angpt, Mapk7, Hand2, Trex1 and Fgf8, whose expression levels were predicted to be dependent on the expression of Osr1 and Tbx5, but were not their direct targets, were located relatively farther from Osr1 and Tbx5 in this network module. When comparing the transcriptional levels across wild-type, Tbx5+/−, Osr1+/−, Tbx5+/−/Osr1+/− for the six genes that did not directly connect to Tbx5, two genes, Dtna and Hand2, were found to have significant transcriptional variations across the four genotypes (P < 0.05), whereas Fgf8 and Mapk7 had marginally significant variations (P < 0.1) (Fig. 5B). The remaining two non-significant genes, Trex1 and Angpt1, as well as Dtna and Tbx5 who are directly linked to Osr1 in the module, were further investigated in Osr1−/− mice. Trex1 and Dtna showed significantly lower transcriptional levels in Osr1−/− compared with wild-type (Fig. 5C).

Figure 5.

A Tbx5–Osr1 network module was verified in mouse tissue and P19CL6 cells. (A) Expression levels of Tll1, Pcsk6 and Osr1 in the micro-dissected pSHF and heart of E9.5 mouse embryos were evaluated using real-time PCR analysis. Expression of the detected genes of Tbx5+/− versus wild-type embryos was calculated by the ΔΔCt method (26). Double asterisks indicate a significant difference (P < 0.05, n = 3). (B) Expression levels of Dtna, Hand2, Trex1, Mapk7, Fgf8 and Angpt1 in the micro-dissected pSHF of E9.5 mouse embryos were evaluated using real-time PCR analysis. The ΔCt value is calculated versus Gapdh. Double asterisks indicate a significant difference among the embryos of the four genotypes (P < 0.05, n = 3). Single asterisk indicates a difference among the four genotypes (P < 0.1, n = 3). (C) Expression levels of Trex1, Angpt1, Dtna and Tbx5 in the micro-dissected pSHF of E9.5 mouse embryos were evaluated by real-time PCR analysis. Expression of the detected genes of Osr1−/− versus wild-type embryos was calculated by the ΔΔCt method (26). Asterisk (*) indicates significant difference (P < 0.05, n = 3).

Figure 5.

A Tbx5–Osr1 network module was verified in mouse tissue and P19CL6 cells. (A) Expression levels of Tll1, Pcsk6 and Osr1 in the micro-dissected pSHF and heart of E9.5 mouse embryos were evaluated using real-time PCR analysis. Expression of the detected genes of Tbx5+/− versus wild-type embryos was calculated by the ΔΔCt method (26). Double asterisks indicate a significant difference (P < 0.05, n = 3). (B) Expression levels of Dtna, Hand2, Trex1, Mapk7, Fgf8 and Angpt1 in the micro-dissected pSHF of E9.5 mouse embryos were evaluated using real-time PCR analysis. The ΔCt value is calculated versus Gapdh. Double asterisks indicate a significant difference among the embryos of the four genotypes (P < 0.05, n = 3). Single asterisk indicates a difference among the four genotypes (P < 0.1, n = 3). (C) Expression levels of Trex1, Angpt1, Dtna and Tbx5 in the micro-dissected pSHF of E9.5 mouse embryos were evaluated by real-time PCR analysis. Expression of the detected genes of Osr1−/− versus wild-type embryos was calculated by the ΔΔCt method (26). Asterisk (*) indicates significant difference (P < 0.05, n = 3).

Pcsk6 is confirmed as a downstream gene of Tbx5 both in vitro and in vivo

We hypothesized that Tbx5 may directly regulate the transcriptional activation of Pcsk6 or Tll1. Therefore, we interrogated their genomic loci to identify potential Tbx5-responsive elements. Though we did not identify any potential Tbx5-responsive genomic region for Tll1, we found three conserved Tbx5-binding sites within the intron regions of Pcsk6 by mapping the ChIP-Seq peaks in the HL-1 cardiomyocyte cell line (31) onto the evolutionary conservative regions in mouse genome. We, therefore, subcloned two Pcsk6 genomic regions containing these binding sites into pGL3 basic luciferase reporter vectors: Pcsk6-Fr1 (chr7: 73047356-73047737) and Pcsk6-Fr2 (chr7: 73182155-73182432) (Table 2 and Fig. 6A). The luciferase reporter assay showed that Tbx5 significantly increased firefly luciferase expression for Pcsk6-Fr2 (P = 0.0005 versus pGL3), but not for Pcsk6-Fr1. We subsequently investigate whether the mutations in the Tbx5 binding motifs of Pcsk6-Fr2 would repress the transactivation of luciferase expression. For the two Tbx5 binding motifs in the Pcsk6-Fr2 region, we introduced a mutation for each of them: Pcsk6-Fr2m1 and Pcsk6-Fr2m2. Both Pcsk6-Fr2m1 and Pcsk6-Fr2m2 lost Tbx5-induced transactivation (Fig. 6B: For Pcsk6-Fr2m1 P = 0.2702 versus pGL3; for Pcsk6-Fr2m2 P = 0.0553 versus pGL3).

Table 2.

Genomic regions assessed in luciferase reporter assay and ChIP assay

 Luciferase assay
 
 Cloned region Tbx5 binding site Deleted binding site 
Pcsk6-Fr1 chr7:73047356-73047737 chr7:73047700-73047708  
Pcsk6-Fr2 chr7:73182155-73182432 chr7:73182205-73182213
chr7:73182219-73182227
chr7:73182226-73182234 
 
Pcsk6-Fr2m1 chr7:73182219-73182234 chr7:73182205-73182213 
Pcsk6-Fr2m2 chr7:73182205-73182213 chr7:73182219-73182234 
 ChIP 
 Amplified region Tbx5 binding site  
Pcsk6-Fr2 chr7:73182155-3182321 chr7:73182219-73182227
chr7:73182226-73182234 
 
Neg ctl chr7:73104575-73104723   
 Luciferase assay
 
 Cloned region Tbx5 binding site Deleted binding site 
Pcsk6-Fr1 chr7:73047356-73047737 chr7:73047700-73047708  
Pcsk6-Fr2 chr7:73182155-73182432 chr7:73182205-73182213
chr7:73182219-73182227
chr7:73182226-73182234 
 
Pcsk6-Fr2m1 chr7:73182219-73182234 chr7:73182205-73182213 
Pcsk6-Fr2m2 chr7:73182205-73182213 chr7:73182219-73182234 
 ChIP 
 Amplified region Tbx5 binding site  
Pcsk6-Fr2 chr7:73182155-3182321 chr7:73182219-73182227
chr7:73182226-73182234 
 
Neg ctl chr7:73104575-73104723   

All genomic coordinates are shown in mouse genome build mm9.

Figure 6.

Pcsk6 was identified as a Tbx5 downstream target in pSHF for atrial septation. (A) Schematic of mouse Pcsk6 genomic locus including Tbx5 binding regions (31) and the cloned genomic fragments Pcsk6-fr1 and Pcsk6-fr2 used for Tbx5 regulation assays: luciferase reporter assay and ChIP-qPCR. (B) Tbx5 stimulated firefly luciferase expression in wild-type Pcsk6-Fr2, but not in wild-type Pcsk6-Fr1 and mutated Pcsk6-Fr2 fragments Pcsk6-Fr2m1 and Pcsk6-Fr2m2). Results are presented as mean ± SEM (*P < 0.05 when compared with pGL3, #P < 0.05 when compared with Pcsk6-Fr2, n = 3). (C) Enrichment of Tbx5 responsive Pcsk6 genomic fragment in the SHF by Tbx5 ChIP-qPCR. Results are presented as mean ± SEM (**P < 0.01, n = 3). (D–G) Expression of Pcsk6 was shown in pSHF of wild-type and Tbx5+/− mouse embryos by IHC staining. Red arrow heads indicate typical staining of Pcsk6 in cytosol. Magnification: (D and F, 100×; E and G, 200×). (H) Basic pedigree structure of an ASD family. ASD subjects are represented by black; a subclinical phenotype, septal aneurysm, is indicated by gray. Asterisks show the subjects that were analyzed using SNP arrays. (I) LOD scores of Pcsk6 SNPs. The loci of the SNPs and the location of Pcsk6 exons are shown on chromosome 15.

Figure 6.

Pcsk6 was identified as a Tbx5 downstream target in pSHF for atrial septation. (A) Schematic of mouse Pcsk6 genomic locus including Tbx5 binding regions (31) and the cloned genomic fragments Pcsk6-fr1 and Pcsk6-fr2 used for Tbx5 regulation assays: luciferase reporter assay and ChIP-qPCR. (B) Tbx5 stimulated firefly luciferase expression in wild-type Pcsk6-Fr2, but not in wild-type Pcsk6-Fr1 and mutated Pcsk6-Fr2 fragments Pcsk6-Fr2m1 and Pcsk6-Fr2m2). Results are presented as mean ± SEM (*P < 0.05 when compared with pGL3, #P < 0.05 when compared with Pcsk6-Fr2, n = 3). (C) Enrichment of Tbx5 responsive Pcsk6 genomic fragment in the SHF by Tbx5 ChIP-qPCR. Results are presented as mean ± SEM (**P < 0.01, n = 3). (D–G) Expression of Pcsk6 was shown in pSHF of wild-type and Tbx5+/− mouse embryos by IHC staining. Red arrow heads indicate typical staining of Pcsk6 in cytosol. Magnification: (D and F, 100×; E and G, 200×). (H) Basic pedigree structure of an ASD family. ASD subjects are represented by black; a subclinical phenotype, septal aneurysm, is indicated by gray. Asterisks show the subjects that were analyzed using SNP arrays. (I) LOD scores of Pcsk6 SNPs. The loci of the SNPs and the location of Pcsk6 exons are shown on chromosome 15.

We next assessed whether Tbx5 physically occupies the Tbx5-responsive region in Pcsk6 locus in the SHF in vivo during atrial septum progenitor specification. The ChIP-qPCR was performed for Tbx5 using micro-dissected 15 SHF tissues of wild-type mouse embryos at E9.5. ChIP-qPCR demonstrated significant Tbx5-dependent enrichment for Pcsk6-Fr2 fragment, but not for the upstream or downstream control fragments subcloned from the neighboring genomic regions (Table 2 and Fig. 6C).

We next performed immunohistochemical staining to detect the expression changes in Pcsk6 in pSHF of Tbx5+/− mouse embryos comparing to that of wild-type at E9.5. Strong Pcsk6 expression was found within the dorsal mesocardium at the inflow region (pSHF) of the wild-type mouse embryo (Fig. 6D and E), however, its expression significantly decreased within the pSHF of the Tbx5+/− mouse embryos (Fig. 6F and G).

Human familial study verified the role of PCSK6 in ASD

Our data have suggested that the network module genes are potentially involved in atrial septation, therefore, we sought to verify the gene associations in a Spanish family with five ASD patients and four members having interatrial septal aneurysm. The genomic DNAs were collected from the whole blood samples of four ASD patients, three septal aneurysm patients and three healthy people (Fig. 6H). The DNAs were assessed by SNP array experiment using the Illumina Human660W-Quad chip. The genotypes of a total of 657 366 SNPs were generated. We assessed a total number of 254 SNPs that were located within 50 Kbp of the 10 network module genes. Three SNPs were found significant for Hardy–Weinberg equilibrium tests, and were eliminated from further analysis. Only one SNP, rs4965381, was found to have a logarithm of odds (LOD) score >2 (Fig. 6I). Table 3 listed the results for all PCSK6 SNPs. SNP rs4965381 is located in the sixth intron of the PCSK6 gene. The P-values were calculated using χ2 tests. SNP rs4965381 has a P-value of 0.004 (Table 3). Thus, we found one SNP for PCSK6 was significantly associated with the ASD phenotype in this Spanish family.

Table 3.

Summary of SNPs in the Pcsk6 gene in the human ASD family study

SNP ID SNP GenomeBuild Chr MapInfo LOD P-value 
rs12437510 [A/G] 37.1 15 1.02E + 08 0.004636 0.023019 
rs12900435 [T/G] 37.1 15 1.02E + 08 0.262085 0.056467 
rs12910524 [A/G] 37.1 15 1.02E + 08 −0.13109 0.084163 
rs12916267 [A/G] 37.1 15 1.02E + 08 0.261162 0.023019 
rs1532364 [T/C] 37.1 15 1.02E + 08 0.261162 0.023019 
rs3784485 [T/C] 37.1 15 1.02E + 08 0.004636 0.023019 
rs4965373 [T/C] 37.1 15 1.02E + 08 −0.13109 0.084163 
rs4965381 [T/C] 37.1 15 1.02E + 08 2.10721 0.004087 
rs6598464 [T/C] 37.1 15 1.02E + 08 −0.13109 0.230693 
rs7169026 [A/G] 37.1 15 1.02E + 08 −0.13109 0.230693 
SNP ID SNP GenomeBuild Chr MapInfo LOD P-value 
rs12437510 [A/G] 37.1 15 1.02E + 08 0.004636 0.023019 
rs12900435 [T/G] 37.1 15 1.02E + 08 0.262085 0.056467 
rs12910524 [A/G] 37.1 15 1.02E + 08 −0.13109 0.084163 
rs12916267 [A/G] 37.1 15 1.02E + 08 0.261162 0.023019 
rs1532364 [T/C] 37.1 15 1.02E + 08 0.261162 0.023019 
rs3784485 [T/C] 37.1 15 1.02E + 08 0.004636 0.023019 
rs4965373 [T/C] 37.1 15 1.02E + 08 −0.13109 0.084163 
rs4965381 [T/C] 37.1 15 1.02E + 08 2.10721 0.004087 
rs6598464 [T/C] 37.1 15 1.02E + 08 −0.13109 0.230693 
rs7169026 [A/G] 37.1 15 1.02E + 08 −0.13109 0.230693 

Discussion

Recent work has helped define the morphogenesis of the atrial septum. However, the molecular pathways underlying molecular ontogeny of atrial septation remain unclear. Long-standing genetic studies of large families have shown that haploinsufficiency of Gata4, Nkx2-5 and Tbx5 causes ASDs in human (19–22). Our previous work has shown that Osr1 is a downstream target of Tbx5 and the two genes co-operatively regulate the proliferation of cardiac progenitors in the pSHF (8,24). However, studies aimed to identify the important genes and gene–gene interactions in atrial septal development have made little progress. One major challenge is that cardiovascular genesis follows not a linear hierarchical model but an integrative, mechanistic model. Human heart development is a dynamic and multiple-step process regulated by a complex gene regulatory network. Therefore, it is difficult to study the mechanism of heart development by the reductionist approach. In our study, a network module involving Tbx5 and Osr1 interaction was identified using a systems biology approach, which provides a more efficient way to decode the gene network for a complex human disease, such as ASD. The transcriptional changes in the module genes and the gene–gene interactions were validated by real-time PCR analysis in pSHF of Tbx5 and Osr1 transgenic mouse embryos, luciferase reporter assay, ChIP-qPCR, time-course gene expression profiles in P19CL6 cells and SNP array data of human ASD patients. Our finding implicates a molecular network involving Tbx5/Osr1/Pcsk6 in the pSHF for atrial septation.

One major challenge of gene networks analysis remains to be the quantification of the gene–gene associations. Gene networks analysis is usually based on gene co-expression such that the genes showing similar expression patterns or trends are linked together (41,42). The standard approach to detect gene co-expression is to use correlation coefficients, such as Pearson and Spearman correlations (42,43). However, correlation coefficients only reflect the linear relationship in expression between two genes, and Pearson correlation tests further require that expression levels follow normal distributions. These assumptions are unlikely held in many gene–gene interactions. In this study, we used an advanced non-parametric statistic called DC that was recently proposed to measure the dependency between two random variables. The DC is able to detect both linear and non-linear relationships between genes, and it is a non-parametric statistic so that it makes no assumption about data distribution.

The 10 interactive genes in the pSHF of mouse embryos at E9.5 were predicted to play a role in atrial septation. In this study, we have confirmed that transcription of those module genes, except Angpt, were dependent on Tbx5 and Osr1 interaction in E9.5 pSHF of transgenic mouse embryos. As summarized in Table 4, eight genes have been reported to be involved in heart development. Among them, mutations of six genes including Tbx5 (8,21,24,44–46), Osr1 (8,25), Tll1 (50,51), Hand2 (54), Dtna (60) and Pcsk6 (52,53) either cause ASDs in the mouse or are associated with higher risks of ASDs in human. These sporadic human genetic and mouse studies provide indirect, but important evidence that our predicted network module is potentially involved in pSHF for atrial septation. Importantly, our human ASD study identified the important role of a network module gene, Pcsk6, in the ASD onset. Pcsk6, also known as Pace4 or Spc4, is a member of subtilisin-like proprotein convertases. Pcsk6 have been shown to regulate transforming growth factor (TGF)-β-signaling network through cleavage of BMP4 and Nodal, and therefore maintains the balance between Nodal and bone morphogenetic protein (BMP) signaling during embryonic development (53,64–66). According to a previous report, Pcsk6 knockout mouse embryos displayed various heart defects including common atrium (53). The results of our DC analysis combined with in vivo expression and transactivation data of Pcsk6, and supported by this human study, suggested that Pcsk6 is a downstream target of Tbx5 involved in atrial septation. Considering the important roles of BMP signaling and TGF-β signaling in the development of atrial septum (66–69), it is possible Tbx5 interact with BMP or TGF-β signaling via Pcsk6 in the SHF for atrial septation. Future studies should focus on elucidating the important role and functional interactions between Tbx5 and Pcsk6 during atrial septation and on identifying whether this Spanish ASD family carries a mutation of Pcsk6.

Table 4.

Summary of network module genes in heart development

Gene ID Functions (heart development) Associated with CHD? Associated with ASD References 
Tbx5 Key transcription factor for heart development Yes Yes (8,19,21,24,44–49
Osr1 Atrial septation and valve development Yes Yes (8,24,25
Tll1 Required for atrial and ventricular normal septation and positioning of heart Yes Yes (50,51
Pcsk6 Axis formation Yes Yes (52,53
Hand2 Key transcription factor function in SHF for cardiogenesis; SHF marker Yes Unclear (54–58
Fgf8 Required in SHF; outflow tract formation and remodeling; left-right patterning Yes Unclear (Reviewed in 59
Dtna Involved in left ventricular trabeculation and atrial septation Yes Yes (60–62
Mapk7 Development of myocardium and coronary vasculature Yes Unclear (63
Angpt1 Angiogenesis Unclear Unclear  
Trex1 Unclear Unclear Unclear  
Gene ID Functions (heart development) Associated with CHD? Associated with ASD References 
Tbx5 Key transcription factor for heart development Yes Yes (8,19,21,24,44–49
Osr1 Atrial septation and valve development Yes Yes (8,24,25
Tll1 Required for atrial and ventricular normal septation and positioning of heart Yes Yes (50,51
Pcsk6 Axis formation Yes Yes (52,53
Hand2 Key transcription factor function in SHF for cardiogenesis; SHF marker Yes Unclear (54–58
Fgf8 Required in SHF; outflow tract formation and remodeling; left-right patterning Yes Unclear (Reviewed in 59
Dtna Involved in left ventricular trabeculation and atrial septation Yes Yes (60–62
Mapk7 Development of myocardium and coronary vasculature Yes Unclear (63
Angpt1 Angiogenesis Unclear Unclear  
Trex1 Unclear Unclear Unclear  

The resulting network hub contained 10 genes where Tbx5 directly connects with Osr1, Pcsk6 and Tll1, suggesting that Tbx5 regulates the pSHF development via multiple interacting pathways. It is noted that the penetration of ASD in Osr1+/−, Tbx5+/−, Tbx5+/−/Osr1+/− and Osr1−/− transgenic mouse is 0, 50, 100 and 75% (in mouse embryos surviving through E13.5), respectively (24). This penetration rank from the highest to the lowest did not inversely follow the level of Osr1: Osr1−/− mouse embryos have the lowest level of Osr1 in pSHF, but ASD incidence is lower, although not significant (P = 0.0547) than that of Tbx5+/−/Osr1+/− mouse embryos. This observation suggests that Tbx5 regulates multiple signaling pathways required for atrial septation. One of the pathways functions through Osr1, which has been molecularly and functionally validated to be required in pSHF for atrial septation (8,24).

Our pathway study suggested that Osr1 expression influenced several signaling pathways in pSHF, such as cell-cycle signaling, Hh-signaling pathway and Wnt-signaling pathway. It is defined that Tbx5–Osr1 interaction controls cell-cycle progression in pSHF for atrial septation (24). A Tbx5-Hh-signaling cascade is found to play a parallel role with Tbx5–Osr1 in the pSHF in that it regulates proliferation and cell-cycle progression during atrial septation (8,70). In addition, connection between Osr1- and Hh-signaling in pSHF for atrial septation has recently been suggested (24). Wnt signaling is known to play an important role in the SHF for heart development (9,71–74). Here, we showed that integrity of Wnt signaling is dependent on the level of Osr1 in pSHF. To be noted, except for Wnt2, expressions of all other tested genes involved in Wnt-signaling cascade were unchanged in either Tbx5 heterozygotes or Tbx5 and Osr1 compounded heterozygotes. Considering that Osr1 is located downstream of Tbx5, this result suggested that integrity of Wnt-signaling pathway is probably more sensitive to the level of Osr1 than to the level of Tbx5. A potential role of Osr1- and Wnt-signaling in pSHF for heart development needs to be further addressed.

In summary, this is a systems biology study aiming to elucidate the genetic mechanism underlying atrial septation and identify novel causative genes involved in ASDs. This study is the integration of mouse studies of the molecular and genetic mechanisms of heart development and human genetic studies combined with recently developed high-throughput bioinformatics technology and classical molecular approaches. The reported network module genes have been validated in pSHF of our transgenic mouse model of ASD and in human ASD study. Future study will focused on the functional study of the network module in atrial sepation.

Materials and Methods

Mouse lines

Mice were maintained in a mixed B6/129/SvEv background. The Tbx5 heterozygous mouse line was obtained from Dr Ivan Moskowitz (University of Chicago), and the Osr1+/− mouse line was obtained from the Jackson Laboratory. Mouse experiments were completed according to a protocol reviewed and approved by the Institutional Animal Care and Use Committee of the University of North Dakota in compliance with the US Public Health Service Policy on Humane Care and Use of Laboratory Animals.

Microdissection of pSHF and RNA extraction

To obtain the SHF splanchnic mesoderm for use in RNA-Seq and real-time PCR analyses, E9.5 embryos were dissected as described previously (8,24,70). Briefly, embryos were dissected above the outflow tract and below the inflow tract to obtain a thoracic section including the heart. The neural tube was removed by cutting through the foregut. Cutting between the outflow and inflow tracts separated the pSHF and aSHF. The heart, aSHF and pSHF were collected separately in RNA-later solution. RNA was extracted using a Qiagen RNeasy Mini Kit. DNA was removed from RNA samples by incubation with ribonuclease-free deoxyribonuclease I (RNase-free DNase I, Qiagen) at room temperature for 15 min.

Real-time PCR analysis

Two hundred nanograms of total RNA underwent reverse transcription by using a SuperScriptTM III Reverse Transcriptase kit (Invitrogen). Real-time PCR analysis was performed with a POWER SYBER Green PCR master mix (Applied Biosystems) and an iQ5 thermal cycler (Bio-Rad). The data were analyzed by using the delta–delta Ct (ΔΔCt) method with Gapdh as a normalization control (26).

RNA-Seq data analysis

Illumina next-generation sequencing was performed at University of Chicago Genomics Core (Chicago, USA) using Illumina HiSeq2000 instruments. RNA-Seq short reads were mapped to mouse reference genome GRCm38 using Bowtie2 package (75). The gene expression level was then quantified by Cufflink package (76). PCA and hierarchical clustering analysis (HCA) was used to assess the global gene-expression variations. The distance between samples in HCA was calculated by Pearson dissimilarity and the agglomeration of clusters was performed by the Ward linkage method. Regression analysis was used to identify the genes that were correlated with the expression of Osr1. The significantly correlated genes were selected by controlling the false discovery rate at 10%. The activated gene pathways were analyzed using a hypergeometric distribution-based method, which was developed in house using R programing language (77,78). The biochemical pathway information was obtained from the KEGG (www.genome.jp/kegg).

Gene-network analysis

A gene co-expression model was developed for regulatory networks analysis. The gene co-expression was quantified using a non-parametric distance metric, DC, which is able to account for both the linear and non-linear dependency between genes (37,38). The resulting adjacent matrix for the networks model was further analyzed using hierarchical clustering in order to identify network modules. Network modules consisted of genes sets with intensive interactions between genes. The resulting gene network model was compared with the 255 CHD-associated genes that were curated from literature (79). The network analysis was conducted using our in-house developed C++ and R programs.

Microarray and SNP array analysis

The Affymetrix microarray data for the P19CL6 cells were processed using the Bioconductor software packages. ANOVA was used to identify the genes that had significant transcriptional changes during P19CL6 differentiation. The gene–gene interaction was measured using DC tests. SNP array data for a Spanish family with six ASD patients and four healthy individuals were collected from whole blood samples using the Illumina Human660W-Quad chip. The SNPs were first tested for Hardy–Weinberg equilibrium in 46 healthy people. A familial analysis was conducted using a logarithm of the odds method. The SNP associated with ASD were further tested using χ2 tests.

Luciferase reporter assay

Regulatory regions were amplified and cloned upstream of a minimal promoter in pGL3 Basic vector (Promega). The pcDNA3-Tbx5 plasmid was obtained from Dr Xiaohua Zhang (Nemours Research Institute, Orlando, FL, USA). Reporter vector and pRL-TK (Promega) were co-transfected into HEK293T cells, with or without Tbx5 expression vector, using FuGene HD Transfection Reagent (Promega). Cells were lysed and assayed 24 h after transfection using the Dual-Luciferase Reporter Assay System (Promega). Firefly luciferase activity was normalized to Renilla luciferase activity followed by endogenous signal subtraction.

Site-directed mutagenesis

Mutant reporter vectors were generated by deleting Tbx5-binding sites using QuikChange Lightning Site-Directed Mutagenesis Kit (Agilent Technologies). Primers used are as follows: Pcsk6-Fr2m1 forward (5′-CAAACCTTACTTAATTCATCCACCCACACCTG-3′), reverse (5′-CAGGTGTGGGTGGATGAATTAAGTAAGGTTTG-3′); Pcsk6-Fr2m2 forward (5′-CTTACTTAATTCATCCACGAGCCCACAAGTTGAG-3′), reverse (5′-CTCAACTTGTGGGCTCGTGGATGAATTAAGTAAG-3′).

Chromatin immunoprecipitation

E9.5 embryos were micro-dissected in cold phosphate buffered saline (PBS) containing Protease Inhibitor Cocktail (Roche) to isolate the primitive heart region. Approximately 20 tissues were pooled into three samples. Tissues were cross-linked with 1% formaldehyde for 15 min at room temperature and terminated with glycine. After PBS washes, tissues were dissociated by shaking at 37°C for 1–2 h at 100 rpm in Collagenase, Type II (Gibco) solution. Sonication was performed using a Covaris S220 sonicator to generate average fragment size of 600 bp. Samples were incubated with Tbx5 antibody (Santa Cruz, sc-17866) overnight at 4°C then incubated with Dynabeads Protein G (Life Technologies) for 2 h, washed and reverse-cross-linked. Quantitative PCR data of regulatory region were expressed as fold enrichment relative to negative control region.

Immunohistochemical staining

E9.5 embryos were fixed in 4% paraformaldehyde overnight at 4°C, and then embedded in paraffin and sectioned sagittally at 5 μm. Immunohistochemistry was then performed on sections using anti-PACE4 antibody (1:50; Abcam, ab140934) and VECTASTAIN Elite ABC Kit (Vector Labs, PK-6101) per instructions.

Authors’ Contributions

L.X. and K.Z. designed the study; K.Z. and M.X. performed bioinformatics and statistical analysis; L.Z., M.X., J.L. Z.X. and N.C. performed experiments; D.H.S. and P.G. collected the human data; L.X., K.Z. and Q.W. wrote the paper.

Funding

This project was supported by grants from the National Institutes of Health (NIH-1R15HL117238 to L.X., National Center for Research Resources, 5P20RR016471-12/8 P20 GM103442-12 to L.X. and K.Z. and NIH-R01HL121358 and NIH-R01HL126729 to Q.W.) and the American Heart Association (13SDG14650009 to L.X. and 15GRNT25700195 to K.Z.).

Acknowledgements

We thank Drs Sergei Nechaev and Archana Dhasarathy for their help on ChIP-qPCR experiment.

Conflict of Interest statement. None declared.

References

1
Snarr
B.S.
,
Wirrig
E.E.
,
Phelps
A.L.
,
Trusk
T.C.
,
Wessels
A.
(
2007
)
A spatiotemporal evaluation of the contribution of the dorsal mesenchymal protrusion to cardiac development
.
Dev. Dyn.
 ,
236
,
1287
1294
.
2
Wessels
A.
,
Anderson
R.H.
,
Markwald
R.R.
,
Webb
S.
,
Brown
N.A.
,
Viragh
S.
,
Moorman
A.F.
,
Lamers
W.H.
(
2000
)
Atrial development in the human heart: an immunohistochemical study with emphasis on the role of mesenchymal tissues
.
Anat. Rec.
 ,
259
,
288
300
.
3
Anderson
R.H.
,
Webb
S.
,
Brown
N.A.
,
Lamers
W.
,
Moorman
A.
(
2003
)
Development of the heart: (2) Septation of the atriums and ventricles
.
Heart
 ,
89
,
949
958
.
4
Therrien
J.
,
Webb
G.
(
2003
)
Clinical update on adults with congenital heart disease
.
Lancet
 ,
362
,
1305
1313
.
5
Blom
N.A.
,
Ottenkamp
J.
,
Wenink
A.G.
,
Gittenberger-de Groot
A.C.
(
2003
)
Deficiency of the vestibular spine in atrioventricular septal defects in human fetuses with down syndrome
.
Am. J. Cardiol.
 ,
91
,
180
184
.
6
Goddeeris
M.M.
,
Rho
S.
,
Petiet
A.
,
Davenport
C.L.
,
Johnson
G.A.
,
Meyers
E.N.
,
Klingensmith
J.
(
2008
)
Intracardiac septation requires hedgehog-dependent cellular contributions from outside the heart
.
Development
 ,
135
,
1887
1895
.
7
Webb
S.
,
Brown
N.A.
,
Anderson
R.H.
(
1998
)
Formation of the atrioventricular septal structures in the normal mouse
.
Circ. Res.
 ,
82
,
645
656
.
8
Xie
L.
,
Hoffmann
A.D.
,
Burnicka-Turek
O.
,
Friedland-Little
J.M.
,
Zhang
K.
,
Moskowitz
I.P.
(
2012
)
Tbx5-hedgehog molecular networks are essential in the second heart field for atrial septation
.
Dev. Cell
 ,
23
,
280
291
.
9
Tian
Y.
,
Cohen
E.D.
,
Morrisey
E.E.
(
2010
)
The importance of Wnt signaling in cardiovascular development
.
Pediatr. Cardiol.
 ,
31
,
342
348
.
10
Snarr
B.S.
,
O'Neal
J.L.
,
Chintalapudi
M.R.
,
Wirrig
E.E.
,
Phelps
A.L.
,
Kubalak
S.W.
,
Wessels
A.
(
2007
)
Isl1 expression at the venous pole identifies a novel role for the second heart field in cardiac development
.
Circ. Res.
 ,
101
,
971
974
.
11
Mommersteeg
M.T.
,
Soufan
A.T.
,
de Lange
F.J.
,
van den Hoff
M.J.
,
Anderson
R.H.
,
Christoffels
V.M.
,
Moorman
A.F.
(
2006
)
Two distinct pools of mesenchyme contribute to the development of the atrial septum
.
Circ. Res.
 ,
99
,
351
353
.
12
Goddeeris
M.M.
,
Schwartz
R.
,
Klingensmith
J.
,
Meyers
E.N.
(
2007
)
Independent requirements for Hedgehog signaling by both the anterior heart field and neural crest cells for outflow tract development
.
Development
 ,
134
,
1593
1604
.
13
Galli
D.
,
Dominguez
J.N.
,
Zaffran
S.
,
Munk
A.
,
Brown
N.A.
,
Buckingham
M.E.
(
2008
)
Atrial myocardium derives from the posterior region of the second heart field, which acquires left-right identity as Pitx2c is expressed
.
Development
 ,
135
,
1157
1167
.
14
Verzi
M.P.
,
McCulley
D.J.
,
De Val
S.
,
Dodou
E.
,
Black
B.L.
(
2005
)
The right ventricle, outflow tract, and ventricular septum comprise a restricted expression domain within the secondary/anterior heart field
.
Dev. Biol.
 ,
287
,
134
145
.
15
Heidt
A.B.
,
Black
B.L.
(
2005
)
Transgenic mice that express Cre recombinase under control of a skeletal muscle-specific promoter from mef2c
.
Genesis
 ,
42
,
28
32
.
16
Dodou
E.
,
Verzi
M.P.
,
Anderson
J.P.
,
Xu
S.M.
,
Black
B.L.
(
2004
)
Mef2c is a direct transcriptional target of ISL1 and GATA factors in the anterior heart field during mouse embryonic development
.
Development
 ,
131
,
3931
3942
.
17
Lin
L.
,
Bu
L.
,
Cai
C.L.
,
Zhang
X.
,
Evans
S.
(
2006
)
Isl1 is upstream of sonic hedgehog in a pathway required for cardiac morphogenesis
.
Dev. Biol.
 ,
295
,
756
763
.
18
Hoffmann
A.D.
,
Peterson
M.A.
,
Friedland-Little
J.M.
,
Anderson
S.A.
,
Moskowitz
I.P.
(
2009
)
Sonic hedgehog is required in pulmonary endoderm for atrial septation
.
Development
 ,
136
,
1761
1770
.
19
Basson
C.T.
,
Bachinsky
D.R.
,
Lin
R.C.
,
Levi
T.
,
Elkins
J.A.
,
Soults
J.
,
Grayzel
D.
,
Kroumpouzou
E.
,
Traill
T.A.
,
Leblanc-Straceski
J.
et al
. (
1997
)
Mutations in human TBX5 [corrected] cause limb and cardiac malformation in Holt-Oram syndrome
.
Nat. Genet.
 ,
15
,
30
35
.
20
Garg
V.
,
Kathiriya
I.S.
,
Barnes
R.
,
Schluterman
M.K.
,
King
I.N.
,
Butler
C.A.
,
Rothrock
C.R.
,
Eapen
R.S.
,
Hirayama-Yamada
K.
,
Joo
K.
et al
. (
2003
)
GATA4 mutations cause human congenital heart defects and reveal an interaction with TBX5
.
Nature
 ,
424
,
443
447
.
21
Li
Q.Y.
,
Newbury-Ecob
R.A.
,
Terrett
J.A.
,
Wilson
D.I.
,
Curtis
A.R.
,
Yi
C.H.
,
Gebuhr
T.
,
Bullen
P.J.
,
Robson
S.C.
,
Strachan
T.
et al
. (
1997
)
Holt-Oram syndrome is caused by mutations in TBX5, a member of the Brachyury (T) gene family
.
Nat. Genet.
 ,
15
,
21
29
.
22
Schott
J.J.
,
Benson
D.W.
,
Basson
C.T.
,
Pease
W.
,
Silberbach
G.M.
,
Moak
J.P.
,
Maron
B.J.
,
Seidman
C.E.
,
Seidman
J.G.
(
1998
)
Congenital heart disease caused by mutations in the transcription factor NKX2-5
.
Science
 ,
281
,
108
111
.
23
Xie
L.
,
Weichel
B.
,
Ohm
J.E.
,
Zhang
K.
(
2011
)
An integrative analysis of DNA methylation and RNA-Seq data for human heart, kidney and liver
.
BMC Syst. Biol.
 ,
5
,
S4
.
24
Zhou
L.
,
Liu
J.
,
Olson
P.
,
Zhang
K.
,
Wynne
J.
,
Xie
L.
(
2015
)
Tbx5 and Osr1 interact to regulate posterior second heart field cell cycle progression for cardiac septation
.
J. Mol. Cell. Cardiol.
 ,
85
,
1
12
.
25
Wang
Q.
,
Lan
Y.
,
Cho
E.S.
,
Maltby
K.M.
,
Jiang
R.
(
2005
)
Odd-skipped related 1 (Odd 1) is an essential regulator of heart and urogenital development
.
Dev. Biol.
 ,
288
,
582
594
.
26
Schmittgen
T.D.
,
Livak
K.J.
(
2008
)
Analyzing real-time PCR data by the comparative C(T) method
.
Nat. Protoc.
 ,
3
,
1101
1108
.
27
Zhang
Z.
,
Cerrato
F.
,
Xu
H.
,
Vitelli
F.
,
Morishima
M.
,
Vincentz
J.
,
Furuta
Y.
,
Ma
L.
,
Martin
J.F.
,
Baldini
A.
et al
. (
2005
)
Tbx1 expression in pharyngeal epithelia is necessary for pharyngeal arch artery development
.
Development
 ,
132
,
5307
5315
.
28
Watanabe
Y.
,
Zaffran
S.
,
Kuroiwa
A.
,
Higuchi
H.
,
Ogura
T.
,
Harvey
R.P.
,
Kelly
R.G.
,
Buckingham
M.
(
2012
)
Fibroblast growth factor 10 gene regulation in the second heart field by Tbx1, Nkx2-5, and Islet1 reveals a genetic switch for down-regulation in the myocardium
.
Proc. Natl Acad. Sci. USA
 ,
109
,
18273
18280
.
29
Rana
M.S.
,
Theveniau-Ruissy
M.
,
De Bono
C.
,
Mesbah
K.
,
Francou
A.
,
Rammah
M.
,
Dominguez
J.N.
,
Roux
M.
,
Laforest
B.
,
Anderson
R.H.
et al
. (
2014
)
Tbx1 coordinates addition of posterior second heart field progenitor cells to the arterial and venous poles of the heart
.
Circ. Res.
 ,
115
,
790
799
.
30
Cao
J.
,
Zhang
S.
(
2014
)
A Bayesian extension of the hypergeometric test for functional enrichment analysis
.
Biometrics
 ,
70
,
84
94
.
31
He
A.
,
Kong
S.W.
,
Ma
Q.
,
Pu
W.T.
(
2011
)
Co-occupancy by multiple cardiac transcription factors identifies transcriptional enhancers active in heart
.
Proc. Natl Acad. Sci. USA
 ,
108
,
5632
5637
.
32
Clifford
R.L.
,
Deacon
K.
,
Knox
A.J.
(
2008
)
Novel regulation of vascular endothelial growth factor-A (VEGF-A) by transforming growth factor (beta)1: requirement for Smads, (beta)-CATENIN, AND GSK3(beta)
.
J. Biol. Chem.
 ,
283
,
35337
35353
.
33
Xu
X.
,
Mao
W.
,
Chen
Q.
,
Zhuang
Q.
,
Wang
L.
,
Dai
J.
,
Wang
H.
,
Huang
Z.
(
2014
)
Endostar, a modified recombinant human endostatin, suppresses angiogenesis through inhibition of Wnt/beta-catenin signaling pathway
.
PLoS One
 ,
9
,
e107463
.
34
Huang
C.L.
,
Liu
D.
,
Nakano
J.
,
Ishikawa
S.
,
Kontani
K.
,
Yokomise
H.
,
Ueno
M.
(
2005
)
Wnt5a expression is associated with the tumor proliferation and the stromal vascular endothelial growth factor—an expression in non-small-cell lung cancer
.
J. Clin. Oncol.
 ,
23
,
8765
8773
.
35
Lindsley
R.C.
,
Gill
J.G.
,
Kyba
M.
,
Murphy
T.L.
,
Murphy
K.M.
(
2006
)
Canonical Wnt signaling is required for development of embryonic stem cell-derived mesoderm
.
Development
 ,
133
,
3787
3796
.
36
Gradl
D.
,
Kuhl
M.
,
Wedlich
D.
(
1999
)
The Wnt/Wg signal transducer beta-catenin controls fibronectin expression
.
Mol. Cell. Biol.
 ,
19
,
5576
5587
.
37
Szekely
G.J.
,
Rizzo
M.L.
(
2009
)
Brownian distance covariance
.
Ann. Appl. Stat.
 ,
3
,
30
.
38
Szekely
G.J.
,
Rizzo
M.
,
Bakirov
N.K.
(
2007
)
Measuring and testing dependence by correlation of distance
.
Ann. Stat.
 ,
35
,
36
.
39
Lage
K.
,
Greenway
S.C.
,
Rosenfeld
J.A.
,
Wakimoto
H.
,
Gorham
J.M.
,
Segre
A.V.
,
Roberts
A.E.
,
Smoot
L.B.
,
Pu
W.T.
,
Pereira
A.C.
et al
. (
2012
)
Genetic and environmental risk factors in congenital heart disease functionally converge in protein networks driving heart development
.
Proc. Natl Acad. Sci. USA
 ,
109
,
14035
14040
.
40
Barth
J.L.
,
Clark
C.D.
,
Fresco
V.M.
,
Knoll
E.P.
,
Lee
B.
,
Argraves
W.S.
,
Lee
K.H.
(
2010
)
Jarid2 is among a set of genes differentially regulated by Nkx2.5 during outflow tract morphogenesis
.
Dev. Dyn.
 ,
239
,
2024
2033
.
41
Stuart
J.M.
,
Segal
E.
,
Koller
D.
,
Kim
S.K.
(
2003
)
A gene-coexpression network for global discovery of conserved genetic modules
.
Science
 ,
302
,
249
255
.
42
Carter
S.L.
,
Brechbuhler
C.M.
,
Griffin
M.
,
Bond
A.T.
(
2004
)
Gene co-expression network topology provides a framework for molecular characterization of cellular state
.
Bioinformatics
 ,
20
,
2242
2250
.
43
Zhang
B.
,
Horvath
S.
(
2005
)
A general framework for weighted gene co-expression network analysis
.
Stat. Appl. Genet. Mol. Biol.
 ,
4
,
Article17
.
44
Bruneau
B.G.
,
Nemer
G.
,
Schmitt
J.P.
,
Charron
F.
,
Robitaille
L.
,
Caron
S.
,
Conner
D.A.
,
Gessler
M.
,
Nemer
M.
,
Seidman
C.E.
,
Seidman
J.G.
(
2001
)
A murine model of Holt-Oram syndrome defines roles of the T-box transcription factor Tbx5 in cardiogenesis and disease
.
Cell
 ,
106
,
709
721
.
45
Bruneau
B.G.
,
Logan
M.
,
Davis
N.
,
Levi
T.
,
Tabin
C.J.
,
Seidman
J.G.
,
Seidman
C.E.
(
1999
)
Chamber-specific cardiac expression of Tbx5 and heart defects in Holt-Oram syndrome
.
Dev. Biol.
 ,
211
,
100
108
.
46
Bruneau
B.G.
(
2003
)
The developing heart and congenital heart defects: a make or break situation
.
Clin. Genet.
 ,
63
,
252
261
.
47
Takeuchi
J.K.
,
Bruneau
B.G.
(
2009
)
Directed transdifferentiation of mouse mesoderm to heart tissue by defined factors
.
Nature
 ,
459
,
708
711
.
48
Moskowitz
I.P.
,
Wang
J.
,
Peterson
M.A.
,
Pu
W.T.
,
Mackinnon
A.C.
,
Oxburgh
L.
,
Chu
G.C.
,
Sarkar
M.
,
Berul
C.
,
Smoot
L.
et al
. (
2011
)
Transcription factor genes Smad4 and Gata4 cooperatively regulate cardiac valve development. [corrected]
.
Proc. Natl Acad. Sci. USA
 ,
108
,
4006
4011
.
49
Mori
A.D.
,
Zhu
Y.
,
Vahora
I.
,
Nieman
B.
,
Koshiba-Takeuchi
K.
,
Davidson
L.
,
Pizard
A.
,
Seidman
J.G.
,
Seidman
C.E.
,
Chen
X.J.
et al
. (
2006
)
Tbx5-dependent rheostatic control of cardiac gene expression and morphogenesis
.
Dev. Biol.
 ,
297
,
566
586
.
50
Stanczak
P.
,
Witecka
J.
,
Szydlo
A.
,
Gutmajster
E.
,
Lisik
M.
,
Augusciak-Duma
A.
,
Tarnowski
M.
,
Czekaj
T.
,
Czekaj
H.
,
Sieron
A.L.
(
2009
)
Mutations in mammalian tolloid-like 1 gene detected in adult patients with ASD
.
Eur. J. Hum. Genet.
 ,
17
,
344
351
.
51
Clark
T.G.
,
Conway
S.J.
,
Scott
I.C.
,
Labosky
P.A.
,
Winnier
G.
,
Bundy
J.
,
Hogan
B.L.
,
Greenspan
D.S.
(
1999
)
The mammalian Tolloid-like 1 gene, Tll1, is necessary for normal septation and positioning of the heart
.
Development
 ,
126
,
2631
2642
.
52
Flaquer
A.
,
Baumbach
C.
,
Pinero
E.
,
Garcia Algas
F.
,
de la Fuente Sanchez
M.A.
,
Rosell
J.
,
Toquero
J.
,
Alonso-Pulpon
L.
,
Garcia-Pavia
P.
,
Strauch
K.
et al
. (
2013
)
Genome-wide linkage analysis of congenital heart defects using MOD score analysis identifies two novel loci
.
BMC Genet.
 ,
14
,
44
.
53
Constam
D.B.
,
Robertson
E.J.
(
2000
)
SPC4/PACE4 regulates a TGFbeta signaling network during axis formation
.
Genes Dev.
 ,
14
,
1146
1155
.
54
Tamura
M.
,
Hosoya
M.
,
Fujita
M.
,
Iida
T.
,
Amano
T.
,
Maeno
A.
,
Kataoka
T.
,
Otsuka
T.
,
Tanaka
S.
,
Tomizawa
S.
,
Shiroishi
T.
(
2013
)
Overdosage of Hand2 causes limb and heart defects in the human chromosomal disorder partial trisomy distal 4q
.
Hum. Mol. Genet.
 ,
22
,
2471
2481
.
55
Tsuchihashi
T.
,
Maeda
J.
,
Shin
C.H.
,
Ivey
K.N.
,
Black
B.L.
,
Olson
E.N.
,
Yamagishi
H.
,
Srivastava
D.
(
2011
)
Hand2 function in second heart field progenitors is essential for cardiogenesis
.
Dev. Biol.
 ,
351
,
62
69
.
56
Xu
W.
,
Ahmad
A.
,
Dagenais
S.
,
Iyer
R.K.
,
Innis
J.W.
(
2012
)
Chromosome 4q deletion syndrome: narrowing the cardiovascular critical region to 4q32.2-q34.3
.
Am. J. Med. Genet.
 ,
158A
,
635
640
.
57
Barnes
R.M.
,
Firulli
B.A.
,
VanDusen
N.J.
,
Morikawa
Y.
,
Conway
S.J.
,
Cserjesi
P.
,
Vincentz
J.W.
,
Firulli
A.B.
(
2011
)
Hand2 loss-of-function in Hand1-expressing cells reveals distinct roles in epicardial and coronary vessel development
.
Circ. Res.
 ,
108
,
940
949
.
58
McFadden
D.G.
,
Barbosa
A.C.
,
Richardson
J.A.
,
Schneider
M.D.
,
Srivastava
D.
,
Olson
E.N.
(
2005
)
The Hand1 and Hand2 transcription factors regulate expansion of the embryonic cardiac ventricles in a gene dosage-dependent manner
.
Development
 ,
132
,
189
201
.
59
Zaffran
S.
,
Kelly
R.G.
(
2012
)
New developments in the second heart field
.
Differentiation
 ,
84
,
17
24
.
60
Chen
C.P.
,
Huang
M.C.
,
Chen
Y.Y.
,
Chern
S.R.
,
Wu
P.S.
,
Chen
Y.T.
,
Su
J.W.
,
Wang
W.
(
2013
)
Prenatal diagnosis of de novo interstitial deletions involving 5q23.1-q23.3 and 18q12.1-q12.3 by array CGH using uncultured amniocytes in a pregnancy with fetal interrupted aortic arch and atrial septal defect
.
Gene
 ,
531
,
496
501
.
61
Chang
B.
,
Nishizawa
T.
,
Furutani
M.
,
Fujiki
A.
,
Tani
M.
,
Kawaguchi
M.
,
Ibuki
K.
,
Hirono
K.
,
Taneichi
H.
,
Uese
K.
et al
. (
2011
)
Identification of a novel TPM1 mutation in a family with left ventricular noncompaction and sudden death
.
Mol. Genet. Metab.
 ,
102
,
200
206
.
62
Xia
S.
,
Wang
H.
,
Zhang
X.
,
Zhu
J.
,
Tang
X.
(
2008
)
Clinical presentation and genetic analysis of a five generation Chinese family with isolated left ventricular noncompaction
.
Intern. Med.
 ,
47
,
577
583
.
63
Bhattacharya
S.
,
Macdonald
S.T.
,
Farthing
C.R.
(
2006
)
Molecular mechanisms controlling the coupled development of myocardium and coronary vasculature
.
Clin. Sci.
 ,
111
,
35
46
.
64
Birsoy
B.
,
Berg
L.
,
Williams
P.H.
,
Smith
J.C.
,
Wylie
C.C.
,
Christian
J.L.
,
Heasman
J.
(
2005
)
XPACE4 is a localized pro-protein convertase required for mesoderm induction and the cleavage of specific TGFbeta proteins in Xenopus development
.
Development
 ,
132
,
591
602
.
65
Constam
D.B.
,
Robertson
E.J.
(
1999
)
Regulation of bone morphogenetic protein activity by pro domains and proprotein convertases
.
J. Cell .Biol.
 ,
144
,
139
149
.
66
Beck
S.
,
Le Good
J.A.
,
Guzman
M.
,
Ben Haim
N.
,
Roy
K.
,
Beermann
F.
,
Constam
D.B.
(
2002
)
Extraembryonic proteases regulate Nodal signalling during gastrulation
.
Nat. Cell. Biol.
 ,
4
,
981
985
.
67
McCulley
D.J.
,
Kang
J.O.
,
Martin
J.F.
,
Black
B.L.
(
2008
)
BMP4 is required in the anterior heart field and its derivatives for endocardial cushion remodeling, outflow tract septation and semilunar valve development
.
Dev. Dyn.
 ,
237
,
3200
3209
.
68
Briggs
L.E.
,
Phelps
A.L.
,
Brown
E.
,
Kakarla
J.
,
Anderson
R.H.
,
van den Hoff
M.J.
,
Wessels
A.
(
2013
)
Expression of the BMP receptor Alk3 in the second heart field is essential for development of the dorsal mesenchymal protrusion and atrioventricular septation
.
Circ. Res.
 ,
112
,
1420
1432
.
69
Nakajima
Y.
,
Yamagishi
T.
,
Nakamura
H.
,
Markwald
R.R.
,
Krug
E.L.
(
1998
)
An autocrine function for transforming growth factor (TGF)-beta3 in the transformation of atrioventricular canal endocardium into mesenchyme during chick heart development
.
Dev. Biol.
 ,
194
,
99
113
.
70
Hoffmann
A.D.
,
Yang
X.H.
,
Burnicka-Turek
O.
,
Bosman
J.D.
,
Ren
X.
,
Steimle
J.D.
,
Vokes
S.A.
,
McMahon
A.P.
,
Kalinichenko
V.V.
,
Moskowitz
I.P.
(
2014
)
Foxf genes integrate tbx5 and hedgehog pathways in the second heart field for cardiac septation
.
PLoS Genet.
 ,
10
,
e1004604
.
71
Gessert
S.
,
Kuhl
M.
(
2010
)
The multiple phases and faces of Wnt signaling during cardiac differentiation and development
.
Circ. Res.
 ,
107
,
186
199
.
72
Ai
D.
,
Fu
X.
,
Wang
J.
,
Lu
M.F.
,
Chen
L.
,
Baldini
A.
,
Klein
W.H.
,
Martin
J.F.
(
2007
)
Canonical Wnt signaling functions in second heart field to promote right ventricular growth
.
Proc. Natl Acad. Sci. USA
 ,
104
,
9319
9324
.
73
Klaus
A.
,
Saga
Y.
,
Taketo
M.M.
,
Tzahor
E.
,
Birchmeier
W.
(
2007
)
Distinct roles of Wnt/beta-catenin and Bmp signaling during early cardiogenesis
.
Proc. Natl Acad. Sci. USA
 ,
104
,
18531
18536
.
74
Cohen
E.D.
,
Wang
Z.
,
Lepore
J.J.
,
Lu
M.M.
,
Taketo
M.M.
,
Epstein
D.J.
,
Morrisey
E.E.
(
2007
)
Wnt/beta-catenin signaling promotes expansion of Isl-1-positive cardiac progenitor cells through regulation of FGF signaling
.
J. Clin. Inv.
 ,
117
,
1794
1804
.
75
Langmead
B.
,
Trapnell
C.
,
Pop
M.
,
Salzberg
S.L.
(
2009
)
Ultrafast and memory-efficient alignment of short DNA sequences to the human genome
.
Genome Biol.
 ,
10
,
R25
.
76
Trapnell
C.
,
Roberts
A.
,
Goff
L.
,
Pertea
G.
,
Kim
D.
,
Kelley
D.R.
,
Pimentel
H.
,
Salzberg
S.L.
,
Rinn
J.L.
,
Pachter
L.
(
2012
)
Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks
.
Nat. Protoc.
 ,
7
,
562
578
.
77
Zhang
K.
,
Wang
H.
,
Bathke
A.C.
,
Harrar
S.W.
,
Piepho
H.P.
,
Deng
Y.
(
2011
)
Gene set analysis for longitudinal gene expression data
.
BMC Bioinformatics
 ,
12
,
273
.
78
Garrett
S.H.
,
Clarke
K.
,
Sens
D.A.
,
Deng
Y.
,
Somji
S.
,
Zhang
K.K.
(
2013
)
Short and long term gene expression variation and networking in human proximal tubule cells when exposed to cadmium
.
BMC Med. Genomics
 ,
6
(
Suppl 1
),
S2
.
79
Lage
K.
,
Mollgard
K.
,
Greenway
S.
,
Wakimoto
H.
,
Gorham
J.M.
,
Workman
C.T.
,
Bendsen
E.
,
Hansen
N.T.
,
Rigina
O.
,
Roque
F.S.
et al
. (
2010
)
Dissecting spatio-temporal protein networks driving human heart development and related disorders
.
Mol. Syst. Biol.
 ,
6
,
381
.