Molecular Atlas of Postnatal Mouse Heart Development

Background The molecular mechanisms mediating postnatal loss of cardiac regeneration in mammals are not fully understood. We aimed to provide an integrated resource of mRNA, protein, and metabolite changes in the neonatal heart for identification of metabolism‐related mechanisms associated with cardiac regeneration. Methods and Results Mouse ventricular tissue samples taken on postnatal day 1 (P01), P04, P09, and P23 were analyzed with RNA sequencing and global proteomics and metabolomics. Gene ontology analysis, KEGG pathway analysis, and fuzzy c‐means clustering were used to identify up‐ or downregulated biological processes and metabolic pathways on all 3 levels, and Ingenuity pathway analysis (Qiagen) was used to identify upstream regulators. Differential expression was observed for 8547 mRNAs and for 1199 of 2285 quantified proteins. Furthermore, 151 metabolites with significant changes were identified. Differentially regulated metabolic pathways include branched chain amino acid degradation (upregulated at P23), fatty acid metabolism (upregulated at P04 and P09; downregulated at P23) as well as the HMGCS (HMG‐CoA [hydroxymethylglutaryl‐coenzyme A] synthase)–mediated mevalonate pathway and ketogenesis (transiently activated). Pharmacological inhibition of HMGCS in primary neonatal cardiomyocytes reduced the percentage of BrdU‐positive cardiomyocytes, providing evidence that the mevalonate and ketogenesis routes may participate in regulating the cardiomyocyte cell cycle. Conclusions This study is the first systems‐level resource combining data from genomewide transcriptomics with global quantitative proteomics and untargeted metabolomics analyses in the mouse heart throughout the early postnatal period. These integrated data of molecular changes associated with the loss of cardiac regeneration may open up new possibilities for the development of regenerative therapies.

Remarkably, certain amphibians and fish as well as neonatal rodents can fully regenerate their hearts after an injury. [7][8][9] This occurs predominantly through proliferation of remaining cardiomyocytes in the areas adjacent to the injury. 8,10 However, rodents lose this regenerative capacity within 7 days after birth because of cardiomyocyte cell cycle withdrawal, 8 which is considered to represent a major hurdle for the development of regeneration-inducing therapies. [11][12][13] In line with cardiac regeneration in neonatal rodents, full functional recovery of a newborn baby with a massive MI suggests that humans possess a similar intrinsic capacity to regenerate their hearts at birth. 14 Consequently, it is crucial to elucidate the molecular mechanisms mediating the postnatal loss of cardiac regeneration.
At birth, the resistance of pulmonary circulation drops dramatically because of the expansion of lung alveoli, causing an increase in systemic blood pressure and shunt closure. This in turn induces a steep increase in arterial blood oxygen content. Increased cardiac workload and altered metabolic environment induce a switch from anaerobic glycolysis, which is the main source of energy in the embryonic heart, to mitochondrial fatty acid b-oxidation soon after birth. 15 This switch is controlled by signaling that involves HIF1a (hypoxia-inducible factor 1a) and HAND1 (heart and neural crest derivatives expressed 1). 16 At birth, HAND1 expression is rapidly downregulated, allowing cardiomyocytes to shut down glycolysis and initiate lipid oxidation. Failure to downregulate HAND1 and thereby initiate lipid oxidation is fatal, emphasizing that oxidative metabolism is a prerequisite for sufficient energy production in the postnatal heart. Moreover, PPAR (peroxisome proliferator-activated receptor) signaling plays an important role in activating lipid metabolism. 17 The profound changes in the energy metabolism of cardiomyocytes are associated with alterations in mitochondria: the fetal-type mitochondria undergo mitophagy and are replaced with mature adult-type mitochondria in a Parkindependent process to allow more efficient ATP production. 18 Increased oxidative metabolism promotes reactive oxygen species production and thereby induces a DNA damage response, which is believed to contribute to cardiomyocyte cell cycle arrest. 19 On a molecular level, endogenous mechanisms known to regulate cardiomyocyte proliferation include neuregulin-ERBB2 (erb-b2 receptor tyrosine kinase 2) signaling, the Hippo-YAP (YY1 associated protein 1) pathway, and the transcription factor GATA4 (GATA binding protein 4), among others. 11,[20][21][22] Moreover, even though metabolic pathways are controlled by the same signals that regulate cell proliferation, 23 their significance in cardiac regeneration and cardiomyocyte proliferation is unknown.
In search of regenerative therapies, a number of studies have utilized transcriptomics to identify mechanisms regulating cardiac regeneration (for recent examples, see Natarajan et al, 24 Quaife-Ryan et al, 25 and Kang et al 26 ). However, proteomic and metabolomic changes have not been thoroughly investigated. In this article, we report the first integrated study combining genomewide RNA sequencing (RNAseq), global proteomics, and untargeted metabolomics to characterize in detail the metabolic changes occurring throughout the early postnatal heart development.

Methods
The data, analytic methods, and study materials have been made available to other researchers for purposes of reproducing the results or replicating the procedure. Detailed methods are available in Data S1, and detailed protocols are available from the corresponding authors on reasonable request. The RNAseq data have been made publicly available at the National Center for Biotechnology Information gene expression and hybridization array data repository GEO. 27 The proteomics raw data have been made available at the MassIVE repository. 28 All other omics data are available in Appendices S1 through S5. The images and raw quantification data for cell viability and proliferation experiments are available from the corresponding authors on reasonable request.
The experimental design is summarized in Figure 1. Two separate sets of mouse (strain C57BL/6JOlaHsd, both sexes) ventricular tissue samples were used. Animal handling and all procedures were carried out in accordance with University of Helsinki institutional guidelines. Set 1 encompassed samples from postnatal day 1 (P01), P04, and P09, representing hearts with full regenerative capacity (P01), partial regenerative capacity (P04), and negligible regenerative capacity (P09), and was used for proteomics and metabolomics. To validate the results, a second set of samples was analyzed. In set 2, an additional sample group from 23-day-old mice (P23) was included to discriminate phenomena related to heart growth. The set 2 samples were subjected to transcriptomics using pooled samples (3 hearts per sample) and proteomics and metabolomics analyses (without pooling). The results from both sample sets are presented for proteomics and metabolomics, as the combination of both sample sets in metabolomics analyses using the linear mixed effects model produced more reliable results than either of the 2 sample sets alone. The methods and bioinformatics analyses are shown in Figure 1B. For transcriptomics analyses, statistically significant differential expression was defined as fold change >1.5 and q<0.01. For proteomics, metabolomics and bioinformatics analyses, q<0.01 with no fold change limit was considered statistically significant.

