This article is mentioned in:

Introduction

Myocardial infarction is one of the most prevalent medical conditions worldwide (1). Pathological damage caused by myocardial infarction often results in irreversible injury to myocardial cells (2). This damage grossly affects cardiac function, due to the markedly decreased number of cardiomyocytes. The low capacity of adult cardiomyocytes to undergo proliferation and differentiation is the main reason why myocardial infarction cannot be repaired. The concept that hypertrophy, but not proliferation, of cardiomyocytes occurs under some pathophysiological conditions after birth has long been accepted (3). However, recent evidence has challenged this notion. Poss et al investigated myocardial generative capacity using a ventricular resection model, and observed that adult myocardial proliferation may exist in adult teleosts, including zebrafish (4). In addition, Porrello et al reported that excising partial cardiomyocytes at the apical portion of the heart in 24-h-old mice could stimulate cardiomyocytes to repair the damaged region (5), whereas Naqvi et al demonstrated that a surge of thyroid hormones during preadolescence may serve an important role in myocardial proliferation (6). These phenomena indicate that adult mammalian myocardial proliferation may be driven through physiological or pathological stimulation. Furthermore, Mollova et al observed that cardiomyocytes exhibited increased proliferative potential in 10-year-old children compared with in adults (7). Therefore, analyzing and investigating the gene expression profile of cardiomyocytes in the postnatal period may enhance our current understanding of the molecular mechanisms underlying the limited proliferative capacity of cardiomyocytes in myocardial growth arrest.

Gan et al conducted an integrative analysis of the developing postnatal mouse heart transcriptome by comparing 2- and 13-day-old postnatal mouse hearts (8), which reflect the dividing and growth-arresting phases of cardiomyocytes, respectively. The results revealed that GATA binding protein 4, myosin heavy chain 7 and insulin like growth factor 1 receptor (IGF1R) may act on the gene interaction network, whereas taspase 1, transducer of ERBB2, 1, chromosome 1 open reading frame 61, allograft inflammatory factor 1, Rho associated coiled-coil containing protein kinase 1, trefoil factor 2 and microRNA-503-5p were the key drivers inhibiting cardiomyocyte proliferation. Another study by Naqvi et al reported that mice exhibited a burst of cardiomyocyte proliferation during early preadolescent development, between 7 and 15 days of life, which was regulated by the thyroid hormone/IGF1/IGF1R/protein kinase B (AKT) pathway (6). Therefore, based on the physiological characteristics of mouse cardiomyocytes, 13-day-old mice may not fully reflect either the growth-arrested stage or the proliferative burst phase. To better elucidate the mechanisms and pathways involved in myocardial growth and arrest, the present study analyzed the differentially expressed genes in mouse hearts obtained from 24-h-old mice, which contain proliferative cardiomyocytes (5,6); 7-day-old mice, in which cardiomyocytes are undergoing a proliferative burst (6); and 10-week-old mice, which contain growth-arrested cardiomyocytes (6).

Furthermore, Gene Ontology (GO) annotation, pathway enrichment and series test of cluster (STC) were conducted to investigate the proliferative and growth-arresting phases of cardiomyocytes from numerous perspectives. Furthermore, the signal transduction pathway associated with the cell cycle was identified, in order to characterize the network involved in the phenomenon of myocardial cell cycle arrest. Elucidating the molecular-genetic alterations involved in the proliferative, proliferative burst and growth-arrest phases of cardiomyocytes may enable the identification of potential signaling factors that can induce the regeneration of cardiomyocytes in heart disease.

Materials and methods

Animal sample collection

Pregnant C57BL/6 mice (n=9; age, 10–11 weeks; weight, 25–30 g) were purchased from the National Institutes for Food and Drug Control (Beijing, China) and were quality inspected by the National Quality Inspection Center of Experimental Animals (Beijing, China). All mice were maintained at 22°C and 55% relative humidity under a 12-h light/dark cycle with free access to normal food and water. Mice were maintained according to the protocols of the National Science Council of the People's Republic of China (Beijing, China). In the present study, all of the mice used at 24 h (n=9; weight, 1–2 g) and 7 days (n=9; weight, 3.5–5 g) of age were delivered from 9 pregnant mice (10-week group) through eutocia. Each time-point included three parallel samples, and each sample consisted of three mouse hearts. The 18 hearts from 24-h- and 7-day-old mice were obtained from the 9 10-week-old pregnant mice. The mice were sacrificed by cervical dislocation following mild sedation with isoflurane. Whole hearts were harvested from 24-h-old, 7-day-old and 10-week-old mice. Ventricular myocardial samples were immediately dissected from the whole hearts for total RNA isolation following several washes with D-Hanks solution. The present study was approved by the Tianjin Medical University Animal Care and Use Committee (Tianjin, China), and all animal experiments were performed according to the Guidelines for Animal Experiments of the Tianjin Medical University, 2014 revision.

