In vivo injury and regeneration of the zebrafish heart

Using a cryoinjury procedure26, we induced damage in adult zebrafish hearts. To dynamically monitor the regeneration process, we recovered samples at different post-injury times: 4 hours (hpi), 1, 3, 7, 14 and 90 days (dpi). Injured hearts were compared to healthy hearts from control fish in 3 independent experiments. We extracted RNA from heart ventricles for microarray experiments. We also recovered whole hearts to visualize healthy cardiomyocytes, apoptosis and fibrotic scar formation (Fig. 2). Healthy cardiomyocytes are observed in the whole ventricle and the atrium of control hearts, and injured hearts are totally devoid of such cells at 3 dpi. However, the size of the injury decreases all along the regeneration process, while new healthy cardiomyocytes are added from the border zone of the injured region, until 90 dpi when there are no visible (anatomical/cellular) differences with controls. Whereas massive cell death (in green, Fig. 2) is visible at the injury site as soon as 4 hpi and until 1 dpi, control hearts lack apoptotic cells. Lastly, while blood accumulates in the infarcted area at 1 dpi, formation of the fibrotic scar (in blue) is observed in the injury site at 3 dpi. Scar size decreases during the regeneration process and the injured area is gradually occupied by newly formed cardiomyocytes. At 90 dpi, ventricles are indistinguishable from controls. These results are concordant with previous work5,25 and corroborate the relevance of our model.

Figure 2: The different stages of heart regeneration in the zebrafish. Sagittal sections of adult zebrafish heart: anterior is towards the top and ventral towards the right. Hearts were cryoinjured and recovered at the indicated time points post-injury, and compared to healthy hearts from control fish. Sections were immunostained for tropomyosin in red, or TUNEL in green (nuclei are stained with DAPI in blue). Fibrosis was monitored by Masson-Goldner trichrome staining: healthy myocardial tissue (in red) and fibrotic areas (in blue). A: atrium; B: bulbus arteriosus; V: ventricle; hpi: hours post-injury; dpi: days post-injury. The injured area is indicated by a dotted circle. Scale bar is 100 μm. Full size image

Time-specific changes of gene expression during heart regeneration

Using whole-genome microarrays, we obtained expression profiles at 4 hpi, 1, 3, 7, 14 and 90 dpi. Independent measurements were obtained in triplicates at each time point, and from 3 control samples. The differential expression of genes at each time in relation to controls was statistically estimated. The largest numbers of differentially-expressed genes (DEGs, FDR < 0.001) were obtained at the early stages of regeneration (from 4 hpi to 3 dpi). Although at later times the number of DEGs progressively decreased, thousands of statistically detectable changes were still observed (FDR < 0.001, Fig. 3A). The DEGs obtained at each time are enriched in functional annotations relevant to heart regeneration, such as apoptosis and angiogenesis (Fig. 3B, Supplementary Table S1). As early as 4 hpi, enrichments in genes involved in apoptosis (consistent with our staining, Fig. 2), angiogenesis and cell migration are statistically detectable. Previous work by our group has also reported angiogenesis from 3 days post-injury onwards. The expression of angiogenesis-involved genes precedes the development of new vessels. The fact that this activation can occur early cannot be explained solely by the time needed for gene transcription and translation. Notably, many angiogenic marker genes are also expressed in circulating hematopoietic cells. This is the case, for example for fli1a, kdlr or erg27,28. Given the fact that the infiltration of inflammatory cells is an early response to injury, it is possible that hematopoietic cells expressing endothelial cell markers can also home to the injured hearts.

Figure 3: Significant changes in gene expression during heart regeneration. (A) Numbers of DEGs at each time point. (B) Summary of statistically enriched functional terms of DEGs at each time point. (C) Visualization of control and post-injury samples based on PCA (genes with FDR < 0.001). (D) Hierarchical clustering of samples based on gene expression data (genes with FDR < 0.001). Full size image

At 1 dpi, transcriptomic alterations mostly impact genes implicated in energy metabolism, amino-acid biosynthesis and DNA replication, which could be an indication of enhanced cell proliferation. Changes are also observed in genes involved in proteolytic activities, indicating the beginning of the regeneration process29,30. At 3 and 7 dpi, the regeneration activity is boosted, as shown by the enrichment in peptidase activity, together with processes linked to cell proliferation such as DNA metabolism and replication. Also the extracellular matrix is highly implicated at 3 dpi, when fibrosis becomes visible in our histology staining. Processes related to cell adhesion are mainly enriched at later times. Altogether these results reflect crucial steps of the heart regeneration process.

