A metabolic switch from OXPHOS to glycolysis is essential for cardiomyocyte proliferation in the regenerating heart

The capacity to regenerate damaged tissues, such as the heart, various enormously amongst species. While heart regeneration is generally very low in mammals 1–3, it can regenerate efficiently in certain amphibian and fish species 4,5. Zebrafish has been used extensively to study heart regeneration, resulting in the identification of proliferating cardiomyocytes that drive this process 5–7. However, mechanisms that drive cardiomyocyte proliferation are largely unknown. Here, using a single-cell mRNA-sequencing approach, we find a transcriptionally distinct population of dedifferentiated and proliferating cardiomyocytes in regenerating zebrafish hearts. While adult cardiomyocytes are known to rely on mitochondrial oxidative phosphorylation (OXPHOS) for energy production, these proliferating cardiomyocytes show reduced mitochondrial gene expression and decreased OXPHOS activity. Strikingly, we find that genes encoding rate-limiting enzymes of the glycolysis pathway are induced in the proliferating cardiomyocytes, and inhibiting glycolysis impairs cardiomyocyte cell cycle reentry. Mechanistically, glycolytic gene expression is induced by Nrg1/Erbb2 signaling, and this is conserved in a mouse model of enhanced regeneration. Moreover, inhibiting glycolysis in murine cardiomyocytes abrogates the mitogenic effects of Nrg1/ErbB2 signaling. Together these results reveal a conserved mechanism in which cardiomyocytes undergo metabolic reprogramming by activating glycolysis, which is essential for cell cycle reentry and heart regeneration. This could ultimately help develop therapeutic interventions that promote the regenerative capacity of the mammalian heart.


Summary
The capacity to regenerate damaged tissues, such as the heart, various enormously amongst species. While heart regeneration is generally very low in mammals [1][2][3] , it can regenerate efficiently in certain amphibian and fish species 4,5 . Zebrafish has been used extensively to study heart regeneration, resulting in the identification of proliferating cardiomyocytes that drive this process [5][6][7] . However, mechanisms that drive cardiomyocyte proliferation are largely unknown. Here, using a single-cell mRNAsequencing approach, we find a transcriptionally distinct population of dedifferentiated and proliferating cardiomyocytes in regenerating zebrafish hearts. While adult cardiomyocytes are known to rely on mitochondrial oxidative phosphorylation (OXPHOS) for energy production, these proliferating cardiomyocytes show reduced mitochondrial gene expression and decreased OXPHOS activity. Strikingly, we find that genes encoding rate-limiting enzymes of the glycolysis pathway are induced in the proliferating cardiomyocytes, and inhibiting glycolysis impairs cardiomyocyte cell cycle reentry. Mechanistically, glycolytic gene expression is induced by Nrg1/Erbb2 signaling, and this is conserved in a mouse model of enhanced regeneration. Moreover, inhibiting glycolysis in murine cardiomyocytes abrogates the mitogenic effects of Nrg1/ErbB2 signaling. Together these results reveal a conserved mechanism in which cardiomyocytes undergo metabolic reprogramming by activating glycolysis, which is essential for cell cycle reentry and heart regeneration. This could ultimately help develop therapeutic interventions that promote the regenerative capacity of the mammalian heart.