Transcriptomics
To analyze postnatal gene expression changes, RNAseq was carried out with ventricular tissue samples using 3 pooled samples from 9 hearts at each time point. Principal component analysis of RNAseq data showed clear grouping of samples and separation of sample groups ( Figure S1). Figure 2A presents hierarchical clustering and a heat map of 1000 genes with most significant changes between P01 and P04. In total, 8547 individual protein-coding genes were up-or downregulated statistically significantly (q<0.01, fold change >1.5); their numbers for each time-point comparison are presented in Figure 2B. The greatest changes in gene expression were observed between time-point comparisons P01 to P23 and P04 to P23. The top 10 up-and downregulated genes are presented in Figure S2.
The expression levels of cardiomyocyte-specific structural proteins exhibited an anticipated pattern with a switch from Myh7 to Myh6 (myosin heavy chain 6 and 7, respectively) and The postnatal loss of cardiac regenerative capacity is illustrated for comparison, and numbers of animals in each sample group are presented in the table. *Total sample sizes are indicated for metabolomics (5) and proteomics (14). B, Analysis techniques and bioinformatics analyses used in the study. GCxGC-MS, 2dimensional gas chromatography-mass spectrometry; GO, gene ontology; LC-MS, liquid chromatographymass spectrometry; LC-MS/MS, liquid chromatography-tandem mass spectrometry; RNAseq, RNA sequencing. Figure 2. Gene expression changes in the neonatal mouse heart. A, Heat map of the top 1000 genes with the smallest q values between postnatal day 1 (P01) and P04. B, The numbers of up-and downregulated genes (q<0.01, fold change >1.5). C, Expression patterns of selected cardiomyocyte-specific structural proteins. The data were normalized to P01 and are expressed as meanAESEM (n=3 pooled samples, each from 3 hearts). D, Selected significantly enriched (q<0.01) biological process gene ontology (GO) terms for each time point comparison from gene set enrichment analysis. Blue indicates downregulation, and red indicates upregulation. All gene symbol explanations are available in Appendix S1. from Tnni1 to Tnni3 (troponin I1 and troponin I3, respectively) within the early postnatal period ( Figure 2C). 29,30 In addition, the expression profiles of the main cardiomyocyte ion channels ( Figure S3A) were in line with previous reports. 31 The expression patterns of control genes are presented in Figure S3B. Expression of Actb (actin beta), Rpl4 (ribosomal protein L4), Rpl32 (ribosomal protein L32), Tbp (TATA-box binding protein), Oaz1 (ornithine decarboxylase antizyme 1), and Pgk1 (phosphoglycerate kinase 1) did not change significantly, whereas several other genes either generally used or recommended as control genes, 32 such as Gapdh (glyceraldehyde-3-phosphate dehydrogenase), were up-or downregulated (q<0.01 and fold change >1.5) in at least 1 time-point comparison. Furthermore, genes known to play a role in cardiomyocyte proliferation (Erbb2, Gata4), metabolic switch from glycolysis to fatty acid oxidation (Hand1, Hif1a, Ppar isoforms), and mitochondrial maturation (Park2, parkin RBR E3 ubiquitin protein ligase; Pink1, PTEN induced kinase 1) exhibited expected changes (Table S1). Differential expression analysis of all protein-coding genes is available in Appendix S1.
To investigate biological processes linked to the observed gene expression changes, we carried out gene ontology (GO) 33,34 enrichment analysis, and changes with q<0.01 were considered significant. Changes in biological processes linked to cell proliferation, cardiac muscle, cell adhesion and motility, ion transport, and development of immune response were highlighted among the enriched GO terms ( Figure 2D and Appendix S2). In accordance with increased oxidative metabolism, genes linked to oxidative stress were upregulated at P23 compared with P04 or P09. DNA damage response was upregulated from P01 to P04, whereafter it was downregulated in all other comparisons. Of the GO terms associated with cellular metabolism, oxidation-reduction processes were upregulated at P23 compared with P01, P04, and P09; carbohydrate metabolism was upregulated from P04 to P23; and pyruvate metabolism was upregulated from P01 to P23 and P09 to P23. Sterol biosynthesis was downregulated in all other time-point comparisons except P01 to P04 and P09 to P23, indicating that the transcript-level changes take place mainly between P04 and P09. Amino acid metabolism and branched chain amino acid (BCAA) catabolism were upregulated from P09 to P23.

Proteomics
To quantify protein abundances in ventricular tissue, we used a shotgun proteomics approach. On average, 1140 and 1337 distinct proteins were quantified in individual samples in sample sets 1 and 2, respectively ( Figure 3A). The numbers and fold changes of differentially expressed proteins are shown in Figure S4A (set 1) and Figure 3B (set 2). The corresponding lists of differentially expressed proteins, as well as all quantified proteins and their label-free quantification intensities, are presented in Appendix S3. Hierarchical clustering of the data shows that the sample groups were separated from each other in both sample sets ( Figure 3C and Figure S4B); this was also seen in principal component analysis ( Figure S5). GO enrichment analysis and KEGG 35 pathway analysis were used to identify up-and downregulated processes and activated or inactivated pathways (Appendix S3). Oxidative metabolism and lipid metabolism were highlighted as upregulated phenomena at P04, P09, and P23 compared with P01, whereas in the downregulated GOs, nucleic acid metabolism and glycolysis were highlighted from P01 to P04, peptide metabolism and protein synthesis were highlighted from P01 to P09, and RNA processing and peptide metabolism were highlighted from P01 to P23.
The postnatal increased activity of energy metabolism and the shift from carbohydrates to fatty acids as the primary source of ATP were well represented in the proteomics data. Several components (ENO1 (enolase 1), PFKL (phosphofructokinase, liver type), PGK1 (phosphoglycerate kinase 1)) of the HIF1 signaling pathway that promotes anaerobic metabolism in the embryonic heart were downregulated after birth (Table S2). Furthermore, many key proteins related to fatty acid metabolism, degradation, and oxidative phosphorylation were upregulated after P01. This was paralleled with upregulation of most of the detected components of PPAR signaling, which promotes fatty acid oxidation, after P01. In line with the increasing mitochondrial content and maturation in cardiomyocytes, most mitochondrial proteins were highly upregulated from P01 to P23. These included not only enzymes linked to oxidative metabolism but also structural (IMMT (inner membrane mitochondrial protein)) and regulatory (PHB (prohibitin), DNAJA3 (DnaJ heat shock protein family (Hsp40) member A3)) proteins (Table S2). Furthermore, the mitochondrial isoforms of SOD2 (superoxide dismutase) and peroxiredoxins 3 and 5 (PRDX3, PRDX5) were strongly upregulated at P09 and P23 (Table S2), reflecting a response to oxidative stress.
Analysis of key proteins of cell cycle regulation and signaling pathway activation with proteomics is challenging because of the low abundance of these proteins. Nevertheless, expression of proteins related to microtubule dynamics during mitosis (DYNC1H1 (dynein cytoplasmic 1 heavy chain 1), PAFAH1B1 (platelet activating factor acetylhydrolase 1b regulatory subunit 1), RHOA (ras homolog family member A), and BUB3 (BUB3, mitotic checkpoint protein)) decreased with increasing postnatal age (Table S2). Most of the detected 14-3-3 proteins (YWHAB, YWHAE, YWHAH, YWHAQ, and YWHAZ (tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation proteins beta, epsilon, eta, theta and zeta)), some of which regulate cell cycle and cardiomyocyte proliferation, 36  to P23 in set 2 samples (Table S2). Components of the canonical Wnt pathway were also detected (Table S2), and the decreased abundance of b-catenin (CTNNB1) and its nuclear interaction partner Pontin52 (RUVBL) 38 indicate attenuation of canonical Wnt pathway activity from P01 to P09.