By applying principal component analysis (PCA) to our gene expression data (DEGs, FDR < 0.001), it was possible to visualize well-defined time-specific clusters that clearly mirrored the ordered sequence of regeneration events (Fig. 3C). A hierarchical clustering of the data also distinctly separated samples according to expression patterns and times (Fig. 3D). These results indicate that our gene expression data reflects biologically-meaningful dynamic changes associated with heart regeneration.

A co-expression network characterizes dynamic heart regeneration states

To characterize heart regeneration at a systems-level, we focused on significant transcriptional changes that distinguish injured from healthier hearts. This analysis was also needed to reduce the large number of expressed genes and facilitate network generation and analysis. First, we identified DEGs in the “injured group” (samples at 4 hpi, 1, 3, 7 and 14 dpi) vs. those in the “healthier group” (control and 90 dpi regenerated hearts). Because 90 dpi samples are anatomically very similar to healthy hearts, and at 90 dpi the hearts are actually at the end stages of full heart regeneration, their assignment to the “healthier group” is both reasonable and supported by our analysis. This procedure and the processing of genes with multiple probes resulted in a set of 3467 genes (FDR < 0.005) that were considered for subsequent analyses. Note that, although overlaps are observed, this set of genes differs from those represented in Fig. 3A, which only correspond to DEGs in the specific times vs. controls.

Next, we computed the expression (Pearson) correlations between these genes, and focused on gene associations with relatively high (absolute) correlation coefficients, │r│>0.8. This filtering was necessary to reduce potential spurious associations and to focus on those correlations likely to be biologically informative. This selection allowed us to achieve a balance between network scale-free fit and connectivity properties as recommended previously31. Moreover, this choice made subsequent network visualization interpretable and annotation tasks manageable.

The combination of all the selected gene-gene associations resulted in a global co-expression network of heart regeneration. The network consists of 3467 genes (nodes) and 436,803 associations (edges). Color-coding of the nodes on the basis of their expression fold-changes (in relation to controls) gives an overall view of the dynamic changes of the network at different times (Supplementary Fig. S1). This highlights systems-level response patterns that underlie the regeneration process: from prominent changes at the early stages of regeneration to more subtle, fine-tuned responses at later stages. It also allows us to appreciate the gradual regression from massive heart damage to a network state that resembles that exhibited by control samples (Supplementary Movie S1).

Network modularity is linked to heart regeneration responses

Biological networks can be organized into modules of highly interconnected genes, which are implicated in the same biological processes or associated with specialized functions32. Here we identified modules of highly co-expressed genes through the application of network clustering algorithms. We focused on two techniques with well-established module detection capacity: WGCNA (Weighted Gene Co-Expression Network Analysis)33 and ClusterONE (Clustering with Overlapping Neighborhood Expansion)34. These techniques have enabled several biological network investigations, and represent complementary approaches to detecting network modules. We applied WGCNA to the global co-expression network and detected 17 modules (modules 1A to 17A), with sizes ranging from 38 (module 17A) to 707 genes (module 1A). With ClusterONE we identified 11 statistically significant modules (modules 1B to 11B), with sizes ranging from 8 (module 11B) to 491 genes (module 1B). A pairwise comparison of WGCNA and ClusterONE modules corroborated that the two techniques offer partially overlapping, complementary views of modularity (Supplementary Fig. S2). The largest overlap is represented by modules 11A (WGCNA) and 6B (ClusterONE) with Jaccard similarity coefficient of 0.706 (Supplementary Fig. S2). Apart from the expected high intra-module co-expression, there is a diversity of strong inter-module associations (Fig. 4, online resource). Moreover, a GO enrichment analysis shows that many of these modules (12 WGCNA and 4 ClusterONE modules) are significantly associated with biological processes relevant to heart regeneration (FDR < 0.05) (Figs 4 and S3, numbers of associations indicated in color legend). This includes cardiac cell differentiation, migration and embryonic development (Fig. 5).

Figure 4: Modular architecture of the gene co-expression network in the zebrafish heart regeneration. Circular plots of modules detected with WGCNA: Internal links (grey) represent the intra- and inter-module connectivity, whereas the color of the outer bar represents the number of functional terms significantly enriched in a given module (Fig. 5). Numbers shown next to the colored circles (left panel) indicate the number of functional terms associated with each color in the plot. Full size image

Figure 5: Summary of functional enrichments of modules. Dotted line represents threshold of statistical significance at FDR = 0.05 (Supplementary Methods). Full size image