Results
After injury, the zebrafish heart regenerates by proliferation of a small population of existing cardiomyocytes located mostly adjacent to the injury site, also known as the border zone [6][7][8] . These cycling cardiomyocytes exist within a heterogeneous cell population including non-proliferating cardiomyocytes, endothelial cells and immune cells, hampering the discovery of genetic programs specific for these proliferating cardiomyocytes [9][10][11] . Thus, a single cell approach is required to isolate them for detailed analysis. We generated a transgenic zebrafish nppa reporter line (TgBAC(nppa:mCitrine)) in which mCitrine expression recapitulates endogenous expression of the border zone stress gene nppa (which encodes a natriuretic peptide) (Extended Data Fig. S1). Histochemical analysis of injured adult hearts revealed that 75% (± 7%, n=3 hearts) of the nppa:mCitrine expressing cells were proliferating cardiomyocytes (Fig. 1a). To isolate these proliferating border zone cardiomyocytes, nppa:mCitrine hearts were cryo-injured, followed by cell dissociation and FACS sorting of mCitrine high cells (Fig. 1b). In addition, we FACS sorted mCitrine low cells to obtain non-proliferating remote cardiomyocytes from the same tissue for comparative analysis. Individual, living cells were sorted, followed by single-cell mRNAsequencing using the SORT-seq (SOrting and Robot-assisted Transcriptome SEQuencing) platform 12 (Supplementary Data 1). To identify the cardiomyocytes amongst the other cell types, we first identified the different cell types based on their transcriptomes. k-medoids clustering of the single cell transcriptomes by the RaceID clustering algorithm was used 13 (Fig. 1c and Supplementary Data 2), and visualized in two dimensions using t-distributed stochastic neighbor embedding (t-SNE) (Fig. 1d).
A total of 12 cell clusters were identified, including a large group of cardiomyocytes (clusters 1,2,4,7 and 9), a smaller group of endothelial cells (clusters 5,6,8,10 and 12), and some fibroblasts (cluster 3) and immune cells (cluster 11) using the expression of marker genes for specific cell types (Extended Data Fig. S2). Based on the transcriptome clustering, the cardiomyocytes fell into four main transcriptionallydefined clusters (1, 2, 4 and 7), indicating that the injured heart contained subgroups of cardiomyocytes. To address whether the border zone cardiomyocytes were enriched in one of the four cardiomyocyte clusters we compared the mCitrine fluorescence intensity (recorded during FACS sorting) of the cardiomyocyte and found that the average intensity was highest in cluster 7 and lowest in cluster 2 (Extended Data Fig.   S4). In addition we analysed the single-cell transcriptome data for the expression of known border zone genes (nppa, vmhc and mustn1b) and again found that cells expressing these genes were mostly in cluster 7 ( Fig. 1e and Extended Data Fig. 3a).
Together, these results indicate two things: first, the border zone cardiomyocytes (grouped in cluster 7) can be identified as a separate group in the single-cell RNA-seq data. Secondly, these border zone cardiomyocytes are transcriptionally distinct from remote cardiomyocytes (grouped in cluster 2), while two intermediate cardiomyocyte clusters lie in between.
Cardiomyocytes in the border zone disassemble sarcomeric structures and re-express markers of embryonic cardiomyocytes suggesting their dedifferentiation [6][7][8]14 . We therefore wanted to address the level of dedifferentiation of cluster 7 cardiomyocytes by comparing their transcriptome with embryonic cardiomyocytes. To obtain embryonic cardiomyocytes we performed FACSorting on embryos expressing the cardiomyocyte specific marker Tg(myl7:GFP). Single-cell mRNA-sequencing was performed and combined with the single-cell data from the injured adult hearts (Supplementary Data 3, 4 and 5). The RaceID algorithm identified several cell clusters (Extended Data Fig. S5) with separate clusters for the embryonic and adult cardiomyocytes (Fig. 2a). Importantly, the cluster 7 cardiomyocytes identified in the adult data analysis had a transcriptome that was highly similar to embryonic cardiomyocytes, as shown by pairwise correlation of the differentially expressed genes between the cardiomyocyte clusters: only 184 genes (p-value <0.01), out of 23,707 total detected genes, were differentially expressed between the embryonic and the cluster 7 (border zone) adult cardiomyocytes (Supplementary Data 8). In contrast, over 1800 genes (p-value <0.01) were differentially expressed between embryonic and cluster 2 (remote zone) adult cardiomyocytes.
To identify a possible dedifferentiation route we used part of the RaceID algorithm (StemID) that uses the single cell transcriptome data and cell clustering to derive a branched lineage tree 15 . The algorithm is based on the premise that stem cells and less differentiated cells tend to exhibit more uniform transcriptomes than differentiated cells, which express smaller numbers of genes at higher rates 16 . Using this approach, we found large differences in transcriptome entropy, resulting in low (cluster 2), intermediate (clusters 1 and 4) and high (cluster 7) StemID scores (Fig. 2b). This gradual increase suggests a dedifferentiation axis from cells in cluster 2 (remote myocardium) to cells in cluster 7 (border zone myocardium) and is in good agreement with our finding that the transcriptome of cluster 7 cardiomyocytes resembles an embryonic cardiomyocyte transcriptome (Fig. 2c). Together, these results indicate that clusters 4 and 7 are enriched for dedifferentiated and proliferative border zone cardiomyocytes while clusters 1 and 2 are enriched for differentiated remote cardiomyocytes.
While cardiomyocytes undergo a well-defined sequence of morphological and transcriptional changes during differentiation, very little is known about the reverse process. Ordering whole-transcriptome profiles of single cells with an unsupervised algorithm can resolve the temporal resolution during differentiation by identifying intermediate stages of differentiation without a priori knowledge of marker genes 17 . In this manner, the single-cell mRNA-seq experiment will constitute an in-silico time series, with each cell representing a distinct state of differentiation along a continuum.
To analyze the transcriptional changes occurring during this apparent dedifferentiation, the most likely dedifferentiation path, based on the StemID scores, was chosen starting at cluster 2 and progressing through clusters 1, 4 and 7. Next, gene expression profiles along this pseudo-temporal order were computed for all detected genes using the singlecell transcriptomes. These gene expression profiles were grouped into modules of coexpressed genes using self-organizing maps (SOMs), resulting in 14 modules (Fig. 2d, Supplementary Data 7). Corroborating our hypothesis of varying differentiation states, we observed that gene expression within these modules changed smoothly over pseudo time. We next analyzed the temporally-ordered expression profiles and identified four groups of genes that shared the same dynamics of expression during this differentiation trajectory. The first group (modules 1, 2) contained genes that were most highly Adult cardiomyocytes rely on fatty acid metabolism and aerobic mitochondrial OXPHOS rather than glycolysis for their energy production 18 . It was therefore surprising to see that the pseudotime analysis shows a concomitant down regulation of genes transcribed from mitochondrial DNA and upregulation of glycolysis genes. To address mitochondrial activity directly, we measured succinate dehydrogenase (SDH) enzyme activity, located in the inner mitochondrial membrane that functions in both the citric acid cycle and electron transport chain. We observed a 40% reduction in SDH activity specifically in the border zone cardiomyocytes as compared to the remote cardiomyocytes ( Fig. 3a). In agreement with the reduced mitochondrial OXPHOS activity, transmission electron microscopy (TEM) imaging revealed more immature mitochondria in border zone cardiomyocytes evidenced by their altered morphology and reduced cristae density (Fig. 3b), which is consistent with previous reports linking mitochondrial function with morphology 19,20 . Since the pseudotime analysis suggested Together these observations indicate that, upon cardiac injury, border zone cardiomyocytes switch their metabolism from aerobic mitochondrial OXPHOS to glycolysis, which is necessary for cell cycle reentry.
Next, we investigated the upstream signals that drive the observed metabolic switch in border zone cardiomyocytes during cell cycle reentry. Injury-induced Neuregulin 1 (Nrg1) expression is a potent mitogen that induces cardiomyocyte dedifferentiation and cell cycle reentry by activating ErbB2 receptor signaling 21 . Furthermore, Nrg1 can induce glucose metabolism in skeletal muscle cells and cardiomyocytes in vitro 22,23 .
Therefore, we addressed whether Nrg1 administration can induce glycolysis in vivo. To test this we used a previously described transgenic zebrafish model in which Nrg1 overexpression (OE) can be induced in a heart specific manner, tg(cmlc2:CreER; β-act2:BSNrg1) 21 . After tamoxifen injection into tg(cmlc2:CreER; β-act2:BSNrg1) fish, we observed a profound and consistent upregulation of glycolysis genes in the myocardium (Fig. 4a). Expression of hk1, encoding the rate limiting glycolytic enzyme hexokinase, was strongly induced in the ventricular wall of Nrg1 OE hearts, which correlates well with the observed induction of cardiomyocyte dedifferentiation and proliferation in this region 21 . From these results, we conclude that glycolysis gene expression can be induced by Nrg1 signaling even in the absence of cardiac injury.
While the mammalian neonatal heart regenerates very efficiently, this regenerative capacity is lost soon after birth 24,25 . The loss of regenerative capacity in the first week after birth correlates well with the observed change in cardiomyocyte metabolism from glycolysis to mitochondrial OXPHOS 18 . This has been attributed to an increase in oxygen availability resulting in accumulating reactive oxygen species (ROS) and subsequently DNA damage and reduced proliferation in adult cardiomyocytes [26][27][28] . A direct role for glycolysis and its role on cardiomyocyte proliferation and mammalian heart regeneration has not been addressed. To test the hypothesis that glycolysis might be beneficial to cardiomyocyte proliferation and regeneration, we first determined whether a metabolic switch from OXPHOS to glycolysis occurs in a murine model in which regeneration can be induced, such as the transgenic constitutively active ErbB2 receptor overexpression system (caErbB2 OE) 29 . Consistent with our model, we observed that critical glycolysis genes (e.g. Pfkp, Pdk3 and Pkm2) including glucose and lactate transporters (Slc16A3 and Slc2A1) were significantly upregulated in caErbB2 OE cardiomyocytes while genes transcribed from mitochondrial DNA were downregulated (Fig. 4b). Pdk3 encodes a pyruvate dehydrogenase kinase, which phosphorylates pyruvate dehydrogenase (PDH). PDH is a mitochondrial multi-enzyme complex that converts pyruvate to Acetyl-CoA and provides a primary link between glycolysis and the TCA cycle. Upon phosphorylation by PDK, p-PDH is inactivated and pyruvate is diverted away from the TCA cycle resulting in enhanced lactate production. Consistent with the increase in pdk3 expression, an increased phosphorylation of PDH was observed in the myocardium of caErbB2 OE hearts (Fig.   4c).
Secondly, we addressed whether the observed switch in metabolic gene expression correlated with enhanced regenerative capacity after injury. Patients and animal models with myocardial infarction (MI) show reduced fatty acid metabolism and increased glucose uptake and glycolysis 30,31 , which is explained by the reduced oxygen availability in the ischemic region since glycolysis is oxygen independent 32 33 .
Therefore, we performed MI in wild type and caErbB2 OE hearts. Even though glucose uptake and glycolytic enzyme activity in the ischemic area are increased by MI 30 31 , we observed a stronger upregulation of glycolytic gene expression and decreased mitochondrial gene expression in caErbB2 OE hearts with MI compared to wild type hearts with MI (Extended Data Fig. S8). This stronger and consistent upregulation of glycolytic gene expression in caErbB2 OE hearts correlates with the reported enhanced cardiomyocyte proliferation and improved regeneration 29 . These findings imply that the enhanced glucose uptake and glycolysis observed after MI might be sufficient for cell survival, but that a stronger induction of glycolysis gene expression is required to fully stimulate cardiomyocyte proliferation and regeneration.
Thirdly, we addressed whether the observed metabolic switch to glycolysis in caErbB2 OE cardiomyocytes is required for their reentry into the cell cycle. Corroborating our model, we indeed observed that treating caErbB2 OE cardiomyocytes in vitro with the glycolysis inhibitors 2-DG or lonidamine (inhibiting hexokinase activity) strongly and significantly impaired cell cycle reentry (Fig 4e) and cytokinesis (Extended Data Fig.   9). Together these results demonstrate that in the murine cardiomyocytes, ErbB2 signaling drives a metabolic switch from OXPHOS to glycolysis, which is required for their cell cycle reentry. These results suggest that this metabolic switch in cardiomyocytes is beneficial for heart regeneration.
In conclusion, the use of single cell transcriptomics has allowed us to identify and characterize the proliferating cardiomyocyte population that arises following injury in the regenerating zebrafish heart. Most striking is that these proliferative cardiomyocytes exhibits a metabolic reprogramming hallmark, whereby genes encoding regulators for glucose transport and glycolysis are strongly induced, while simultaneously reducing mitochondrial activity and undergoing dedifferentiation. This metabolic reprogramming phenomenon is necessary for cell cycle reentry of adult cardiomyocytes, and Nrg1/ErbB2 signaling mechanistically acts upstream of this event.
More importantly, metabolic reprogramming appears to be a conserved mechanism for cardiomyocyte cell cycle reentry since mammalian models of adult cardiomyocyte proliferation also display this reprogramming hallmark. Likewise, alteration of metabolic states directly influenced proliferative capacity, further demonstrating its crucial role for mammalian cardiomyocyte cell cycle reentry.
This raises the question why cardiomyocytes need to switch their metabolism to enter the cell cycle? Interestingly, a similar metabolic shift from aerobic OXPHOS to glycolysis occurs in proliferating tumor cells, known as the Warburg effect 34 . Furthermore, progenitor cells in the developing embryo as well as induced pluripotent stem cells, also depend on glycolysis to maintain proliferation and their potency 35,36 . Glycolysis provides essential metabolites that are needed for cellular reprogramming and proliferation 37 , but glycolytic enzymes can also directly interact with cell cycle regulators to promote proliferation 38 39 . Our work suggests that therapeutic targeting of metabolic reprogramming in cardiomyocytes may be beneficial to repair damaged hearts.