Metabolomics
For the untargeted metabolomics analysis, we used 2 complementary techniques to increase the metabolite coverage, namely, liquid chromatography-mass spectrometry (LC-MS) and 2-dimensional gas chromatography coupled to MS (GC9GC-MS), and combined the data from the 2 sample sets using the linear mixed effects model. Quality control results are available in Data S1. In total, 805 and 162 metabolic features changed statistically significantly (q<0.01) in at least 1 of the group comparisons (P01?P04, P01?P09, P01?P23) for analyses with LC-MS and GC9GC-MS, respectively (Appendix S4). Based on the MS/MS analysis of metabolic features, 151 significantly up-or downregulated metabolites were identified and are shown in Figure 4.
In line with the transcriptomics and proteomics results, the metabolomics data show a distinct transition from carbohydrates to lipids as the main source of energy. The abundances of glucose and sugar derivatives (glucose, glucose-6-phosphate, fructose-6-phosphate, galactose-6-phosphate) were downregulated at P23 compared with P01 ( Figure 4B), whereas the abundances of most fatty acids and components of glycerolipid metabolism (glycerol-3-phosphate, glycerol-2phosphate, glycerol, glyceric acid) increased after P01 ( Figure 4A and 4B). The increased fatty acid b-oxidation was also reflected as increased abundance of acylcarnitines at P04 and P09 ( Figure 4A). However, most metabolites of the citric acid cycle (citric acid, fumaric acid, pyruvic acid, malic acid, lactic acid), pentose-phosphate pathway (D-seduloheptulose-7-phosphate), and glycolysis (lactic acid, pyruvic acid) exhibited a constant rise with increasing postnatal age ( Figure 4B), reflecting a total increase in energy metabolism. No changes were detected in the levels of citric acid cycle metabolites succinic acid and a-ketoglutarate (Appendix S4). These data indicate a transition from carbohydrates to fatty acids as the main source of energy and a total increase in energy metabolism over the early postnatal period.
In addition to the general increase in fatty acid abundance, various lipid species (phospholipids, lysophospholipids, monoacylglycerides, and fatty acids) displayed interesting changes in abundance ( Figure 4A). For phospholipids, there was a general trend toward increased saturation level along with increasing postnatal age. All phospholipids showing a constant decrease from P01 to P23 were unsaturated, and many of them were polyunsaturated. The phospholipid species with saturated medium-to-long-chain fatty acids exhibited an increase at P04 and P09 followed by a decrease at P23. Several lysolipid species (LysoPC and LysoPE) increased at P04 and remained constant or increased further through P09 to P23, with the exception of lysolipids with fatty acids 16:1 and 18:1, which decreased after P01, or 14:0, which decreased after P04. Interestingly, myristic acid (fatty acid 14:0) exhibited a comparable pattern both as a free fatty acid and when incorporated into other lipid species (eg, ceramide, sphingomyelin, or phosphocholine): The abundance increased at P04 and P09 and decreased at P23. A similar pattern was also observed for other medium-chain saturated fatty acid species ( Figure 4A).
The levels of most amino acids displayed an initial increase at P04 and/or P09 followed by a decrease at P23 ( Figure 4B). In contrast to other (proteinogenic) amino acids, the abundance of glutamic acid, alanine and histidine increased initially but remained significantly higher at P23 compared with P01 ( Figure 4B). The increase in amino acid abundance from P01 to P09 likely reflects the active protein synthesis required during cardiomyocyte growth and maturation.
We also observed significant changes in purine metabolism, particularly among the metabolites of AMP catabolism. The levels of AMP had decreased already at P04, paralleled by an increase in the abundance of its degradation pathway metabolites inosine, hypoxanthine, and xanthine ( Figure 4B). This likely reflects the alterations in the energy metabolism, as high AMP abundance immediately after birth would promote fatty acid boxidation through AMPK (AMP-activated protein kinase). 39 In line with the postnatal increase in cardiac workload, we also observed increased abundance of creatine and creatinine at P09 and P23 compared with P01 ( Figure 4B). Other interesting metabolite findings include the increase of the ascorbic acid metabolites L-threonic acid and threonic acid 1,4-lactone on P09 and P23, respectively, reflecting increased oxidative stress with increasing postnatal age.

Multiomics Integration
For comprehensive multiomics integration, we utilized fuzzy cmeans clustering, in which the transcripts, proteins, and metabolites were assigned into ≥1 clusters based on their abundance patterns ( Figure 5A, Figure S6). We compared transcriptomics and proteomics clusters on transcript/protein levels but also for enriched GO terms and KEGG pathways (q<0.01). The percentage of proteomics clusters covered by the RNAseq clusters are shown in Figure 5B for the proteins/ transcripts and biological process GO terms and in Figure S7A for cellular component and molecular function GO terms and KEGG pathways. Selected enriched GO terms and KEGG pathways and their enrichment in each transcriptomics and proteomics cluster are presented in Figure 5C (cellular component GO terms in Figure S7B). Because fuzzy clustering could only be carried out separately for the 2 sample sets, it provided only a little more insight to the metabolomics data compared with the linear mixed effects model. Consequently, instead of analyzing the individual metabolomics clusters, we performed metabolite set enrichment analysis on all significantly changed and identified metabolites to investigate the metabolic pathways with which the significantly up-or downregulated metabolites were associated. Based on metabolite set enrichment analysis, the 3 most significantly enriched metabolic pathways were protein biosynthesis, urea cycle, and glycerolipid metabolism ( Figure 5D). Enrichment of multiple individual amino acid metabolic pathways was also observed. The high enrichment of urea cycle metabolites correlated with the changes in individual metabolites, such as arginine and fumarate ( Figure 4B). Not all identified urea cycle metabolites, however, exhibited statistically significant changes (eg, aspartate and urea; Appendix S4).
Multiomics integration for the KEGG pathway "glycolysis and gluconeogenesis" is presented in Figure S8 as an example. In line with the data from the individual omics analyses, glycolysis was not uniformly activated or inactivated. Instead, the genes and proteins in the beginning and the end of the pathway were either upregulated at later time points or exhibited variable abundance patterns, whereas the genes and proteins in the middle of the pathway were either downregulated with increasing postnatal age or displayed variable expression patterns. Most enzymes of the related pyruvate metabolism pathway were also differentially expressed on mRNA and/or protein levels within the early postnatal period ( Figure S9).
To identify potential upstream regulators of gene and protein expression changes in each RNAseq and proteomics cluster, we subjected the individual clusters to Ingenuity pathway analysis (Qiagen). 40 Upstream regulators belonging to the classes transcription regulator, microRNA, and chemical-endogenous mammalian were included in the analyses. Of transcriptional regulators, HAND1 was identified as a potential upstream regulator for RNAseq cluster 5 and proteomics cluster 3 (Appendix S5), which is in line with its known role in the regulation of postnatal energy metabolism. In the microRNA analyses, miR-1, miR-21, miR-122, and let-7 were identified as potential upstream regulators of various mRNA and protein clusters (Appendix S5). Furthermore, the analysis of endogenous chemicals identified several interesting metabolites as potential regulators of up-and downregulated mRNAs and proteins-palmitic acid, cholesterol, fatty acid, amino acids, and butyric acid (Appendix S5)-correlating well with the metabolomics data. This multiomics integration and upstream regulator analysis prompted us to investigate postnatal changes in cardiac amino acid metabolism, fatty acid synthesis, mevalonate pathway (cholesterol synthesis), and ketogenesis in more detail.

BCAA Catabolism
To gain deeper insight into the significant enrichment of the KEGG pathway "valine, leucine, and isoleucine degradation" in transcriptomics and proteomics clusters showing upregulation throughout the postnatal period, we focused on BCAA catabolism in more detail. The concentrations of valine, leucine, and isoleucine increased from P01 to P09, after which they decreased to P01 levels or lower at P23 ( Figure 6A). Most enzymes in the BCAA degradation pathway were significantly upregulated at P23 compared with P01 on either transcript or protein levels, or both ( Figure 6B). On the mRNA level, most differentially expressed genes were upregulated between P09 and P23, thus correlating directly with the BCAA concentrations. Of the 37 quantified individual proteins on this pathway, 19 and 31 proteins were upregulated at P09 and P23, respectively, compared with P01 (set 2 samples; Appendix S3). The rate-limiting step of BCAA catabolism is mediated by the branched-chain a-ketoacid dehydrogenase (BCKDC or BCKDH) complex, the activity of which is controlled by inactivating phosphorylation and activating dephosphorylation. Both the a subunit of the complex, Bckdha, and the protein phosphatase responsible for BCKDH activation (Ppm1k, protein phosphatase, Mg2+/Mn2+ dependent 1K) were significantly upregulated on the mRNA level at P23, and BCKDHA protein abundance was upregulated at P23 (PPM1K was not detected; Table S2).