Total RNA isolation and gene expression analysis by microarray

RNA was extracted from the cardiac tissues using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA), according to the manufacturer's protocol, and subsequent purification was conducted using an RNeasy kit (Qiagen Sciences, Inc., Germantown, MD, USA). cDNA was generated using One-Cycle Target labeling and control reagents (first-strand cDNA was incubated in the thermocycler at 25°C for 1 h, 42°C for 1 h, 4°C for 2 min; second-strand cDNA was incubated in the thermocycler at 16°C for 1 h, 65°C for 1 min, 4°C for 2 min), and cRNA was generated using a GeneChip IVT labeling kit (both Affymetrix; Thermo Fisher Scientific, Inc.). Biotin-labeled, fragmented (≤200 nt) cRNA was hybridized for 16 h at 45°C to Affymetrix GeneChip Mouse Gene 2.0 arrays (Affymetrix; Thermo Fisher Scientific, Inc.). The GeneChips were washed and stained using the Affymetrix Fluidics Station 450 (Affymetrix; Thermo Fisher Scientific, Inc.), and were scanned using Affymetrix® GeneChip Command Console, installed in the GeneChip® Scanner 3000 7G (Affymetrix; Thermo Fisher Scientific, Inc.). The data were analyzed with the Robust Multichip Analysis (RMA) algorithm using Affymetrix default analysis settings and global scaling as the normalization method. Values are presented as log2 RMA signal intensity.

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

RT-qPCR was used to verify microarray results using the same RNA samples as used for the microarray analyses. The primer sequences used in the present study are listed in Table I. Briefly, total purified RNA was reverse-transcribed using the SuperScript III First-Strand Synthesis system (Invitrogen; Thermo Fisher Scientific, Inc.). RNA was incubated in the thermocycler at 44°C for 1 h, 92°C for 10 min to inactivate the reverse transcriptase. RT-qPCR was performed using SYBR-Green PCR mix and was run on the ABI 7300 Real-Time PCR system (both Applied Biosystems; Thermo Fisher Scientific, Inc.). The reactions were performed in triplicate in a total volume of 25 μl. The cycling conditions consisted of an initial, single cycle of 5 min at 95°C, followed by 40 cycles of 30 sec at 95°C, 30 sec at 54°C, 15 sec at 72°C, and fluorescence acquisition at 83°C for 1 sec. A melting curve analysis (60–95°C) was used to assess amplification specificity. The relative expression levels of each gene were determined using the 2−ΔΔCq method (9). Statistical analyses were performed as described for the microarray data. Pearson's correlation coefficient was further calculated for each gene using the normalized data to quantify the consistency between microarray experiments and RT-qPCR (P<0.05 and R>0.9).

Table I Primer sequences.

Table I

Primer sequences.

Gene	Primer sequence (5′ to 3′)
BMP1	F: ATCCAGGTGTCTGAGGGTTTCC
	R: GACTTCCAGGTAGTCGTAGGCA
CCNE2	F: GTACTGTCTGGAGGAATCAGCC
	R: CCAAACCTCCTGTGAACATGCC
IGF1	F: GTGGATGCTCTTCAGTTCGTGTG
	R: TCCAGTCTCCTCAGATCACAGC
CDC20	F: GATGGACGACATCTGGCAAGTG
	R: GTTGCCAGGATATTGGACTGCC
PLK1	F: CCATCTTCTGGGTCAGCAAGTG
	R: CCGTCATTGTAGAGAATCAGGCG
CDC25C	F: AACTCCACGACTCGGCAAACCT
	R: GTAGTAAGCGGAGAGGCAGACA
MYC	F: TCGCTGCTGTCCTCCGAGTCC
	R: GGTTTGCCTCTTCTCCACAGAC
CCNG1	F: GCGTTGGAGATCCAAGCACTGA
	R: GGAAACAAGCTCTTGCCAGAAGG
FGF1	F: CCAAGGAAACGTCCACAGTCAG
	R: ACGGCTGAAGACATCCTGTCTC

[i] BMP1, bone morphogenetic protein 1; CCNE2, cyclin E2; CCNG1, cyclin G1; CDC20, cell division cycle 20; CDC25C, cell division cycle 25C; F, forward; FGF1, fibroblast growth factor 1; IGF1, insulin like growth factor 1; MYC, MYC proto-oncogene, bHLH transcription factor; PLK1, polo like kinase 1; R, reverse.

Multi-class classification