Transgenic zebrafish lines and cryoinjury
All procedures involving animals were approved by the local animal experiments committees and performed in compliance with animal welfare laws, guidelines and policies, according to national and European law.
The following fish lines were used: myl7:GFP twu34Tg 40 . The Tg(cmlc2:CreER; β-act2:BSNrg1) line was used as described before 21 . The TgBAC(nppa:mCitrine) line was generated essentially as described previously 41  EMBRYONIC: Live embryos were immobilized using ms222 and embedded in nitrocellulose + E3 to be mounted on a Leica SPE confocal microscope, followed by a Z-stack maximum projection (step size 2 µm).

MAMMALIAN P7 CARDIAC CULTURES: For immunofluorescence, cardiac
cultures were fixed with 4% PFA for 10 minutes on room temperature on the shaker, followed by permeabilization with 0.5% Triton X-100 in PBS for 5min, and blocking with 5% bovine serum albumin (BSA) in PBS containing 0.1% Triton for 1h at room temperature.

Quantitative PCR
RNA from whole hearts was isolated using the MiRNeasy kit (Qiagen, 217004), according to the manufacturer's instructions. RNA was quantified using a NanoDrop

In situ hybridization
PARAFFIN: After o/n fixation in 4% PFA, hearts were washed in PBS twice, dehydrated in EtOH, and embedded in paraffin. Serial sections were made at 10 µm. In situ hybridization was performed on paraffin-sections as previously described 45 except that the hybridization buffer used did not contain heparin and yeast total RNA. CRYOSECTIONS: Sections were obtained as described earlier. In situ hybridization was performed as for paraffin, however sections were pre-fixed for 10 minutes in 4% PFA + 0.25% glutaraldehyde before Proteinase K treatment. Moreover, slides were fixed for 1 hour in 4% PFA directly after staining.