Fatty Acid Metabolism
The main enzymes regulating the abundance of free fatty acids are FASN (fatty acid synthase); ACSLs (long-chain acyl-CoA [acyl-coenzyme A] synthetases), which attach CoA to free fatty acids; and ACOTs (acyl-CoA thioesterases), which remove the CoA from acyl-CoA, thus releasing free fatty acid ( Figure 7A). The relative mRNA levels of FASN, ACOTs, and ACSLs are presented in Figure 7B. Downregulation of Acot1, Acot2, Acsl3, Acsl4, and Fasn was observed after P01, whereas Acot13 was upregulated at P23. In general, there was relatively poor correlation between protein ( Figure 7C) and mRNA expression. The abundance of ACOT isoforms increased, except for ACOT1, which was first upregulated at P04 and P09, followed by downregulation at P23. The abundance of ACSL1 increased over the postnatal period. FASN, however, was strongly downregulated at P23 to undetectable levels, correlating with the mRNA expression pattern. The products of the cytosolic ACOT1 and mitochondrial ACOT2, saturated and monounsaturated medium-to-long-chain fatty acids, 41 exhibited almost identical patterns with an initial increase peaking at P04 or P09, followed by strongly diminished levels at P23 ( Figure 7D) and thus correlating with ACOT1 abundance. Furthermore, most enzymes in the KEGG pathway "fatty acid degradation" were differentially expressed on mRNA and/or protein levels ( Figure S10). These data reflect the complex regulation of fatty acid levels in the postnatal heart in response to the increased fatty acid abundance from nutrients.

Mevalonate Pathway and Ketogenesis
Because HMGCS2 (HMG-CoA [hydroxymethylglutaryl-CoA] synthase 2) was the most significantly upregulated protein from P01 to P04 and cholesterol biosynthesis was one of the enriched metabolism-related biological processes in GO enrichment analysis, the HMGCS-mediated mevalonate pathway and ketogenesis were investigated in more detail. The HMGCS-catalyzed synthesis of HMG-CoA serves as a substrate for HMGCR (HMG-CoA reductase) in the mevalonate pathway and HMGCL (HMG-CoA lyase) in ketogenesis ( Figure 8A). Both the cytosolic Hmgcs1 and the mitochondrial Hmgcs2 isoforms were downregulated at the mRNA level at P09 and P23 compared with P01, whereas there was no change in the mRNA levels of the upstream enzyme ACAT1 (acetyl-CoA acetyltransferase; Figure 8B). Of the enzymes regulating ketone body abundance, expression of Bdh1 was significantly lower at P09 compared with P01 or P23, and the expression of Oxct1, which oxidizes ketone bodies, was upregulated after P09. In the mevalonate pathway, Hmgcr expression was downregulated at P09 compared with P01 or P04 and at P23 compared with P04. Similarly, Mvk (mevalonate kinase), Pmvk (phosphomevalonate kinase), and Idi1 (isopentenyl-diphosphate d-isomerase 1) were downregulated with increasing postnatal age. Of these enzymes, ACAT1, HMCSS2, HMGCL, and OXCT1 (3-oxoacid CoA-transferase 1) were reliably quantified in proteomics ( Figure 8C). The abundances of ACAT1 and HMGCL increased throughout the early postnatal period, whereas HMGCS2 peaked at P04 and was downregulated to undetectable levels by P23 and OXCT1 was upregulated from P09 to P23. At the metabolite level, the end product of ketogenesis, 3-hydroxybutyric acid (b-hydroxybutyrate), peaked at P09 and was downregulated at P23 ( Figure 8D), correlating with changes in the abundance of HMGCS2, which is the rate-limiting enzyme of ketogenesis, and OXCT1. Metabolites of the mevalonate pathway were not detected or identified; however, the concentrations of cholesterol, which is produced from mevalonate, increased at P09 and decreased thereafter ( Figure 8D). Collectively, these results show that the HMGCS-mediated mevalonate pathway and ketogenesis are activated transiently after birth in the mouse heart. To evaluate the role of mevalonate pathway and ketogenesis in cardiomyocyte proliferation, we investigated the effects of pharmacological HMGCS and HMGCR inhibition on the viability and proliferation of neonatal rat ventricular cardiomyocytes. The HMGCS inhibitor hymeglusin 42 induced a %15% decrease in cell viability at 100 lmol/L but had no effect at lower concentrations ( Figure 8E and 8F). However, simultaneous inhibition of the mevalonate and ketogenesis routes with hymeglusin decreased the percentage of BrdU-positive cardiomyocytes in both serum-free and serum-stimulated conditions already at smaller, nontoxic concentrations (10 and 1-10 lmol/L, respectively; Figure 8E and 8F). Inhibition of the mevalonate pathway alone using the HMGCR inhibitor simvastatin, however, had no effect on the percentage of BrdU-positive cardiomyocytes at nontoxic concentrations ( Figure 8E). These results indicate that the HMGCS-mediated ketogenesis may participate in regulating cardiomyocyte cell cycle activity.