The random variance model (RVM) F-test, which is commonly used to compare >2 groups, was applied to filter the differentially expressed genes for the control and experimental groups, since the RVM F-test effectively increases the degrees of freedom in cases of small samples. After analysis of significance and false discovery rate (FDR) analysis using BRB-ArrayTools version 3.0 [http://linus.nci.nih.gov/BRB-ArrayTools.html (10)], the differentially expressed genes were selected according to the P-value threshold (P-value <0.05, FDR <0.05, fold change ≥1.5) (10–12).

Cluster analysis

Genesis suite was used for cluster analysis (13). The hierarchical clustering tab was used for hierarchical clustering of the data in the present study. Hierarchical clustering is a powerful and useful method used to analyze large genomic datasets. Cluster analysis performs four types of binary, agglomerative and hierarchical clustering. The basic aim is to assemble a set of items (genes or arrays) into a tree, where similar items are joined through short branches and increasingly longer branches as the similarity between items decreases. The first step in hierarchical clustering is to calculate the distance matrix between the gene expression data. Once this matrix of distances is computed, clustering is initiated. Agglomerative hierarchical processing includes repeated cycles, where the two closest remaining items (those with the smallest distance) are joined by a node/branch of a tree, with the length of the branch set to the distance between the joined items. The two joined items are removed from a list of items processed and replaced by an item that represents the new branch. The distances between this new item and all other remaining items are computed, and the process is repeated until only a single item remains.

GO analysis

GO analysis was applied to analyze the main function of the differentially expressed genes according to GO, which is the key functional classification of the National Centers for Biotechnology Information that organizes genes into hierarchical categories and can be used to identify the gene regulatory network on the basis of biological process and molecular function (14,15).

Specifically, two-side Fisher's exact test and χ2 test were used to classify the GO category, and the FDR (16) was calculated to correct the P-value; the smaller the FDR, the smaller the error in judging the P-value. The FDR was defined as FDR = 1 − Nk / T, where Nk refers to the number of Fisher's test P-values less than χ2 test P-values. P-values were computed for the GO terms of all the differentially expressed genes. Enrichment provides a measure of the significance of the function; as the enrichment increases, the corresponding function is more specific, which helps to identify the GO terms with more concrete functions in the experiment. Within the significant category, the enrichment Re was given by: Re = (nf / n) / (Nf / N) where 'nf' is the number of flagged genes within the particular category, 'n' is the total number of genes within the same category, 'Nf' is the number of flagged genes in the entire microarray, and 'N' is the total number of genes in the microarray (17).

Pathway analysis

Pathway analysis was used to identify the significant pathways of the differentially expressed genes according to Kyoto Encyclopedia of Genes and Genomes [KEGG (18–20)], Biocarta (19,20) and Reactome (20). Fisher's exact test and the χ2 test were used to select the significant pathways, and the threshold of significance was defined according to the P-value and FDR. The enrichment Re was calculated as aforementioned (18–20).

STC analysis

According to the RVM corrective analysis, differentially expressed genes were selected at a logical sequence. In accordance with the different signal density alteration tendencies of genes under different situations, a set of unique model expression tendencies were identified. The raw expression values were converted into a log2 ratio. Using a strategy for clustering short time-series gene expression data, some unique profiles were defined. The expression model profiles are associated with the actual or expected number of genes assigned to each model profile. Significant profiles have higher probability than expected according to Fisher's exact test and multiple comparison tests (21,22).

Signal-net analysis

The gene-gene interaction network was constructed based on the data for differentially expressed genes. Using Java (23), which enables users to build and analyze molecular networks, the network maps were constructed. For example, if there was confirmatory evidence that two genes interacted, then an interaction edge was assigned between the two genes. The considered evidence included the source of the interaction database from KEGG. The networks were stored and presented as graphs, where the nodes primarily represent genes (protein, compound, etc.) and the edges represent relation types between the nodes, e.g., activation or phosphorylation. The networks were examined using the powerful tools implemented in R (24).

To investigate the global network, the most important nodes were identified. To this end, connectivity, also known as degree, was defined as the sum of connection strengths between the network genes: where Ki is the sum of connection strengths between the network genes; aui represents the weight relations between gene i and gene j.

In gene networks, connectivity measures how correlated a gene is with the other network genes. For a gene in the network, the number of source genes of a gene is known as the in-degree and the number of target genes of a gene is known as the out-degree. The character of genes is described using 'betweenness' and centrality measures, which reflect the importance of a node in a graph relative to other nodes. For a graph G:(V, E) with n vertices, the relative 'betweenness' centrality C′B(V) is defined according to the following equation: where σst is the number of shortest paths from s to t, and σst(V) is the number of shortest paths from s to t that pass through a vertex v (24–28).

Results

Global gene expression analysis of 24-h-, 7-day- and 10-week-old mouse hearts

Affymetrix Mouse Gene 1.0 ST microarray chips were used to generate the transcriptome profiles of 24-h-, 7-day- and 10-week-old mouse hearts, and 35,557 transcripts were presented in the global gene expression profiles. To reduce selection bias, the RVM was used to modify the F-test to calculate the level of significant difference (P-value), and the FDR was used to obtain the differentially expressed genes. A total of 2,736 differentially expressed genes were selected according to the threshold of P<0.05, FDR<0.05, and fold-change ≥1.5 in 24-h-, 7-day- and 10-week-old mouse hearts. A comparison between 24-h- and 7-day-old mouse hearts demonstrated a larger proportion of upregulated genes (1,445 upregulated genes vs. 1,291 downregulated genes), whereas the other pairwise comparisons of the groups (24-h- vs. 10-week-old hearts and 7-day- vs. 10-week-old hearts) revealed that the downregulated genes accounted for the majority of expression alterations (Fig. 1A). Using the hierarchical clustering tab, a clear separation between 24-h-, 7-day- and 10-week-old mouse heart samples was observed among the 2,736 differentially expressed transcripts (Fig. 1B). These results suggested that the differences between the three groups were predominant, and not strongly affected by biological variations and differences within the groups. Therefore, 2,736 differentially expressed genes, accounting for 7.69% of the 35,557 identified transcripts, were selected for further analysis.

Figure 1 Gene expression profiling of hearts obtained from 24-h-, 7-day- and 10-week-old mice. (A) Differentially expressed genes between the three time-points; there were 1,138 upregulated genes and 1,598 downregulated genes in the 24-h group vs. the 10-week group; 1,107 upregulated genes and 1,629 downregulated genes in the 7-day group vs. the 10-week group; and 1,445 upregulated genes and 1,291 downregulated genes in the 24-h group- vs. the 7-day group. (B) Hierarchical clustering of differentially expressed genes. The abscissa represents the sample, and the ordinate represents the different genes. Red indicates high gene expression, whereas green indicates low gene expression.

Significant functional annotation of the differentially expressed genes

To understand the significant functional alterations in the differentially expressed genes, a GO analysis was conducted, which matches the functional annotation to each gene based on cellular component, molecular function and biological process. Subsequently, using Fisher's test and multiple comparison tests, the P-value (significance level) of each functional annotation was obtained. The standard of significant screening was set at P<0.001. A total of 357 significant categories were identified from the GO analysis. The top 10 categories are presented in Fig. 2A. The results demonstrated that with increasing age, the differentially expressed genes were involved in 'cell cycle', 'cell division', 'mitosis', 'metabolic process', 'oxidation-reduction process', 'regulation of transcription, DNA-dependent', and 'protein phosphorylation'. These results indicated that the progressive loss of proliferative potential of cardiomyocytes reflects multifunctional and multisystem co-regulation.

Figure 2 Top 10 significant GO functions (P<6E-29, FDR<3E-26) and pathways (P<10E-21, FDR <2.5E-19). (A) Representative GO annotations. With increasing age, the following enriched GO categories were identified: 'cell cycle', 'cell division', 'mitosis', 'metabolic process', 'oxidation-reduction process', 'transport', 'cell adhesion', 'phosphorylation', 'regulation of transcription, DNA-dependent' and 'protein phosphorylation'. (B) Representative pathways. With increasing age, the following enriched pathway categories were identified; 'systemic lupus erythematosus', 'metabolic pathways', 'cell cycle', 'alcoholism', 'PI3K-AKT signaling pathway', 'viral carcinogenesis', 'phagosome', 'pathways in cancer', 'focal adhesion' and 'ECM-receptor interaction'. The -lgp is displayed on the x-axis. The greater the -lgp, the higher the significance. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Signal transduction pathways associated with the differentially expressed genes

Using KEGG, Biocarta and Reactome databases, the pathways were identified according to systematic relationships, functions and genomic information of the differentially expressed genes. Fisher's test and the χ2 test were used to calculate the P-values, in order to identify the significant signal transduction pathways. A total of 161 significant pathways were identified in the pathway analysis (P<0.05). The 10 most significant pathways are listed in Fig. 2B, including classical 'cell cycle', and proliferation pathways, including 'PI3K-AKT signaling pathway'. 'metabolic pathways' and 'alcoholism'. In addition, pathways associated with the immune system were significantly enriched, which included 'systemic lupus erythematosus', 'viral carcinogenesis', 'phagosome' and 'focal adhesion'. In addition, the 'ECM-receptor interaction' pathway, which is associated with structural and biochemical support for cardiomyocytes, was also enriched. These findings indicated that postnatal cardiac development is a complex regulatory network associated with numerous proliferative, metabolic and immune-associated pathways.

Analysis of gene expression tendency

STC was used to analyze the differentially expressed genes, in order to accurately and intuitively visualize the genes with the same expression tendency between 24-h-, 7-day- and 10-week-old mouse hearts. Using a significance threshold of P<0.05, differentially expressed genes were grouped into five significant patterns (Fig. 3). In profile 1, which contained 169 genes, gene expression was significantly increased in 7 day hearts compared with in 24 h hearts, and was markedly decreased compared with in 10 week hearts (Fig. 3A). In profile 2, which contained 898 genes, there was no obvious difference between the 24 h and 7 day hearts; however a significant decrease in expression was observed in 10 week hearts (Fig. 3B). Profile 3, which contained 254 genes, exhibited continuously decreasing expression from 24 h to 7 day to 10 week hearts (Fig. 3C). Conversely, profile 4, which included 280 genes, exhibited continuously increasing expression from 24 h to 7 day to 10 week hearts (Fig. 3D). In profile 5, which contained 625 genes, there was no significant difference in gene expression between the 24 h and 7 day hearts, whereas much higher expression was observed in 10 week hearts (Fig. 3E).

Figure 3 Five significant trends summarized by the model profile. Each profile indicates an expression pattern. The horizontal axis represents time-points, and the vertical axis indicates the time series of gene expression levels for the gene after log-normalized transformation. The value in brackets after the profile represents the variation intensity, the genes assigned represent the gene number in each profile, and the genes expected represent the theoretical gene number in each profile. The P-value represents the significance level between the actual number of genes and the random distribution of genes in theory. (A) Profile 1 shows the contra-regulation pattern. (B) Profile 2 shows the downregulation pattern after 7 days. (C) Profile 3 shows the downregulation pattern. (D) Profile 4 shows the upregulation pattern. (E) Profile 5 shows the upregulation pattern after 7 days.

These five significant patterns indicated that the genes were involved in dynamic activation in 24-h-old mouse hearts, which contain proliferative cardiomyocytes; 7-day-old mouse hearts, which contain cardiomyoctes undergoing a proliferative burst; and 10-week-old mouse hearts, which contain growth-arrested cardiomyocytes. The dynamically regulated genes included the cyclin genes [downregulated, cyclin (CCN) A2, CCNB1, CCNB2, CCND3 and CCNE1; upregulated, CCNG1], the genes encoding cyclin-dependent kinases (CDKs) (downregulated, CDK1, CDK2, CDK4 and CDK6; upregulated, CDK18), downregulated cell cycle-related genes [cell division cycle (CDC)20, CDC25A and CDC25B], and downregulated transcription factor genes [MYC proto-oncogene, bHLH transcription factor (MYC) and E2F transcription factor 1 (E2F1)]. In addition, growth factor genes [downregulated bone morphogenetic protein (BMP)1, BMPR1a, BMP2, BMP5, BMP7, BMP10, IGF1, IGF2 and IGF2R; upregulated BMP6, fibroblast growth factor (FGF)1 and FGF2], downregulated mitotic checkpoint serine/threonine-protein kinase genes [BUB1 mitotic checkpoint serine/threonine kinase (BUB1) and BUB1B] and downregulated CDK-interacting protein/kinase inhibitory protein genes [CDK inhibitor (CDKN)1C, CDKN2D and CDKN3] were also well enriched. The gene encoding checkpoint kinase 1 (CHEK1), which interacts with numerous downstream effectors to induce cell cycle arrest, was also downregulated in 10-week-old group. Notably, the expression of BMP1, BMP10, CCNE2, E2F1 and IGF1, which were involved in profile 1, peaked in hearts obtained from 7-day-old mice, which was consistent with the findings of Naqvi et al (6).

Signal-net analysis and cell cycle signal transduction network

Signal-net analysis, based on the KEGG database of interactions between different gene products and the theory of network biology, was established to illustrate the inter-gene signaling between the differentially expressed genes. To specifically assess gene interaction, 38 genes were used to construct the cell cycle signal transduction network (0.11% of all tested genes; the genes selected were the genes involved in the annotations and pathways related to 'cell cycle', selected based on the GO and pathway analyses) (Fig. 4). The majority of these genes were involved in profile 2, including classical cell cycle genes, such as cyclin genes (CCNA2, CCNB1, CCNB2 and CCNE1), genes encoding CDKs (CDK1, CDK2, CDK4 and CDK6), CDKN2D, genes encoding the cyclin-dependent kinase regulatory subunit [CDC28 protein kinase regulatory subunit (CKS)1b and CKS2], and cell cycle-associated genes (CDC20, CDC25A and CDC25B). The genes encoding cell cycle checkpoints, regulators and growth factors, including BUB1, BUB1B, CHEK1, anaphase promoting complex subunit 1 (ANAPC1), mitotic arrest deficient 2 like 1, polo like kinase 1 (PLK1), RB transcriptional corepressor like 1 (p107), transformation related protein 53, shugosin 1, ribosomal protein S6 kinase A3, transforming growth factor-β2 (TGF-β2) and thrombospondin 1, were also enriched in this network. The genes involved in profile 5, including growth arrest and DNA-damage-inducible protein 45α (GADD45α) and pituitary tumor-transforming 1 (PTTG1; APC substrate); profile 3, including CDC25C, CCND3, CDKN1C, WEE1 G2 checkpoint kinase (Cdk1 inhibitor), MYC (transcription factor), and host cell factor C1 (transcription factor); profile 4 (CCNG1); and profile 1 (CCNE2) were also well enriched. Notably, in the network, the genes with higher degree (interactive ability), including ANAPC1, CDC20, CDK1, MYC and CDC25C, may serve important roles in myocardial cell cycle arrest from 7 days to 10 weeks. In the network, ANAPC1, which inhibits CCNB1, CCNB2, CDK1 and PTTG1 through ubiquitination, was activated through CDC20. CDC20 also inhibited CCNB1 and CCNB2 through ubiquitination. CDC25C achieved its function through the dephosphorylation of CDK1. These data indicated that the regulation of these five genes (ANAPC1, CDC20, CDK1, MYC and CDC25C) may induce cardiomyocytes to maintain proliferative capacity after preadolescence. Furthermore, in the network, GADD45α is likely to inhibit the expression of CCNB1 and CCNB2, and CDK1 may regulate myocardial cell cycle at birth.

Figure 4 Cell cycle signal transduction pathway and hierarchical clustering of associated genes. Genes associated with the cell cycle were identified through signal-net construction. Cycle nodes represent the genes, and the edges between two nodes represent interactions between genes (arrowheads represent targets). The more edges of a gene indicates that more genes are connected with it, thus suggesting it has a more central role within the network.

RT-qPCR verification of gene expression

RT-qPCR was used to examine nine representative genes, including BMP1, CCNE2 and IGF1 from profile 1, CDC20 and PLK1 from profile 2, CDC25C and MYC from profile 3, CCNG1 from profile 4 and FGF1 from profile 5, to validate the microarray data. As presented in Fig. 5, these results were consistent with the normalized microarray data, based on either fold changes in expression and/or potential functions.

Figure 5 Nine representative differentially expressed genes validated using RT-qPCR. The relative expression levels of the nine genes detected using RT-qPCR (red box) were consistent with the results of the normalized microarray (green box) based on Pearson's correlation analysis (R>0.9, P<0.05). The data are expressed as the mean ± standard error of the mean of three independent experiments. BMP1, bone morphogenetic protein 1; CCNE2, cyclin E2; CCNG1, cyclin G1; CDC20, cell division cycle 20; CDC25C, cell division cycle 25C; FGF1, fibroblast growth factor 1; IGF1, insulin like growth factor 1; MYC, MYC proto-oncogene, bHLH transcription factor; PLK1, polo like kinase 1; R, Pearson's correlation coefficient; RT-qPCR, reverse transcription-quantitative polymerase chain reaction.

Discussion

Mammalian cardiomyocytes retain relatively high proliferative potential at birth, which gradually diminishes after a few days as a result of myocardial cell cycle arrest (29–31). Furthermore, during adolescence, a burst of proliferation occurs in cardiomyocytes (6). However, the alterations in gene expression during the critical period from birth to adolescence and then to adulthood remain unknown. To elucidate the molecular mechanism and the gene expression alterations that occur during this key period, the present study conducted global gene expression profiling to generate transcriptomes from 24-h-, 7-day- and 10-week-old C57BL/6 mouse hearts. These time-points reflect the phases of proliferative activation, proliferative burst and proliferative arrest, respectively. In addition, GO annotation, STC, pathway enrichment and network construction were conducted to identify the state of cardiomyocytes from numerous perspectives.

As expected, the GO analysis revealed that, with increasing age, the expression of genes involved in the cell cycle, cell division and mitosis, as well as metabolic process, protein synthesis and protein modification was significantly altered. These results were highly consistent with the findings of the pathway analysis, which indicated that the phosphoinositide 3-kinase (PI3K)-protein kinase B (AKT) signaling pathway was enriched. These findings indicated that the cell cycle and protein modification are important in postnatal mouse heart development. Numerous proliferative stimulators directly activate the PI3K-AKT pathway by promoting nuclear trans-location, and increasing the phosphorylation and degradation of proteins associated with myocardial proliferation (32). The progressive loss of proliferative potential in cardiomyocytes reflects multifunctional and multisystem co-regulation. In addition, pathways associated with the immune system were enriched, reflecting the fact that neonatal mice are undergoing a transition from reliance on the intrauterine physiological environment to self-sufficiency. Following birth, the immune system faces the challenge of transferring from a sterile environment to a microbe-containing environment (33).

Notably, in the present study the differentially expressed genes were grouped into five patterns, according to the relative expression of the differentially expressed genes at three different time-points. Specifically, the majority of cyclin genes (CCNA2, CCNB1, CCNB2 and CCNE1), CDK genes (CDK1, CDK2, CDK4 and CDK6) and cell cycle-related genes (CDC20, CDC25A and CDC25B) were involved in profile 2. In addition, cyclin genes (CCNE2) and growth factors, including BMP1, BMP10, E2F1 and IGF1, were involved in profile 1, which exhibited a burst of increased expression in 7-day-old hearts. These genes may increase mice myocardial proliferative potential within 7 days of age and could potentially elicit myocardial proliferation in damaged adult hearts. These data were consistent with the findings of Naqvi et al, which generated a molecular model of myocardial proliferation involving thyroid hormones and the IGF1/AKT pathway (6). In addition, Rochais et al reported that IGF1 regulated activation of the PI3K/AKT pathway to downregulate the expression of p21 and p27 in proliferating cardiomyocytes (34). BMP1 is essential for normal heart development; also this gene serves non-redundant roles in the genetic and morphogenetic processes in human congenital heart disease (35). Furthermore, BMP10 is essential for maintaining cardiac growth during murine cardiogenesis (36). Although E2Fs promote S-phase entry and DNA synthesis, Ebelt et al reported that E2F1 promotes mitotic cell division by increasing the expression of cyclin B1 and B2 in neonatal cardiomyocytes (37). Therefore, these genes are likely to be involved in the burst of cardiomyocyte proliferation during early preadolescence. Theoretically, stimulating the expression of these genes may increase proliferative capacity in post-adolescent and adult hearts.

To date, some critical genes, including BUB1, dynein light chain Tctex-type 1A (DYNLT1α), tropomodulin 4 (TMOD4) and glutamate metabotropic receptor 1, which were screened in the present study, have not been described in the context of myocardial proliferation. However, in other tissues, these genes serve regulatory roles in cell proliferation. Bieniek et al reported that the marked downregulation of key proteins in kinetochore/centromere assembly, including Cdc20, Ndc80, Bub1 and Plk1, may arrest prostate cancer proliferation (38). Yeh et al reported that DYNLT1α encodes Tctex-1, which is involved in the IGF1R-Gbg-phospho (T94) Tctex-1 signaling pathway and promotes the proliferation of neural progenitors via the modulation of ciliary resorption and G1 length (39). Berger et al reported that loss of the Tmod4 protein may cause cytoplasmic rod formation and muscle weakness reminiscent of nemaline myopathy in zebrafish (40). Therefore, it may be hypothesized that these genes serve important roles in myocardial proliferation.

Notably, genes associated with progenitor cells, pluripotent stem cells and fibroblasts were also well enriched in the present study. A previous study confirmed that cardiac progenitor cells reside in postnatal hearts (41). These genes were identified using the expression of surface markers, such as c-Kit (42,43), which was encoded by KIT and was downregulated in the present study. BMP genes (BMP1, BMP5, BMP7, BMP10 and BMPR1a), which are required for the optimization of efficient cardiac differentiation from pluripotent stem cells, were also downregulated in the present study (44). Fibroblasts and myofibroblasts secrete collagens (encoded by COL genes), which function as primary structural proteins and provide mechanical support (45). Collagens can be degraded by matrix metalloproteinases (MMPs), which are a group of endopeptidases, thus enabling cell migration into the injured area (46). MMPs are inhibited via the tissue inhibitors of metalloproteinases (TIMP) family, which is associated with extracellular matrix (ECM) expansion and cardiac fibrosis (47). In the present study, COL genes, MMP genes (MMP2, MMP14 and MMP16) and TIMP1 were downregulated. The homeostasis of ECM factors, which exhibited a similar expression trend between 24-h- and 10-week-old hearts, provides a perfect scaffold for all cardiomyocytes. In addition, the 'ECM-receptor interaction' pathway was enriched in the pathway analysis. Other factors associated with cardiac fibrosis were also well enriched; TGFβ-2, FGF12, FGFR2 and platelet derived growth factor (PDGF)A were continuously decreased after birth, whereas FGF1, FGF2, FGF7, FGF16 and PDGFD were highly expressed in 10-week-old hearts. Considering the fact that at birth, the heart can completely repair its site of injury without scar, whereas fibrotic scar formation occurs in injured adult hearts, the results of the present study may provide a basis for the physiological condition of progenitor cells, pluripotent stem cell differentiation and fibrotic scar formation.

A signal transduction network analysis of the cell cycle was used to efficiently investigate the phenomenon of myocardial cell cycle arrest. The results suggested that ANAPC1, CDC20, CDK1, MYC, CDC25C and GADD45α, which exhibited the CDK1 and PTTG1 through ubiquitination, and was activated by CDC20, which also inhibited the expression CCNB1 and CCNB2 via ubiquitination. A previous study reported that APC/cyclosome (C) regulates the ubiquitin-dependent proteolysis of specific cell-cycle proteins to coordinate chromosome segregation in mitosis and entry into the G1 phase (48). The catalytic activity of the APC/C and its ability to initiate the destruction of particular proteins at different phases of the cell cycle are controlled through interactions with the structurally related coactivator subunit, Cdc20. In cardiomyocytes, blocking APC/CCdc20 may lead to G2 arrest (49). These findings are also consistent with the data from the GO and pathway analyses. ANAPC1 and CDC20 shared the same functional annotation, which included cell cycle, cell division and mitosis, and both participate in cell cycle and ubiquitin-mediated proteolysis pathways. These genes were involved in profile 2, which exhibited downregulation from 7 days to 10 weeks of age. Together, these results suggested that ANAPC1 and CDC20 may have important roles in the phenomenon of myocardial cell cycle arrest from early puberty to adulthood.

In the network, CDC25C was activated through the dephosphorylation of CDK1, which activates CDC25C by phosphorylation. Consistent with a previous study by Forester et al Cdk1 activity is controlled through the phosphatase Cdc25C, which is turned on at the G2/M transition to catalyze Cdk1 activation (50). This constitutive activation of Cdc25C and Cdk1 leads to a delayed exit from mitosis. Franckhauser et al observed that during mitosis, the nonphosphorylated mutant forms of CDC25C impair mitotic progression in human cells (51). In addition, cell cycle, cell division, mitosis and protein phosphorylation were well enriched in CDK1, and CDK1 and CDC25C were involved in the cell cycle pathway. MYC activates CDK4 and CDK6 in the network. Villa Del Campo et al reported that overexpression of MYC increases the population of cardiomyocytes and enhances the epicardial contribution to the developing heart (52). In addition, GO analysis revealed that MYC has a large role in the regulation of transcription and participates in the cell cycle, and PI3K-Akt, Hippo, mitogen-activated protein kinase (MAPK) and TGF-β signaling pathways in the pathway analysis. These results suggested that MYC may be involved in various functions, and function as an intermediate mediator of various pathways in myocardial cell cycle arrest. MYC and CDC25C were both involved in profile 3, which began to decrease in 24-h-old hearts and may initiate regulatory functions in the cell cycle at birth.

In the network, GADD45α inhibited the expression of CCNB1, CCNB2 and CDK1. GADD45α serves a role in S-phase and G2/M arrest (53). Similarly, Gadd45α binds Cdk1, likely preventing its association with cyclin B1, inhibiting Cdk1 activity, and arresting the cell at the G2/M checkpoint (54). At present, no report has been published regarding the relationship between myocardial cell cycle arrest and GADD45α. However, the results of an in vitro assay of cortical neuron development indicated that overexpression and knockdown of GADD45α transcription suppressed the formation of distal neurite processes and often promoted aberrantly shaped and sized cell bodies. In addition, overexpression did not affect migration but rather led to irregular and hypertrophied cell body development (55). The results of the pathway analysis in the present study indicated that GADD45α may be involved in the p53 and MAPK signaling pathways, and was the only gene from profile 1 with increased expression at 7 days old. These results suggested that the increased expression of GADD45α may serve important roles in myocardial cell cycle arrest from early puberty to adulthood.

Gene expression profile analyses have the advantage of generating huge amounts of biological information, which can be objectively filtered and accurately analyzed. The present study used GO annotation, pathway analysis and STC to analyze the specific functions, pathways and logical patterns of the differentially expressed genes. Notably, a potential cell cycle signal transduction network was constructed in the present study. This network enables the screening of the potential genes associated with a single function, the cell cycle. The present study may have value as an important reference concerning the gene expression profile in the heart during proliferation.

Acknowledgments

The present study was supported by the National Natural Science Foundation of China [grant no. 81570256 (to M.W.)], and the China Scholarship Council [grant no. 201506940022 (to R.W.)]. The authors would like to thank Dongjia Ge for contributing to the design of the figures. The global gene expression data are available at the Gene Expression Omnibus website (https://www.ncbi.nlm.nih.gov/geo/) under accession no. GSE93426.

References