Isolation of single cells from cryoinjured hearts
Cryoinjured hearts (n=13) were extracted at 7 dpi. Cells were dissociated according to 46 . For cell sorting, viable cells were gated by negative DAPI staining and positive YFPfluorescence. In brief, the FACS gating was adjusted to sort cells for nppa:mCitrine high (to enrich for proliferating cardiomyocytes) and nppa:mCitrine low (remote cardiomyocytes and other cell types) cells. In total n=576 mCitrine high and n=192 mCitrine low cells were sorted into 384-well plates and processed for mRNA sequencing as described below.

Single-cell mRNA sequencing
Single-cell sequencing libraries were prepared using SORT-seq 12

Bioinformatic analysis
To analyze the single-cell RNA-seq data we used the previously published RaceID The StemID algorithm were used as previously published 15  StemID score, which is representative to "stemness" of a cell population.

Inference of co-expressed gene modules
To identify modules of co-expressed genes along a specific differentiation trajectory (defined as a succession of significant links between clusters as identified by StemID) all cells assigned to these links were assembled in pseudo-temporal order based on their projection coordinate. Next, all genes that are not present with at least two transcripts in at least a single cell are discarded from the sub-sequent analysis. Subsequently, a local regression of the z-transformed expression profile for each gene is computed along the differentiation trajectory. These pseudo-temporal gene expression profiles are topologically ordered by computing a one-dimensional self-organizing map (SOM) with 1,000 nodes. Due to the large number of nodes relative to the number of clustered profiles, similar profiles are assigned to the same node. Only nodes with more than 3 assigned profiles are retained for visualization of co-expressed gene modules.
Neighboring nodes with average profiles exhibiting a Pearson's correlation coefficient >0.9 are merged to common gene expression modules. These modules are depicted in the final map. Analyses were performed as previously published 15 .