Discussion
Therapeutic strategies to promote regeneration of the adult human hearts are urgently sought. 12,13 Nevertheless, more detailed understanding of the metabolic changes and signaling pathways mediating cardiomyocyte maturation and cell cycle withdrawal is required for the development of regenerative therapies. In this work, we employed an integrated multiomics approach to investigate the metabolic changes occurring in the mouse heart within the early postnatal period to identify metabolic pathways associated with the postnatal loss of regenerative capacity. This study provides an important resource for molecule (RNA, protein, and metabolite) abundances in the neonatal mouse heart. Furthermore, we highlighted examples of metabolic pathways that exhibited correlative changes on all 3 levels. According to previous reports and the present data, the cardiac energy metabolism changes drastically after birth in response to altered nutrient availability and transition to oxygen-rich environment. Increased oxidative metabolism gives rise to reactive oxygen species causing oxidative stress and DNA damage, which is thought to contribute to cardiomyocyte cell cycle withdrawal. 19,43 Remarkably, exposure of adult mice to chronic hypoxia reduces oxidative metabolism and DNA damage and promotes cardiac regeneration after MI. 44 Oxidative DNA damage does not, however, correlate directly with cardiomyocyte cell cycle withdrawal in humans. 45 The fact that human pluripotent stem cell-derived cardiomyocytes continue to respond to mitogenic stimuli despite a normoxic environment may also indicate that other mechanisms contribute to the irreversible cell cycle withdrawal. Nevertheless, the metabolic switch from glycolysis to fatty acid oxidation, achieved by increased palmitic acid availability and insulin depletion, was recently reported to be sufficient for inducing irreversible cell cycle exit of human pluripotent stem cell-derived cardiomyocytes, indicating that oxidative metabolism plays a major role in regulating cardiomyocyte cell cycle. 46 Furthermore, increased fatty acid abundance has been shown to mediate postprandial physiological cardiac hypertrophy in the Burmese python. 47 Administration of a combination of myristic, palmitic, and palmitoleic acid to mice or pythons also induces cardiac hypertrophy without pathological fibrosis or activation of the fetal gene program. In this study, we showed a temporally regulated postnatal increase in the abundance of saturated and monounsaturated medium-chained fatty acids-both as free fatty acids and in various lipid species such as ceramides, sphingomyelins, and phosphocholines. This transient increase could serve as the physiological mechanism regulating nonpathological cardiomyocyte hypertrophy during postnatal heart growth, as reported for the Burmese python. 47 It is also tempting to speculate that the same fatty acids may play a role in driving postnatal cardiomyocyte cell cycle withdrawal, as reported for human pluripotent stem cell-derived cardiomyocytes. 46 Another key finding of this work, not previously described in the context of postnatal heart development, is the temporal regulation of the mevalonate pathway, which was strongly activated immediately after birth and downregulated after the regenerative window. The mevalonate pathway plays a role in cancer cell proliferation and is upregulated by several oncogenic signaling routes. 48 Furthermore, its downregulation has been linked to increased cell size in vivo, 49 and inhibition of HMGCR, the rate-limiting enzyme of mevalonate pathway, attenuates cell proliferation and increases cell size in various cell types in vitro. 50 Under our experimental conditions, specific inhibition of the mevalonate pathway with simvastatin, however, had no effect on neonatal cardiomyocyte proliferation, indicating that the active mevalonate pathway is dispensable for cardiomyocyte proliferation.
Parallel to the mevalonate pathway, we observed transient postnatal activation of ketogenesis in the mouse heart. The circulating levels of b-hydroxybutyrate, which is mainly produced in the liver, increase temporarily after birth and provide an important source of energy for the developing brain. 51 However, the observed strong temporal upregulation of HMGCS2, the rate-limiting enzyme of ketogenesis, indicates that ketone body synthesis is also locally regulated and postnatally activated in the myocardium. In addition to its role as a circulating energy source, b-hydroxybutyrate participates in cellular signaling by acting on cell membrane receptors and by directly inhibiting histone deacetylases, 52 which has been reported to suppress oxidative stress. 53 Increased abundance of ketogenic enzymes has been linked to aggressiveness of prostate cancer, 54 indicating that augmented ketogenesis may provide a proliferative advantage in certain conditions. We further showed that simultaneous inhibition of the mevalonate pathway and ketogenesis attenuates neonatal cardiomyocyte proliferation in vitro. Because inhibition of the mevalonate route with simvastatin had no effect on cardiomyocyte proliferation, our results suggest that the observed transient postnatal upregulation of ketogenesis in the postnatal mouse heart may participate in regulating the cardiomyocyte cell cycle. It is tempting to speculate that by inhibiting HMGCS, hymeglusin could reduce the production of b-hydroxybutyrate, which is known to suppress oxidative stress through inhibition of histone deacetylases, and thereby increase oxidative stress and DNA damage and promote cardiomyocyte cell cycle withdrawal. Clarification of the exact mechanisms requires further investigation.
Unlike other amino acids, the BCAAs valine, leucine, and isoleucine are metabolized mainly in organs other than the liver, such as skeletal and cardiac muscle. They provide nitrogen for maintaining glutamate, alanine, and glutamine pools and function as signaling molecules activating mTOR (mammalian target of rapamycin) signaling, which regulates cardiac homeostasis and plays a crucial role in cardiac pathophysiology. 55,56 In addition, isoleucine has also been reported to inhibit the transport and utilization of fatty acids in skeletal muscle in mice. 57 The observed increase in BCAA degradation in the postnatal heart is thus in line with the decreased protein synthesis at P23 and the increase in fatty acid metabolism after P01. Even though the abundance of BCAAs is altered in experimental models of pressure overload and MI, 58  cardiac regeneration, as they predominantly take place after the regenerative window. With respect to methodological limitations, proteomic and metabolomic analyses can detect only a fraction of proteins and metabolites. By using shotgun proteomics, we were able to perform relative quantification of >2000 proteins. Hydrophobic (eg transmembrane) proteins and proteins with very low abundance, for example, were often below detection. To increase metabolite coverage, we applied 2 complementary MS-based methods for the metabolomics analyses. Metabolite identification was confined to previously reported metabolites in spectral libraries; therefore, numerous detected metabolic features with interesting abundance patterns remain unidentified. In general, the correlation between mRNA and protein levels is only modest, which can be caused by, for example, variable protein turnover rates. Integration of protein turnover analysis to transcriptomics and proteomics was recently shown to increase the yield of identified disease gene candidates in pathological cardiac hypertrophy in mice, 59 highlighting the benefit of adding further layers to omics-based analyses. Qualitative and quantitative assessment of all biological processes has been deemed essential for the investigation of cardiac metabolism by the American Heart Association. 60 In the present study, integration of 3 omics analyses with 4 time points over the early postnatal period provides a comprehensive view to the metabolic changes taking place in the developing heart. Interestingly, in some cases we observed surprisingly long delays in the abundance changes from mRNA to protein and from protein to metabolite, as exemplified by the mevalonate pathway and ketogenesis. This highlights how crucial it is not only to use several omics methods but also to include several time points when analyzing dynamic phenomena with integrative systems biology approaches.
Considering the multicellular composition of the heart, the present study cannot elucidate cardiomyocyte-specific phenomena because the data represent tissue-level changes. In cell volume, cardiomyocytes build up 70% to 80% of the cardiac tissue, 61 whereas in cell numbers, cardiomyocytes constitute roughly one third of the cell population, and endothelial cells represent the most abundant cell type, with a %45% proportion. 62 These 2 cell types play a crucial role in postinfarction cardiac regeneration and, unlike fibroblasts or leukocytes, fail to initiate the neonatal-type gene expression response upon injury in adult hearts. 25 In cardiomyocytes, this was suggested to result from chromatin inaccessibility and thus to represent an epigenetic roadblock that prevents cardiomyocyte cell cycle reentry. Recent evidence suggests that metabolic pathways are intimately connected with epigenetic regulation. 63,64 Whether the postnatal metabolic remodeling is in fact the driving force that induces cardiomyocyte cell cycle exit requires further investigation.
In summary, we used an integrative systems biology approach with 3 levels of omics analyses to characterize changes in molecule abundances throughout the early postnatal period and to identify metabolic pathways associated with cardiac regeneration. To our knowledge, this study is the first combining transcriptomics with untargeted proteomics and global metabolomics analyses over several time points in the early postnatal heart and, as such, provides an extensive resource for molecule abundances for future mechanistic studies. Furthermore, we present several examples of metabolic pathways with correlative changes on all omics levels. These include the well-established metabolic switch from glycolysis to fatty acid b-oxidation as well as many previously unreported changes in cardiac metabolic pathways, such as the mevalonate pathway and ketogenesis. Finally, we identified a biological function for mevalonate and ketone body metabolism in the heart with a potential role in the regulation of neonatal cardiomyocyte proliferation. This integrated molecule-level data may open up new possibilities for the development of regenerative therapies. Data S1.

Preparation of tissue samples
Mouse pups (strain C57BL/6JOlaHsd, both sexes) were acquired from the experimental animal facility of the University of Helsinki with an internal use licence (KEK14-014) and were used on postnatal days 1, 4, 9 and 23 (P01, P04, P09 and P23, respectively) before weaning. Animals were housed in standard conditions and their handling and all procedures were carried out in accordance with University of Helsinki institutional guidelines, which conform to the National Research Council (US) Guide for the Care and Use of Laboratory Animals. 1 The tissue samples were collected at approximately midday without fasting. All P01-P09 pups had been suckling recently and had milk in their stomachs. Although not weaned, the P23 pups had access to general mouse feed in addition to milk. A total of 92 mice were decapitated, and their ventricles were excised and cut into 3-5 pieces, rinsed in phosphate-buffered saline, snap-frozen in liquid N2 and stored at -80 °C.

