Heart morphogenesis gene regulatory networks revealed by temporal expression analysis

Abstract

During embryogenesis the heart forms as a linear tube that then undergoes multiple simultaneous morphogenetic events to obtain its mature shape. To understand the gene regulatory networks (GRNs) driving this phase of heart development, during which many congenital heart disease malformations likely arise, we conducted an RNA-seq timecourse in zebrafish from 30 hpf to 72 hpf and identified 5861 genes with altered expression. We clustered the genes by temporal expression pattern, identified transcription factor binding motifs enriched in each cluster, and generated a model GRN for the major gene batteries in heart morphogenesis. This approach predicted hundreds of regulatory interactions and found batteries enriched in specific cell and tissue types, indicating that the approach can be used to narrow the search for novel genetic markers and regulatory interactions. Subsequent analyses confirmed the GRN using two mutants, Tbx5 and nkx2-5, and identified sets of duplicated zebrafish genes that do not show temporal subfunctionalization. This dataset provides an essential resource for future studies on the genetic/epigenetic pathways implicated in congenital heart defects and the mechanisms of cardiac transcriptional regulation.

Keywords: Gene regulatory network; Heart development; RNA-seq; Timecourse; Zebrafish.

Fig. 1.

RNA-seq timecourse analysis during heart looping morphogenesis. (A) Schematic of zebrafish heart (green) looping during the time period covered by the timecourse. Although several overlapping morphogenetic events are occurring, including cardiomyocyte maturation, initial trabeculation and sino-atrial (SA) node, atrioventricular canal (AVC) and valve formation, animals were staged based on heart looping. (B) Multidimensional scaling (MDS) of the RNA-seq samples. Relative distances between samples indicate their relative similarity (closer indicates more similar). Letters indicate the replicate (A, B or C); numbers (and colors) indicate the hours post-fertilization that the sample was collected. Each sample contained pooled hearts from ∼200 embryos. Outliers A30 and B42 were excluded from subsequent analysis. (C) Volcano plot showing the maximum log₂-scaled fold change of any time point compared with 30 hpf on the x-axis and the Phred-scaled P-value of the negative binomial likelihood ratio test, which tests for differential expression anywhere in the timecourse, on the y-axis. The blue horizontal line shows a Phred-scaled P-value cutoff of 13 (equivalent to a P-value of 0.05). The green vertical lines indicate a log₂ fold change cutoff of 1 and −1. Further analyses were conducted only on the differentially expressed genes in the upper left and upper right quadrants. (D) Heat map and hierarchical clustering of genes showing statistically significant changes in gene expression over the timecourse. Red indicates high expression and blue indicates low expression. Replicates were first averaged to create one column per time point.

Fig. 2.

Self-organizing map (SOM) analysis of differentially expressed genes during heart looping. (A) SOM results for the timecourse data. Each panel represents one cluster identified by SOM analysis. The red line indicates the median expression for the pattern, while the blue and gray regions represent the interquartile range (IQR) and the range of the expression levels, respectively. Numbers in the lower right corner indicate the number of genes assigned to each SOM cluster. (B) Graph of interactions found in the GEA_CLR database from the UCSC Interaction Browser between two genes within cluster A1. Orange lines indicate an activating interaction, blue lines a repressing interaction, and gray lines indicate unknown interaction types. (C) Bootstrap analysis of the assortativity coefficients using randomized cluster assignments (10,000 replicates). The vertical red line indicates the assortativity coefficient for the actual results, which was greater than any of the 10,000 bootstrap replicates.

Fig. 3.

Transcription factors binding motifs enriched in specific SOM clusters. (A) Transcription factor binding motifs enriched in one or more clusters. Gene regulatory interactions predicted by HOMER and found in the GEA_CLR database (‘known') or not found in the GEA_CLR database (‘novel') are shown, as well as the percentage of interactions that are novel. (B) Enriched transcription factor motifs and their respective SOM clusters.

Fig. 4.

Characteristics of enriched transcription factor binding motifs. (A) The number of transcription factor binding motifs enriched in each cluster. (B) The number of clusters that each transcription factor binding motif is enriched in. (C) The ten largest gene batteries. Genes were grouped by shared complement of transcription factor binding sites enriched in the first 1 kb upstream of the transcription start site.

Fig. 5.

Distribution of enriched transcription factor motifs across the SOM. (A) Enriched transcription factor binding motifs in their corresponding SOM cluster locations. Black font indicates transcription factor binding motifs enriched exclusively in one SOM cluster. Colored transcription factor binding motifs are found in more than one cluster. (B) Correspondence of Sox gene expression (diagonal lines) and Sox binding motif enrichment (orange shading) in the SOM clusters. (C) Correspondence of Klf gene expression (diagonal lines) and Klf motif enrichment (green shading) in the SOM clusters.

Fig. 6.

GRNs for cluster A1. (A) Gene regulatory interaction graph for SOM cluster A1 with Tbx transcription factors added. (B) Hierarchical clustering dendrogram based on vertex similarity in the graph shown in A. Putative batteries with identical edges in the graph are marked by a red line.

Fig. 7.

SOM cluster assignment of genes with altered expression in Tbx5 mutant mice and nkx2-5 mutant zebrafish. (A) Comparison of genes in cluster A1 predicted to have Tbx binding sites with those differentially expressed in the Tbx5 null mouse. (B) Location of genes with altered expression in nkx2-5 null zebrafish hearts at 48 hpf in the SOM.

Fig. 8.

Distribution of genes involved in specific processes/cell types within the SOM. A curated list of annotations for genes known to be involved in heart development was generated from the literature. Each panel represents a single annotation, and the location of shaded rectangles within the panel corresponds to SOM clusters where the genes with that annotation are located. The color of the shading indicates the number of genes for the given annotation found in the corresponding SOM cluster.