Transmission Electron Microscopy
Hearts were excised and immediately chemically fixated at room temperature with

Quantification of mitochondrial parameters.
In every heart, each in the borderzone and remote myocardial region, 100 welldeliniated mitochondria with clearly visible outer and inner membranes were selected.
Mitochondrial perimeter and surface were measured using the freehand tool of Image J. The perimeter to surface ratio was calculated and used as a factor that describes the pluriformity of mitochondria. The amount of cristae was estimated by counting the in Gelatin-coated (0.1%, G1393, Sigma) wells with DMEM/F12 medium supplemented with L-glutamine, Na-pyruvate, nonessential amino acids, penicillin, streptomycin, 5% horse serum and 10% FBS ('complete-medium') at 37•C and 5% CO2 for 24h.
Afterwards, medium was replaced with FBS-depleted medium (otherwise same composition) for additional 48hours of culture in either 3mM 2DG (Sigma-Aldrich) or 80µM lonidamine before further processing.

Ex vivo glucose uptake
Fish were euthanized on ice water before hearts were extracted in PBS + heparin and were allowed to bleed out for 15 minutes. Hearts were then transferred into fresh PBS + 10%KCl to stop the heart from beating and mounted directly on a glass bottom cell culture dish in 1% agarose. Thereafter, 2NBDG (Caymanchem #11046, 400µM) was added to the dish and the hearts were taken directly for imaging. Imaging was performed using a Leica SP5 multiphoton microscopy using 930nm laser excitation wavelength. 150µm z-stacks were made with z-step size 5µm every 5 minutes for 2 hours.

Gene set enrichment analysis (GSEA)
GSEA (Genepattern, Broad Institute) was performed to assess enrichment for glycolytic genes upregulated between cluster 7 and 2. A list of genes involved in zebrafish glycolysis was obtained from KEGG. As number of permutations 1000 was used, which means p=0 indicates p<0.001.    Time-lapse multi-photon confocal images of whole heart. The heart was isolated 7 days post injury and incubated with 2-NBDG, a fluorescent glucose analogue, at t=0. Confocal images were taken every 5 min to monitor cellular glucose uptake. Dotted line indicates injury area. Arrows point to regions of the border zone where 2-NBDG uptake is enhanced compared to the remote area. Scale bar represents 100 µm. Figure 8: Activation of ErbB2 signaling in murine MI model induces glycolytic gene expression while repressing mitochondrial genes. (a) Cartoon showing experimental procedure. Myocardial infarction (MI) was induced by left anterior descending coronary artery ligation at the same day when Dox was removed to induce cardiomyocyte specific caErbB2 overexpression (OE). Whole hearts where isolated for mRNA extraction 14 days after the MI and caErbB2 induction. (b) qPCR results with relative mRNA fold changes comparing wild type (WT, blue bars) with caErbB2 OE (OE, orange bars). Note the significant upregulation of glycolytic genes (Pfkp, Pdk3 and Pgk1) and significant downregulation of genes transcribed from the mitochondrial DNA (Mt-Nd1, Mt-Co2, Mt-Co1 and Mt-Co3). * p<0,05; ** p<0,01; *** p<0,001 Supplemental Figure 9: Glycolysis is required for caErbB2 induced cell division. (a) Cartoon showing the experimental procedure to analyse the effects of glycolysis inhibitors (2-DG and lonidamine) on cardiomyocyte proliferation. (e) Confocal images of cells from control and caErbB2 OE hearts stained for the cytokinesis marker AuroraB kinase (green) and myosin heavy chain (aMHC, red) to identify cardiomyocytes and DAPI (Blue). Cells were either control treated or treated with the glycolysis inhibitors Lonidamine or 2-DG. Arrows point to dividing cardiomyocytes. Scale bar represents 50 µm. ** p<0,01