RNA extraction and purification
Ventricular tissue samples from nine mice for each time point from set 2 were randomly pooled together (3 samples / pool) to yield three samples for each time point. The tissue pieces were homogenized in TriZol ® (Thermo Fisher Scientific) using Precellys soft tissue beads (Bertin Instruments, Montigny-le-Bretonneux, France) and Precellys24 Homogenizer (Bertin Instruments). Total RNA was extracted using chloroform extraction, precipitated from the aqueous phase with isopropanol, washed using 75% ethanol and dissolved in RNAse free water, whereafter it was purified using RNeasy ® Mini Kit (Qiagen) according to manufacturer's instructions. The integrity and quality of the purified RNA was analysed using Agilent 4200 TapeStation (Agilent, Santa Clara, CA, USA).

RNA sequencing
The RNA library for sequencing was prepared with SureSelect Strand-Specific RNA Library Prep Kit (Agilent). Ribosomal RNA was depleted from 1 µg of total RNA using Ribo-Zero Complete Gold rRNA Removal Kit (Illumina, San Diego, CA, USA), whereafter the rRNA-depleted RNA was purified with RNeasy mini Elute columns (Qiagen) and fragmented to generate approximately 300 bp segments. Double-stranded cDNA (ds-cDNA) was then synthesized, 3' ends of the ds-cDNA were adenylated and adaptors ligated and the ds-cDNA was purified with Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA). The produced ds-cDNA library was analysed for quantity using Qubit ® dsDNA HS kit (Thermo Fisher Scientific) and Qubit ® Fluorometer (Thermo Fisher Scientific) and for quality with 2100 Bioanalyzer ® (Agilent) using a high-sensitivity chip. Sequencing was then performed as single-end sequencing for read length of 75 bp using Illumina NextSeq Sequencer (Illumina) in one high-output run.