To further determine the biological meaning of these results, we looked deeper into the gene composition and connectivity of these modules. Here we concentrate on module 14A because it is highly interconnected with other modules (Fig. 4) and consists of 61 highly co-expressed genes. This module also includes several genes with multiple connections with genes in other modules (Supplementary Fig. S4), and displays one of the largest numbers of statistically enriched biological processes implicated in heart regeneration, including muscle cell differentiation and cardiac cell differentiation (Fig. 5). It is highly enriched in processes related to embryogenesis, such as “embryonic development and morphogenesis” and “cell fate commitment and specification”. This indicates that genes in module 14A may play critical regeneration regulatory roles, promoting reactivation of the embryonic program which leads to dedifferentiation and further cell proliferation and differentiation (Fig. 5). Different genes found in Module 14A, such as nkx2.5 and csrp1, are indeed key regulators of cardiac or vascular embryogenesis. The transcription factor nkx2.5 (Supplementary Fig. S4) initiates the cardiogenic differentiation program in zebrafish35. This gene is a marker of cardiac progenitor cells and is re-expressed at the resection plane following ventricular amputation in the adult zebrafish14. Likewise csrp1, a member of the Wnt pathway, coordinates cardiac mesoderm cell migration during zebrafish embryonic development and its inactivation leads to cardiac bifida36.

Other genes such as tbx5 and ctgf, may also be relevant to cardiac regeneration. Indeed, tbx5, a transcription factor naturally expressed in the developing heart, is required (together with gata4 and mef2c) to allow reprogramming of cardiac fibroblasts into cardiomyocytes37. Whereas ctgf is a regulator of fibrosis during maladaptive remodeling in humans38. Other genes in module 14A have not yet been implicated in cardiac regeneration. They include dyrk2, LOC100535315 (Supplementary Fig. S4), and zgc:110366 which is still uncharacterized. The dyrk2 kinase negatively regulates cardiomyocyte growth39, and thus may play an important role in the restoration of the correct organ size during regeneration. Future research will be required to validate the relevance of such genes in cardiac regeneration.

Network hubs play controlling roles in heart regeneration

The identification of network hubs is useful to understand network function. Hubs may represent genes with influential biological roles or regulatory activity32. In the case of gene co-expression networks, hubs are genes exhibiting a statistically significant number of strong connections with other genes in the network. We applied the WiPer technique40 to identify hubs in our gene co-expression network of heart regeneration, and identified 425 genes that display a statistically-detectable large number of strong connections (adjusted-P < 0.05). Examination of top hubs located in different network modules highlighted the diversity of molecular functions potentially triggered or mediated by these genes (online resource). Il6st (also known as gp130 receptor), which promotes the differentiation of embryonic stem cells into cardiomyocytes via the gp130/jak2/stat3 pathway41, is overexpressed in proliferating cardiomyocytes following ventricular amputation24. Moreover, impaired regulation of il6st promotes adverse remodeling following myocardial infarction (MI)42. The disintegrin and metalloproteinase adam8 mediates cell adhesion, migration, signaling and angiogenesis via proteolysis of various substrates43. Interestingly, adam8 improves muscle fiber regeneration by regulating inflammatory reactions that are necessary to eliminate injured fibers prior to the regeneration step44, and single nucleotide polymorphisms in this gene are associated with MI45. Among other top-ranked hubs, stx11a and cd63 are markers of intracellular vesicles, while arpc5a and cotl1 are actin-binding proteins that may modulate cell migration and immune response during heart regeneration46,47. Other top-ranked hub genes are regulators of the inflammatory response, such as: mrc1b, tmem154, igsf6 and cd22, which corroborates the importance of the immune response in heart regeneration (online resource). Lastly, other hubs (such as si:ch211-264f5.2) are still uncharacterized and represent interesting candidates for future investigation.

Hubs are relevant to heart regeneration in mammals

We investigated the biological importance of the network hubs in different mammals with limited, but inducible, heart regeneration capacity. First, we determined the levels of homology of the hubs in mouse, rat and human. We found that a large majority of hubs have orthologs in human (78% of hubs), mouse (79%) and rat (78%) (Supplementary Table S2). Furthermore, hubs are statistically enriched in evolutionary conserved genes in comparison to other (non-hub) genes in our network. We detected this significant association in humans (P = 0.02, 2 = 5.67, Chi-square test), mouse (P = 0.0001, 2 = 15.15) and rat (P = 0.0003, 2 = 13.18).

Among our hubs, there are genes with mouse homologous whose importance in neonatal heart regeneration have been previously reported. Some of such hubs (LOC100331505, csf2rb, max, rb1, epas1a) mapped to DEGs or their putative regulators in a recent heart regeneration study by O’Meara et al.48. Others (zgc:123190, zgc:77517, ctssb.1) have homologous genes that are DEGs in fully regenerated hearts in mice as reported by Haubner et al.49.