Analysis of RNA sequencing data
Quality of the sequencing data was analysed using FastQC 2 and was considered excellent with average Phred Quality Scores of >30. Quality trimming was applied to the data with the Trimmomatic software. 3 The sample reads were then aligned against the GENCODE M12 (GRCm38.p5) reference genome. 4 STAR-aligner 5 was used in genelevel output mode to produce sorted alignment map files, from which alignment statistics were collected with Qualimap 2. 6 The counts were produced from the STAR output data with the featureCounts software 7 using strand specificity option of 2 due to the reverse-stranded sequencing library. The successfully assigned reads varied from 2x10 6 to 12x10 6 reads per sample, which was considered adequate as the duplication level was low.
The quality of the data and differential expression were analysed with DESeq2 software 8 in R 9 environment. The count values were normalized using a geometric mean and sample-wise size factors were estimated to correct for variation in library size. The dispersion of gene-wise values between sample groups was estimated and negative binomial linear model and Wald test were used to produce p values. To optimize p value adjustment, removal of lowexpression results and outliers (by Cook's distance) was carried out, whereafter multiple testing adjustment of p values was performed using the Benjamini-Hochberg procedure. Gene ontology (GO) 10,11 enrichment analysis was performed with ranked gene lists using GOrilla 12 (available online at http://cbl-gorilla.cs.technion.ac.il).

Tissue homogenization for protein and metabolite extraction
The ventricular tissue samples (Figure 1) were homogenized using a probe sonicator (Rinco Ultrasonics, Arbon, Switzerland) for sample set 1 and FastPrep-24 5G bead homogenizer (MP Biomedicals, Santa Ana, CA, USA) with 1 mm zirconia beads (BioSpec Products, Bartlesville, OK, USA) for sample set 2. The homogenization was done in Milli-Q water and the volume was adjusted according to tissue weight (approximately 1 µL of water per 0.1 mg of tissue, in total 80 -662 μL). After homogenization, samples were divided into two for metabolomics and proteomics experiments, and stored at -80 °C freezer prior to other pre-treatment steps. Total protein content of the homogenates was measured five times with a Thermo NanoDrop 2000 (Thermo Scientific, Waltham, MA, USA).

Proteomics
The chemicals and solvents were from Sigma-Aldrich (Steinheim, Germany) unless otherwise mentioned. Water was purified with a Milli-Q water purification system (Millipore, Molsheim, France). For liquid chromatography (LC) running buffers, LC-MS grade water was used.
Cell lysis and protein denaturation were done by sonication in 8 mol/L urea (set 1: 200 μL; set 2: 360 μL; Amresco, Solon, OH, USA). Insoluble cell debris was removed by centrifugation at 21000 g for 15 min at room temperature, and the samples were diluted to < 1.5 mol/L urea with water (set 1: 900 μL; set 2: 1640 μL). Cysteine residues were reduced with 50 mmol/L dithiothreitol (set 1: 110 μL; set 2: 222 μL; final concentration 5 mmol/L, 20 min incubation in 60 °C) and carbamidomethylated with 150 mmol/L iodoacetamide (set 1: 121 μL; set 2: 246 μL; 20 min incubation in room temperature in dark), after which the pH was adjusted to 7-8 with 200 μL of 1 mol/L ammonium bicarbonate (Fluka Chemie AG, Buchs, Switzerland), and the proteins were digested with 2.5 μg of sequencing grade modified trypsin (Promega, Madison, WI, USA) overnight at 37 °C with shaking. The resulting peptide solution was acidified to pH < 2 with 200 μL of 10 % trifluoroacetic acid and purified using C18 MicroSpin columns (The Nest Group, Southboro, MA, USA). The columns were conditioned with 3 x 100 μL of acetonitrile, equilibrated with 4 x 100 μL of 1 % acetonitrile and 0.1 % trifluoroacetic acid in water, after which the samples were loaded into the columns in 400 μL aliquots. The washing was done with 4 x 100 μL of 5 % acetonitrile, 0.1 % trifluoroacetic acid in water and the elution with 3 x 100 μL of 50 % acetonitrile and 0.1 % trifluoroacetic acid in water. After the C18 purification, the samples were vacuum centrifuged to dryness and stored in -20 °C.
Before the liquid chromatography-mass spectrometry (LC-MS) analysis, the samples were recovered by sonication in 60 μL of 1 % acetonitrile and 0.1 % trifluoroacetic acid in water. The LC-MS run order was randomized, and pooled quality control samples were used for monitoring system stability throughout the analysis. The injection volume was 4 μL (equivalent to approximately 4 μg of total digested protein, determined based on the total protein content measurement of the homogenate).
Protein identification and label-free quantification were done with Andromeda search engine and MaxQuant software, respectively. [13][14][15] One percent false discovery rate (FDR) filtering was applied on both peptide and protein level. An up-to-date mouse UniProtKB/Swiss-Prot proteome (16,800 proteins restricted to canonical entries, downloaded 2016-08-09) was used as the library, 16 with search mass tolerance 6 ppm in MS1 and 20 ppm in MS2; a maximum of two missed cleavages was allowed. Cysteine carbamidomethylation was set as a static modification, and methionine oxidation as a variable modification. Contaminant and decoy matches were omitted before further data analysis, which was carried out with R 9 using the normalized MS1 intensity (LFQ intensity) values from MaxQuant without further normalization. Protein grouping was applied, and only unique and razor peptides were used for quantification. In addition to base R, the packages ggplot2 17 and FactoMineR 18 were used. The bioinformatic analysis was performed with DAVID 6.8. 19,20 Metabolomics

Sample preparation for GCxGC-MS analysis:
Tissue homogenates were thawed on ice and 25 µL of homogenate was transferred into a new Eppendorf-tube. Methanol (410 µL) containing labeled internal standards (0.9-4 µg/mL) was added to the samples to precipitate the proteins. The samples were vortexed, sonicated for 5 min, and incubated on ice for 30 min, whereafter the extracts were centrifuged for 5 min (14 338 g, 4°C) and 180 µl aliquots of supernatants were transferred into vials and evaporated to dryness under nitrogen flow. The metabolites were then converted into their methoxime and trimethylsilyl (TMS) derivatives by automated two-step derivatization. First, 25 μL of MOX was added to the residue, and the mixture was incubated for 1 h at 45 °C. Next, 25 μL of MSTFA including the retention index standards was added, and the mixture was incubated for 1 h at 45 °C. Finally, hexane including injection standard (50 µL) was added to the mixture immediately prior to injection. Quality control (QC) samples were prepared by pooling P1, P4, and P10 tissue samples, which were prepared similarly as described above. All metabolomics samples and pooled QCsamples were randomized before the sample pretreatment and the analysis.
GCxGC-MS analysis of polar metabolites: Polar metabolites were analyzed by GC×GC-MS at Steno Diabetes Center (Gentofte, Denmark) as described earlier. 21 The analyses were performed with a LECO Pegasus 4D instrument consisting of an Agilent 7890 gas chromatograph equipped with a split/splitless injector (Agilent Technologies, Santa Clara, CA), a cryogenic dual-stage modulator and a time-of-flight mass spectrometer (Leco Corp., St. Joseph, MI, USA). Instrumentation included a multipurpose sampler from Maestro software (Gerstel, Mülheim an der Ruhr, Germany), which was used for derivatization and sample introduction. A 10 m × 0.18 mm I.D. RTXi-5MS column with film thickness of 0.2 μm (Restek Corp., Bellefonte, PA, USA) was used as the first column and a 1.5 m × 0.1 mm I.D. BPX-50 column with a film thickness of 0.1 μm (SGE Analytical Science, Austin, TX, USA) as the second column. A deactivated fused silica retention gap column (1.7 m × 0.53 mm I.D., Agilent technologies, Santa Clara, CA) was connected in front of the first column. The injector temperature was 240 °C and the samples (1 μL) were injected using pulsed splitless mode with a splitless period of 90 s. High-purity helium (99.9999%, Yara Praxair) was used as the carrier gas in a constant-pressure mode with initial pressure of 40 psig. The first column oven temperature program was as follows: 50 °C in isothermal mode for 2 min, from 50 °C to 240 °C at 7 °C/min, from 240 °C to 300°C at 25 °C/min, and kept at 300 ° for 3 min. The second dimension column oven temperature was maintained constantly 20 °C higher than the first column oven and the modulation time was 4 s. The temperature of the transfer line was 260 °C. The temperature of the electron impact ion source was 200 °C and the electron energy 70 eV. The mass range was from 45 to 700 m/z and 100 spectra/s were measured.

LC-MS Analysis:
All experiments were performed with a Thermo Orbitrap Fusion high resolution mass spectrometer equipped with a Thermo Dionex Ultimate 3000 UHPLC (both from Thermo Fisher Scientific). The column was Waters Acquity BEH C18 (2.1 x 100 mm, 1.7 µm). The separation was performed using a 15-min linear gradient from 95% A (0.1% formic acid in Milli-Q water) to 100% B (0.1% formic acid in methanol) followed by a 10-min isocratic period of 100% B and a 5-min equilibration to initial conditions. For the sample set 2, a minor modification was done for the LC method so that the starting eluent was 100% A (0.1% formic acid and 5% methanol in MQ water). The flow rate was 0.3 mL/min, total run time 30 min, the column temperature 25 °C, and the sample tray was cooled to 10 °C. The injection volume was 2 µL for the first sample set and 3 µL for the second sample set. In the FTMS full scan the mass resolution was 120 000 FWHM and the mass range from m/z 100 to m/z 1000. The spectra were collected in centroid mode with a duty cycle time of 0.6 s, applying internal calibration (Easy-IC with fluorantine). Heated electrospray, with a spray voltage of 3.0 kV in the positive ion mode was applied for ionization. The other source parameters were: sheath gas 30 arb, aux gas 10 arb, ion transfer tube temperature 333 °C and vaporizer temperature 317 °C. For identification purposes MS/MS spectra were collected applying quadrupole isolation window of 1, HCD fragmentation energies of 25% and 40% and mass resolution of 30 000 FWHM.
GCxGC-MS data processing and analysis: ChromaTOF version 4.71 (LECO Corporation, St. Joseph, MI) was used for the raw data processing. Automatic peak detection and mass spectrum deconvolution were performed using a peak width of 0.2 s and a peak signal-tonoise (S/N) value threshold of 100. The data files obtained by the ChromaTOF software were exported to text files and Guineu software 22 was used for the compound alignment. The column bleed was removed and only the TMS-derivatized peaks (73 m/z ion in the spectra) were selected for further data processing. The peak areas from total ion chromatograms (TIC) were used for most of the compounds. For a pre-selected set of compounds (see Supplemental results, Table 2), the peak areas of the selected ions were monitored and quantified based on external calibration curves. The linear retention indices (RIs) were calculated based on the retention times of the retention index standards (n-alkanes). The alignment of the data was performed based on RIs, first and second dimension retention times as well as electron ionization (EI) spectral matches. Only the metabolites detected in more than 10% of the analyzed samples were included in the dataset.
The metabolites were identified tentatively by comparing the measured mass spectra and the RIs to those in NIST 2014 or in-house libraries. In addition, Golm metabolome database 23 and Fiehn library 24 was used for the further identification of interesting features. Identification annotations are based on the recommended identification criteria by metabolite standard initiative. 25 Identification level 1 is annotated to compounds identified against standard compounds. Level 2 identifications are annotated to compounds, whose spectral match is > 850 and RI difference < 30 compared to respective values in NIST 2014 or in-house library. Level 3 is compound class specific annotation for the findings, whose spectral match is >800 against NIST 2014 library and are additionally searched against other library (Golm metabolome database or Fiehn library). Level 4 spectral match is below these criteria and has been annotated as unknown.
Data processing was performed with R 3.1.2. 9 The data were filtered further with the criteria that the peak should be found in 75% of the samples in at least one of the sample groups (P01, P04, P09, P23). Principal component analysis (PCA) as well as t-tests with FDRcorrection (Benjamini and Hochberg) were performed to both sample sets separately. The features in the two separate sample sets were aligned and combined according to maximum dot product of spectra, retention index difference ≤ 10 and second retention time ≤ 0.25 s, respectively. Combined data was explored with linear mixed effect model (LME), 26 where the first variable was the postnatal age and the second was the dataset. P01 estimate was set to 0 with variance one and other sample groups were compared to the P01.
LC-MS data processing and analysis: MzMine2 27 was used for the data processing, and the following functions were applied: mass detection (centroid algorithm), chromatogram builder, chromatogram deconvolution (local minimum search algorithm), deisotoping, alignment (joint aligner), filtering (number of found), and cap filling. Peak intensity was normalized to the intensity of internal standard LysoPC(17:0) and to the mass of the tissue sample. Differences between the sample groups of the individual data sets were searched with one-way ANOVA followed by Tukey`s HSD post-hoc test using R and metadar package. Fold change values were calculated and p-values adjusted with the FDR correction (Benjamini and Hochberg). Features were also combined with LME model, with the following criteria: a maximum mass difference of 2 ppm and a maximum retention time difference of 0.35 min. If multiple matching features were found with several retention times, correct alignments were confirmed manually. MS/MS spectra were collected for all the features from the sample set 1 showing statistical difference (q < 0.01) in ANOVA. The features showing significant differences also with LME model were identified by searching the measured spectra against existing spectral databases or repositories; mzCloud (https://www.mzcloud.org/), NIST 2014 MS/MS (http://chemdata.nist.gov/), HMDB 28 , MoNA (http://mona.fiehnlab.ucdavis.edu). In addition, in silico database LipidBlast 29 and molecular structure database search tool CSI:Finger ID 30 or known group-specific fragmentation of lipid species was used in the identification. Phospholipids were identified and annotated based on the exact mass information and known polar head group fragments of phospholipids. These are annotated with phospholipid class, the number of carbons and double bonds in the fatty acid chains (e.g. PC(32:0)).
Metabolite set enrichment analysis (MSEA) 31 was performed to all significantly changed (q<0.01) identified metabolites in LME model. Pathway-associated metabolite set library was applied and all compounds from the reference metabolite library were selected for the enrichment analysis.

Fuzzy clustering
To compare and integrate the results of the three different omics analyses, unsupervised pattern classification was carried out using fuzzy c-means clustering, [32][33][34] in which the transcripts, proteins, and metabolites were assigned into one or more clusters based on their abundance pattern. The analysis was carried out using the C-means clustering package in R 9 with fixed numbers of clusters (7 for transcriptomics and proteomics, 5 for metabolomics) and the fuzziness value was determined according to the method described in 34 . The input data for fuzzy c-means clustering was limited to transcripts, proteins, and metabolites that were significantly different in abundance in one or more of the comparisons in sample set 2.
Fuzzy clusters were analyzed with DAVID 6.8 19,20 for GO term 10,11 and KEGG 35 pathway enrichment using the recommended settings. Upstream regulator analysis was performed with IPA (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathwayanalysis) with the standard settings. Upstream regulators of the classes "transcription regulator", "microRNA", and "chemicalendogenous mammalian" with p < 0.01 were included in the analysis.
In order to analyze cell viability, the cardiomyocyte cultures were treated with hymeglusin or simvastatin (from Santa Cruz and Sigma-Aldrich) at concentrations of 1-100 µmol/L (µM) for 24 h in CSFM and cell viability was analyzed using the 3-(4,5-dimethylthiazol-2-yl)-2,5diphenyltetrazolium bromide (MTT) assay as described previously. 37 Briefly, 0.5 mg/mL MTT was added to the medium for 2 hours (incubation at 37 °C, 5% CO2, humidified atmosphere), whereafter the medium was removed and formazan crystals dissolved in DMSO (200 µL/well). The absorbance was then measured at 550 nm with absorbance at 650 nm subtracted as background. Data are presented as percentage of control, mean+SEM from three independent experiments carried out in triplicate.
Statistical analysis of cell viability and proliferation data was carried out with IBM SPSS Statistics version 24 using Welch ANOVA followed by the Games-Howell post hoc test. The level of statistical significance was set at p < 0.05.

GCxGC-MS analysis:
The repeatability of the analytical method was evaluated by determining the peak areas of the internal and injection standards (Table 1) added to all heart samples (n=44 in sample set 1 and n=39 in sample set 2) and the peak areas of the quantified compounds (Table 1, Figure  2) in pooled heart samples (n=10 in sample set 1 and 2) used as quality control (QC) samples. The samples were pretreated and analyzed with GCxGC-MS in random order and quantified using external calibration. The relative standard deviation (RSD%) values for the internal and injection standards were between 8 and 21% indicating good repeatability of the GCxGC-MS method. The repeatability of the quantified compounds was typically below 30% (RSD%), which is considered acceptable in metabolomics analysis. The metabolites in the heart tissue samples were analyzed as their TMS derivatives and the peak areas were used in determining the relative abundances of the detected features. After the first data preprocessing steps, the GCxGC-MS datasets of the sample set 1 and 2 contained 760 and 771 metabolic features, respectively. After filtering, the total numbers of features were 347 and 443. These features were further analyzed with t-test applying FDR correction (Supplemental Dataset S4). The data from the two separate sample sets were combined based on spectral alignment, retention index and second retention time difference. The combined dataset resulted in 328 features, which were analyzed with LME model resulting in total 162 metabolic features with statistical significance (q <0.01) in at least one of the group comparisons (P01 vs. P04, P01 vs. P09, P01 vs. P23; Supplemental Dataset S4).
Principal component analysis of all samples and pooled quality control samples from the sample set 2 shows some grouping of the samples based on postnatal age (Figure 2). The P01, P04, and P09 sample groups cannot be totally distinguished from each other with the two principal components. However, P23 is totally separated from the other sample groups and the QC samples clustered among all sample groups in the individual factor maps plot ( Figure 2). Even though the QC samples show some variation in the direction of principal component 1, the variation in the direction of principal component 2 is explained more with the postnatal age.

LC-MS analysis
The analytical performance of the LC-MS method using the orbitrap mass spectrometer was evaluated in terms of repeatability of peak intensities (peak area), retention times and mass accuracies of the internal standards added to the heart tissue samples (n=44 in sample set 1 and n=39 in sample set 2). The samples were pretreated and analyzed in random order. The RSD% for the peak intensities without normalization and after normalization against LysoPC(17:0) were below 25%, which is considered acceptable in metabolomics. Furthermore, no systematic trend in signal intensities were observed. The mass errors between exact mass (calculated) and accurate masses (mean, measured) were below 1 ppm indicating excellent mass accuracies. The differences of mass accuracies between the sample sets 1 and 2 were below 0.1 ppm and those of retention times below 0.25 min. Based on these values the maximum mass difference of 2 ppm and maximum retention time difference 0.35 min were used as limiting values for combining the heart tissue sample sets 1 and 2 with LME model. After preprocessing, the LC-MS analysis of the sample sets 1 and 2 comprised 5591 and 9043 metabolic features, respectively. The detected features from the sample sets were analyzed separately with ANOVA. The q values and fold changes are presented in Dataset S4. The features from individual datasets were further combined based on retention time and mass-to-charge ratio difference observed with the internal standards in the sample sets 1 and 2 (Table 2) and the combined dataset containing 3398 features was analyzed with the LME model. In total 805 metabolic features displayed statistical significance (q <0.01) in at least one of the group comparisons (P01 vs. P04, P01 vs. P09, P01 vs. P23; Dataset S4).   Statistically significant changes (q < 0.01) are indicated in bold Upregulated Downregulated q-value < 0.01 q-value < 0.05 Not significant log2 fold q-value log2 fold q-value log2 fold q-value log2 fold q-value log2 fold q-value log2 fold q-value log2 fold q-value log2 fold q-value log2 fold q-value

P17182 Eno1
Alpha-enolase -0. 21 Figure S1. Individuals factor map of the RNA sequencing data principal component (PC) analysis shows perfect separation of sample groups.
The genes were selected based on fold change among all statistically significantly changes genes (p < 0.01) and gene expression changes are presented as log2 fold change.
All gene symbol explanations are available in Dataset S1. A, Relative expression levels of cardiac ion channels in the mouse ventricular tissue within the postnatal period analysed. B, Relative expression levels of control genes in the mouse ventricular tissue within the postnatal period analysed. The data was normalised to P1 and is expressed as mean + SEM from 3 pooled samples (each from 3 hearts). All gene symbol explanations are available in Dataset S1.   For RNA sequencing and proteomics, respectively, transcripts or proteins detected in > 2/3 samples in at least one sample group were included in the clustering. For metabolomics, both LC-MS/MS and GCxGC-MS data of the metabolites with significantly different abundance in any comparison with LME model were included. The transcripts, proteins, and metabolites belonging to each cluster are presented in Dataset S5. Figure S6. The abundance profiles of all individual transcripts, proteins and metabolites from fuzzy c-means clustering of RNA sequencing, proteomics, and metabolomics data.