Next, we investigated whether the homologous genes are implicated in regulatory processes that are crucial for heart regeneration in mammals. We analyzed regulatory relationships between our hubs and 20 microRNAs (miRs) whose capacity to induce heart regeneration in mammals has been previously demonstrated19,50. Specifically, we addressed the question of whether our hubs are potential targets of miRs that are known to function as regeneration drivers. This was done by first generating a list of experimentally-validated miRs and their interactions in mouse, rat and human. Also we gathered putative miR-target interactions predicted by multiple computational techniques. We then searched for orthologs of our network hubs in the resulting datasets (Fig. 6). In humans, two of our hubs (fam49ba and il6st) are known targets of hsa-miR-590-3p, which has the capacity to trigger cardiomyocyte proliferation in (neonatal and adult) mice and rats50. We also identified 18 (unique) computationally-inferred interactions between our hubs and other heart regeneration miRs: hsa-miR-1, hsa-miR-195 and hsa-miR-199a51. In mouse and rat, we did not detect experimentally-validated interactions between regeneration miRs and hub orthologs, but we found hundreds of computationally-inferred hits. The latter included not only putative interactions with the mammalian heart regeneration miRs, but also with other miRs known to be relevant to cardiac cell proliferation and differentiation in mammals (Fig. 6).

Figure 6: Network hubs are functionally important for heart regeneration in mammals. It offers an overall view of biologically important associations between hubs and miRs known to be drivers of regeneration in mammals. Higher resolution views, including their integration with other biological information, are provided on the website. (A) miR-hub interactions in humans. (B) miR-hub interactions in the rat. (C) miR-hub interactions in the mouse. Lines are used to indicate miR-target interactions, and are colored to group miR-specific interactions. Full size image

A closer look at these interactions showed that, for instance, miR-199a targets the homologs of esrrga and rb1. Esrrg is highly expressed in the heart at fetal and postnatal stages, where it coordinates the oxidative metabolic program52. This gene is crucial for promoting the reprogramming of fibroblasts into cells of the cardiac lineage53. Rb1 plays a fundamental role in priming embryonic stem cells toward cardiac cells54. As in the case of miR-199a, miR-195 also targets esrrg. Moreover, mirR-195 targets alox5 and epas1a. The human homolog of alox5 is vital to improve healing after MI through the regulation of inflammation and collagen production55. Epas1 promotes angiogenesis and may support the adaptation of cardiomyocytes to hypoxia during heart failure56. Thus, our analyses show that mammalian homologs of network hubs are targeted by functionally-relevant miRs, which strengthens the biological importance of our predictions. Furthermore, our approach expands knowledge of miR-target associations that may be relevant to understand, and possibly elicit, heart regeneration in mammals.

A web resource enables mining of key cardiac regeneration genes in zebrafish

We integrated the time-course differential expression, module/hub analysis and microRNA data into a self-contained Web resource enabling exploration of our data through an intuitive interface. Figure 7 illustrates how the interface can be used to identify a potential key gene and associated correlation network involved in the regeneration process. Since each module represents a set of correlated genes modulated during the cardiac regeneration process, we suggest using the modules as a starting point for investigating gene signatures. In the example shown in Fig. 7, module 7A, containing 179 genes, is first selected (Fig. 7A-1). Through the interface we can see that this module, along with several other modules, is significantly enriched for genes involved in embryonic development and muscle cell differentiation (Fig. 7A-2). We can select to show only genes that are “Hub” nodes according to our analysis (Fig. 7A-3); this highlights 15 genes. Further selection based on differential expression at early time-points (Fig. 7B-1) identifies Il6st as a “hub” gene that is up-regulated early on in the cardiac regeneration process (from 4 hpi) (Fig. 7A-2). The network of genes most correlated with il6st can be explored and the expression fold changes of the genes viewed side-by-side against il6st (Fig. 7C). This reveals a number of highly correlated genes including Jak1, which is the tyrosine kinase responsible for transducing the signal from the Il6st receptor complex, and Stat3, a transcriptional co-activator of the signaling cascade. Inspection of potential mammalian microRNAs targeting il6st indicates that hsa-miR-590-3p (Fig. 7B-3), which has been previously shown to impact cardiac regeneration, potentially regulates this gene at the post-transcriptional level. This is just one example of how new hypotheses on pathways involved in the regeneration process can start to be built rapidly using our Web resource. A more detailed step-by-step example is available on the website (help section).