Keratin5-BMP4 mechanosignaling network regulates cell cyto- skeletal reorganization and chromatin accessibility during re- vascularization-based blastema regeneration

Heart regeneration after myocardial infarction remains challenging in reconstruction of blood resupply system. Here, we find that in zebrafish heart after resection of the ventricular apex, the local myocardial cells and the clotted blood cells undergo cell remodeling process via cytoplasmic exocytosis and nuclear reorganization within revascularization-based blastema. The regenerative processes are visualized by spatiotemporal expression of three blastema representative factors (alpha-SMAwhich marks for fibrogenesis, Flk1for angiogenesis/hematopoiesis, and Pax3a for remusculogensis),and two histone modifications (H3K9Ac and H3K9Me3 mark for chromatin remodeling). Using the cultured zebrafish embryonic fibroblasts we identify blastema fraction components and show that Krt5 peptide could link cytoskeleton network and BMP4 signaling pathway to regulate the transcription and chromatin accessibility at the blastema representative genes and bmp4 genes. Our study provides new mechanistic insights into the epithelial-dependent and revascularization-based blastema regeneration for potential myocardial infarction therapy.


Introduction
Coronary artery disease (CAD) also known as ischemic heart disease is the most common type of cardiovascular disease. Limitation of blood flow to the heart results in cellular oxygen starvation (ischemia), and cardiomyocyte necrosis (myocardial infarction，MI). The death of cells then triggers a cascade of local inflammatory response and fibrogenesis [1,2]. If the impaired blood flow does not deteriorate, myofibroblasts and endothelial cells proliferate and migrate into the infarct zone, where myofibroblasts deposit a network of collagen to form a hypocellular collagenous scar [3,4]. The resultant myocardial scarring could impede conduction velocity and cause lethal arrhythmias. Besides scarring, another common problem is the second round of ischemia-reperfusion injury (IRI) or reoxygenation injury after a period of ischemia or lack of oxygen and nutrients from blood [5]. Therefore optimal myocardial infarction therapy requires a scarless regeneration, which can modulate inflammatory reaction, blood resupply and myocardial renewal.
Gene edition of periostin-expressing myofibroblasts could reduce collagen production and scar formation [6], and forced viral-mediated gene expressions or chemical cocktails could induce reprogramming of mouse non-muscle cells into cardiomyocytes [7,8]. Intriguingly, inflammatory cytokines, such as interleukin-1beta (IL-1beta) and tumor necrosis factor-alpha (TNF-alpha), could induce the transformation of human dermal microvascular endothelial cells into myofibroblasts in skin fibrogenesis [9,10], and the (myo)fibroblast phenotype could reverse back to the precursor phenotype [11,12]. These data elicit anticipation to convert the inflammatory reaction and inflammatory cells into cardiac tissue regeneration [13].
Several lines of evidence support a notion that zebrafish can avert myofibroblast fibrogenesis into blastema regeneration [14][15][16] . Historically, both myofibroblast and blastemal cells are characterized by fibroblast phenotype with elevated expression of extracellular matrix. They are heterogeneous in the cell origin and plastic in the phenotype. Myofibroblasts can arise from conversion of resident fibroblasts, differentiation of bone marrow-derived progenitors, endothelial-mesenchymal transition (EndMT) and epithelialmesenchymal transition (EMT) [11,17]. Likewise heterogeneous blastema cells in lower jaw regeneration originate from fibroblast cells, nucleated blood cells, fragmented muscle cells, pigment cells, and local vascular system [16,18]. Different from fibrogenesis with the accumulation of scar tissue, blastema proceeds with tissue respecification and blood vessel reconstruction. Similar cellular response and regenerative capacity of the zebrafish heart has been reported after cryoinjury, a procedure that closely models the pathophysiological process undergone by the human heart after MI [19]. Thus, the zebrafish regeneration strategy could provide an excellent model for MI therapy.
Here, we employed transgenic reporter fish and cell identity analysis, traced the fate of the injured myocardial cells and blood clots, and investigated molecular mechanisms by which heart blastema regenerate in zebrafish. Our findings at subcellular and molecular levels revealed a dual regulation of cytoplasmic intermediate filament mechanosignaling and histone modification switch on cell cytoplasmic remodeling and nuclear reorganization during heart blastemal regeneration.

Results
This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

Epicardium and subepicardium derived cells invaded into the resection site and induced fibrinolysis and blastema formation
The previous studies showed that zebrafish can vigorously regenerate the amputated cardiac muscle and restrict scar formation [15,20]. However, these studies did not clarify how zebrafish renovate the thrombus and reconstruct the blood supply system. In this study, we surgically removed up to 20% of the ventricle apex and then studied the progresses of bleeding, blood clot formation, and wound healing during zebrafish heart regeneration ( Fig S1). As proliferation and migration of pericardial cells sealed the wound and the blood clot formed, the resultant pericardial-covered clot mix (hereafter called epiclot) was established a pre-blastema field, which morphologically covering the resected apex (Fig S1A and 1B). From 2 to 3 days post-amputation (dpa), the adjacent subcardial cells invaded into the epi-clot, and triggered tissue lyses, fibrogenesis, blastema revascularization (new vessels invasion) and cell remodeling (Fig 1B-C and S1B). Meanwhile, new types of light staining cells (or pre-blastema cells) emerged from the lysed tissues and clot, and then transformed to heterogeneous blastemal cells (Fig 1B-C and F). By 8dpa, the preblastema field was filled with a variety of blastemal cells as well as the residual epi-clot and hypocellular mus-clot (Fig 1G-H and S1D). During blastema reformation/respecification from 10dpa to 30dpa when the myofibrils remuscularized the blastema and encapsulated the residual clot, the inside nucleated blood cells (NBCs) were induced and converted into various types of hematopoietic-derived cells in morphology (Fig 1I and S1E-G). Here we not only confirmed the previous observations that epicardial and subepicardial cells were actively involved in blastema formation [20], but also showed that blastema formation was initiated from fibrinolysis and musculolysis, followed by fibrogenesis and revascularization/ hematopoiesis.

Figure 1.
Histological analysis of blastema formation and reformation during heart regeneration. After amputation of the ventricle tip, initial response was bleeding and hemostasis. In the following days, the epicardial cells (A) immigrated and incorporated with the outmost fibrined nucleated blood cells (NBCs) to form a pericardium-like outer layer covering the NBCs clot (B). Tissue lyses began along the intersection between the regenerating epicardium-clot interaction (epi-clot int), subcardium-clot interaction (sub-clot int), and the damaged muscle-clot interaction (mus-clot int) (B-E). The light cells refer to those cells that showed hematoxylin-resistant and had a large nucleus with vacuolar cytoplasm. As more subepicardial cells invaded, the light-staining areas progressively expanded toward the mus-clot int (D) and epi-clot int (E). Subsequently, the light cells became basophilic and stem cell-like in morphology (F). With the increase of blood vessels, a typical blastema was formed by the mass of basophilic and fibroblast-like cells with an excess of the extracellular matrix (G and H). During the weeks that follow a ventricular resection, the deposited extracellular matrix and blastema were gradually replaced by regenerated myocardium meanwhile NBCs within the encapsulated residual clots (RCs) were progressively induced (I). As shown in the islets, there were two size types of light cells. The large light cells (B, C and F) had abundant cytoplasm and a variety of vacuoles; the small light cell (D, E and I) contained a small amount of cytoplasm and a relatively large intact nucleus. ECs, epicardial cells; Mf, myofibril; BV, blood vessel (encircled by dashed line); PHSCs, presumptive hematopoietic stem cells; RCs, redisual clots. Scale bar = 200 μm.

Blastema cells displayed cytoplasmic exocytosis and nuclear reorganization
Compared to the histological sections, the transformation processes of the local myocardial cells and clotted blood cells were more discernible at the subcellular level by transmission electron microscopy (TEM). TEM analysis revealed that the transforming myocardial cells displayed cytoplasmic vacuolation and efflux, and reproduction of new cytoplasmic components (Fig 2A). In contrast to the extensive cytoplasmic autolysis, the nuclei were resistant to degradation (Fig 2B-D). As seen within the activated NBCs, nuclear reorganization started from nuclear condensation and heterochromatin formation (Fig 2E-I). Some of the induced blood clot cells exhibited morphological characteristics similar to the presumptive hematopoietic stem cells (PHSCs) described by Rubinstein [21] ( Fig 1F, 2E-I). Together, these histological and ultrastructural observations clearly showed the blastema cell resalvage processes from the local cardiomyocytes and clotted NBCs through cytoplasmic exocytosis and nuclear reorganization.

Figure 2.
Ultrastructural observations on the cytoplasm remodeling and nuclear reorganization of blastema cells by transmission electron microscopy (TEM).During the blastema formation, the large light cells contained numerous transparent vesicles (A). These membrane-bound vesicles were limited by a double or single membrane and usually filled with granules or other inclusions. Some vesicles were fused into bigger autophagosome-like vacuoles, or eliminated from the cells (B and C). The chromatin lining the nuclear periphery throughout the nucleus became condensed, and the nuclear pore-like structure became evident due to the complicated folding of the nuclear envelope (D-G). Endoplasmic reticula and Golgi apparatus aggregated around the poles of the nucleus. Mitochondria and myofibrils were scattered throughout the myocyte organelles (G). Within the NBCs, the nucleus was initially condensed, and left an empty rim around the nucleus (D and E). Then the nucleus was enlarged with increasing thick clumps of chromatin beneath the nuclear membrane (G-I). N, nucleus; Mf, myofibril; Gg, Golgi; V, vesicle, AV, autophagosome-like vacuole. Scale bars: A = 0.5 μm; D and G = 1 μm; the rest = 2 μm.

Three tissue-specific pioneer factors tracked myofibroblast and blastema formation
According to the previous report, zebrafish adult cardiac muscle regeneration is comparable to embryo cardiac morphogenesis from the primitive embryonic structure into three mature myocardial lineages of primordial, trabecular and cortical cardiomyocytes [22]. To trace the origin and phenotypic respecification processes of blastema cells, we utilized isl1:GFP and flk1:GFP transgenic fish to visualize isl1+-and flk1+-GFP fluorescence in the injury site, and also immunostained for myofibroblast marker alpha-smooth muscle actin (alpha-SMA) and Pax3a protein to visualize the proliferating myofibroblasts and myogenic progenitors, respectively. Except signal of isl1-GFP undetectable (data not shown), the spatiotemporal expression patterns of alpha-SMA, Flk1 and Pax3a were roughly paralleled to blastemal fibrogenesis, angiogenesis and myogenesis, respectively ( Fig 3A). These results indicated that the three proteins could mark myofibroblast formation, blood vessel reconstruction and muscular regeneration.
To unveil the genotype-phenotype correlation during the blastema regeneration, we conducted time-point based RNA-seq analyses and identified 626 significant transcripts between three time points (P < 0.05). Functionally, these genes reflected the blastema re-generative progress as above histomorphological descriptions, and exhibited the equivalent functions of kdrl (flk1), acta2, and pax3a in contribution to the extracellular matrix organization, cell migration, blood vessel morphogenesis and muscle tissue development ( Fig 3B, Table S1). In accordance with progressive activation of kdrl, transient activation of acta2 and subsequent pax3a activations, their dynamic expression patterns shown in Fig 3A and Fig 3C, suggested that the three blastema mark genes are temporally expressed in three tissue-specific progenitors along with the blastema regeneration progress. , only a few Flk1-GFP cells and alpha-SMA-expressing cells emerged from the damaged muscle and the subepicardium (a and g), then extensively penetrated to blastemal center and the fibrin clot, and expanded to the whole pericardial field (b and h). By 9dpa, sparse alpha-SMA expression was detected in the outskirt of the blastema, presumptively in regenerating subepicardium (i). Pax3a expression was consistently increased in the regenerating epicardium and outskirt of the blastema, paralleling with masculaturation during regeneration (d-f). Phalloidin-labeled F-actin indicated cardiac muscle staining (jl). Scale bar = 50 μm. (B) Bubble plot showing the GO (gene ontology) enrichment analyses of 626 significantly differential expression genes. Key GO terms were colored in red. Gene ratio refers to the ratio of the number of genes to the total number of genes enriched in the terms. BP, biological processes; CC, cellular component; MF, molecular function. (C) Representative images of H3K9Ac and H3K9Me3 immunohistochemistry. Anti-H3K27Ac antibodies were intensively stained in almost all cells (a and d). In contrast, H3K9Ac and H3K9Me3 were spatiotemporally stained in the regenerative cells (b-c and e-f). At 6dpa, intense H3K9Ac signals (b) mainly aggregated at the outmost epicardium (dashed line) and less in the muscle-clot mix in contrast to the enhanced H3K9Me3 expression (c) in the epi-clot mix. At 9dpa, both H3K9Ac (e) and H3K9Me3 (f) expressing cells increased but differentially distributed over epi-clot mix in blastema region. Scale bar = 50 μm. (D) The WashU Epigenome Browser tracks show the putative H3K9Ac/Me3-specific enhancers and transcription factor (TF) binding sites related to three blastema mark genes using HOMER v4.9. Dashed lines indicate the enriched peaks. Color bubbles indicate the predicated TFs and binding motifs, which were potential blastema formation and respecification regulators during zebrafish heart regeneration.

H3K9Me3 and H3K9Ac switch dynamically marked the blastema cell proliferation and tissue respecification
Considering epigenetic regulation as a bridge between genotype and phenotype, we examined the distributions of 8 common epigenetic marks during heart regeneration using immunohistochemistry ( Fig S2). Among the tested histone modification marks H3K9Ac and H3K9Me3 exceptionally exhibited a dynamic change from the epicardium to the blastemal center. The signal of H3K9Ac is reversely correlated with the signal of H3K9Me3. Specifically, intense H3K9Ac territories along the epicardium showed less staining for H3K9Me3, while dense H3K9Me3-staining foci within the blastema center showed week H3K9Ac signals (Fig 3C). Subnuclear localization showed intense punctual staining of H3K9Ac in the nuclei of mus-clot cells, in contrast to intense H3K9Me3 signal in DAPI-heavy pericentric heterochromatin of the clotted NBCs ( Fig S3). The spatiotemporal staining pattern of H3K9Ac was coincidence with the blastema formation route from the epicardial cells to the inside blastema cells while H3K9Me3 staining patterns accompanied blastema cell expansion from blastema center to outskirt remuscularization ( Fig  3C and S3). When the regenerating heart was pretreated with H3K9Me3-specific and H3K9Ac antagonists, either H3K9Me3 or H3K9Ac reduction impaired blastema formation as indicated by the Flk1-GFP and DAPI staining signals (Fig S4). These results suggest that dynamic changes of H3K9Me3 and H3K9Ac deposition are generally correlated with blastema cell proliferation and respecification.

Genome-wide identification of potential regulatory sequences involved in chromatin state dynamics
The patterns of H3K27Ac and H3K27Me3 profiles along the genome were previously used to identify potential enhancers and promoters in particular cells [23]. In the same way, we performed genome-wide ChIP-seq analyses to identify the potential regulatory elements of the H3K9Ac/Me3 marks ( Figure S5). The distribution of the two chromatinpredicted enhancers was primarily distal to the transcription start sites (TSSs), with a total of 73.36% lying in intergenic regions, and 0.29% falling in at promoter-TSS sites ( Figure  S5C). The intergenic regions enriched 28.03% of H3K9Ac and 79.55% of H3K9Me3, compared with 15.06% of H3K9Ac and 2.57% of the H3K9Me3 poised at promoter-TSS sites ( Fig S5).
Even if the two histone marks were colocalized at the intergenic and promoter-TSS sites (Table S2, Fig 3D), H3K9Me3-dominant sites showed less enrichment for H3K9Ac, whereas those dominantly marked by H3K9Ac were less enriched for H3K9Me3 (Fig 3D  and S6). The reciprocal situation was also confirmed by ChIP-PCR evaluation of the enrichment of the two histone modification marks at the promoter regions of several interest genes, particularly obvious at the pax3a and acta2 genes (Fig S7). Statistical analyses of the top 1000 of the enriched peaks supported reverse correlation between the two histone mark depositions (r=-0.56) but displayed uncertain relationship between the H3K9Ac/Me3 depositions and the targeted gene transcriptions (Fig. S6, Table S3).

Identification of blastema subcellular components responsible for blastema cell nuclear reorganization and cytoplasmic remodeling
Since cell-derived extracellular vesicles (EVs) are powerful conveyors of messengers for intercellular communication [30][31][32], we hypothesized that the vacuolation and exocytosis in the regenerating blastema serve as the vesicle-mediated intercellular communication signals within the blastema. We partitioned the regenerating tissues into three subcellular fractionations, and analyzed the proteins by mass spectrometry. Out of 1561 differential proteins (6dpa vs. 0dpa) identified from three subcellular fractions, 82 proteins were included in the 626 differential transcripts.
We then transferred each fraction into the in vitro cell cultures of the PAC2, and examined the cellular transcription profile. We found that the induced transcription profiles were functionally comparable to the identified 82 proteins and 626 differential transcripts from the regenerating blastema. Table S1 listed the equivalent functional images of the three blastema mark genes and two histone marks.
Further functional clustering analysis of the 82 proteins showed the highest enriched ten extracellular molecules and five intermediate filaments (Table S6). Some of these specific components have been reported to accompany the cell shape change [33], and are implicated in both nuclear and cytoskeletal organization [34,35].

Krt5 targeted cytoskeletal system and non-canonical BMP4 signaling pathway
Since H3K9 methylation is an important epigenetic determinant for chromatin accessibility, and H3K9 methyltransferases are downstream targets of bone morphogenetic proteins (BMPs) [36], and overexpressing BMP-antagonist noggin under the control of the keratin 5 (K5) could reduce apoptosis and cell differentiation in the eyelid epithelium [37], we directly tested the effects of enforced Keratin 5 (Krt5) and noggin on BMPs signaling pathway and on H3K9Ac/Me3-marked chromatin architecture. Because FBS contains abundant BMPs, and noggin is an effective natural BMP antagonist [36], including BMP4, BMP6 and BMP9, we cultured PAC2 cells in an FBS-contained culture medium with or without Krt5 and/or noggin. PAC2 cells did not show detectable Krt5 signals at basal culture conditions. In Krt5-treated PAC2 cells, the Krt5 showed typical coarse dot-like structures scattered through the cytoplasm, particularly displaying the string of beads-like structure surrounding the periphery of the nucleus (Fig 4A, B, Fig S9). In comparison with the untreated control and Noggin treated responses, Krt5 increased F-actin assembly around the nuclei, and decreased the pSmad1/5/8 staining intensity. Krt5 also slightly increased the ratio of H3K9Me3 and H3K9Ac depositions.
Next, we examined whether the enforced Krt5 and noggin impacted the expression of three core blastema genes and two histone modification-related genes. As shown in Fig  4C, Krt5 and noggin upregulated pax3a and acta2 expression, but reversely regulated transcription of bmp4, bmp3 and kdrl. Unlike noggin, Krt5 seldom changed the transcription level of the tested histone modification enzymes. However, Krt5 and noggin mixture enhanced the noggin induction of histone acetyltransferase (esco2), deacetylase (tblxr1a), and demethylase (phf8, kdm7ab) activities.
To clarify the Krt5 transduction network, we took advantage of Krt5 (N-6xHis) and/or noggin (C-Fc) to isolate Krt5-and Noggin-precipitated proteins by specific affinity purifications. The mass spectrometry analyses demonstrated that Krt5 closely interacted with actin-binding protein (Actn1) and non-muscle myosin isoforms (Myh9b, Myh10). The other binding proteins with low affinity contained the heat shock protein 70 (Hspa5, Hspa8), adaptor-related protein complex 3 (Ap3d1), the Golgi apparatus proteins (Golga2), endoplasmic rough endoplasmic reticulum membrane protein (Rpn1), pyrimidine biosynthesis protein (cad), and heterogeneous nuclear ribonucleoprotein (Hnrnpub) ( Table S8). The proteomic profiling of Krt5-binding proteins was matched with the Krt5 distribution map within PAC2 cells. It also indicated that noggin could compete with Krt5 to interact with Myh9b and Myh10.
As BMP signal transduction could be regulated by potent intracellular and extracellular antagonists [38], the differential regulatory effects between Krt5 and Krt5-noggin mixture indicated that Krt5 did not take classic transmembrane receptor pathway to inhibit phosphorylation of Smad 1/5/8, but activated non-canonical BMP4 signaling pathway through intracellular mechanosignal transduction mechanisms (Fig 4E). The hyper-and hypo-accessible regions were shaded with red and blue, respectively. (F) A proposed model that Keratin 5-cytoskeleton mechanosignaling crossly links to BMP4 non-canonical cytoplasmic signal transduction pathways in the regulation of cytoskeletal reorganization and target chromatin accessibility at the loci of blastema cell remodelingrelated genes. Chromatin modifiers were indicated by the H3K9 acetylation/methylation switch, temporal blastema respecification regulator (tBRR), and constitutive blastema respecification regulator (cBRR). The target genes were represented by bmp4 and three blastema mark genes.

Discussion
Cell plasticity studies have proposed a desirable regenerative strategy that relies on the availability of stem cells and/or on the dedifferentiation and/or transdifferentiation of differentiated cells to replenish the missing structure [39]. For heart regeneration after myocardial infarction, various anti-inflammatory agents, and stem cell and gene therapies have been tested to promote myocardial recovery and ameliorate reperfusion injury [5]. However, to date, no conclusively effective strategies have been established. In this study, we found that zebrafish heart ventricle resection also encounters bleeding, hemostasis, and blood clot. Importantly, our data indicate that to achieve regeneration, zebrafish blastema regeneration could integrate inflammatory reaction signals and epithelial-mesenchymal transition (EMT) into cell remodeling networks, thereby salvaging the damaged myocardium and resolving blood clots.
First, zebrafish can faithfully repair the complicated tissues through regeneration blastema [14][15][16]. Several theories have emerged on the origins of the heart blastema cells, such as the proliferation of the reserved undifferentiated cardiac progenitors, or transformation from differentiated myocardial and nonmyocardial cell types [40]. The present study clearly showed that regeneration blastema arose from an advancing cell remodeling, which was marked by three representative cell populations in the blastema (Flk1+ hemoangiogenic progenitor, alpha-SMA+ myofibroblast, and Pax3a+ myogenic progenitor). In essence, the alpha-SMA cells-dominant fibrogenesis provides a temporal tissue scaffold to pave the way for Flk1+ and Pax3a+ cell respecification within the blastema.
A significant observation on blastema cell remodeling was the coexistence of cytoplasmic efflux and nuclear reorganization in the damaged muscle cells. The extensive cytoplasmic degradation is accompanied by vacuole efflux and retention of the condensed nucleus with small cytoplasm. Similar cytoplasmic vacuolation and efflux have once been reported in exocytosis of melanosomes described in cell type conversion of Newt iris epithelial cell type conversion [41]. The cytoplasmic lysis and exocytosis in the blastema seem to be lysosome-dependent, and can get support from the identification of 28 proteolysisrelated genes (Table S1, Fig S10).
Another important observation was that the induced clot compartment showed very similar morphology to the previously described hematon, a multicellular functional unit of hematogenesis in normal mammalian bone marrow [42], which contains the putative hematopoietic stem cells and the stem cell niche [43]. Transcriptomal analyses also revealed that several clusters of genes that are responsible for angiogenesis and hematopoiesis were significantly elevated (Table S1).
Secondly, reverse histone modification patterns (such as methylation of H3K4/H3K27, and H3K27Ac/H3K27Me3) have been reported to regulate chromatin architecture and function in cell lineage specification [23,44,45]. The present results suggested that reciprocal acetylation and methylation of H3K9 could regulate the blastema cell proliferation and respecification through their specific enhancers and enhancer motifs (the key blastema regulators including temporal blastema respecification regulator and constitutive blastema respecification regulator) (Fig 5C). Presumably, those constant motifs were recruited to maintain the cell differentiation status while the temporal-sensitive motifs were responsible for the blastema cell respecification. Our RNA-seq and ATAC-seq analyses further suggested that low-level induced global transcription and chromatin accessibility delimit the blastema cell reprogramming at tissue-specific progenitors rather than pluripotent stem cells.
Finally, the present study can clarify why zebrafish blastema regeneration is epithelial-dependent [16,18], and how they undergo EMT [15]. Our data showed that the epithelial-predominant keratin could form the cytoskeletal network to bridge mechanosignaling between extracellular-cytoplasmic BMP4 signal transductions and histone modifications-mediated chromatin remodeling, thus clarifying how the keratin could modulate IF-cytoskeleton organization, interact with focal adhesions and cell-surface molecules [46,47]. BMP signaling has long been acknowledged to play a pivotal role in mesoderm induction, vascular biology and hematopoietic commitment [48,49]. Consistent with the reports that keratins could regulate yolk sac hematopoiesis and angiogenesis through BMP4 signaling [50], and BMP-Wnt signaling could specify hematopoietic fate through the Cdx-Hox pathway [51], our results further suggested that keratin 5-BMP4 non-canonical pathway could induce mesoderm-derived embryonic fibroblasts undertaking blastema-like angiogenesis and musculature through intracellular mechanosignal transduction mechanisms. Given that the blastema is a mixture of tissues, and PAC2 is a fibroblast cell line, our newly established keratin 5-BMP4 cytoskeleton mechanosignaling network should represent a basic mechanism of mesoderm-derived blastema regeneration, particularly toward blood/vascular system reconstruction.

Animals
Zebrafish TAB lines, Flk1 promoter-derived green fluorescent protein (GFP) expression construct Flk1: GFP transgenic zebrafish (Flk1-transgenic fish) were reared and used for the experiments in accordance with the Animal Care and Use Committee guidelines of Shanghai Ocean University (SHOU-DW-2016-004).

Sample process, histology and immunohistochemistry
After resection of the ventricular apex, the heart tissues were dissected at the day indicated. Histology sectioning, and staining for the immunohistochemistry (IHC) were performed as described previously [18]. For ultrastructural observations, heart tissues were dissected as mentioned above, and fixed 3 hours in Karnovsky's fixative at 4°C. After 6 PBS (phosphate buffer) washes, tissues were again fixed for 2 h in 1% osmium tetraoxide in phosphate buffer at 4°C. After several washes, the tissues were dehydrated in graded acetone solutions and embedded in Epoxy resin using Epoxy embedding medium kit (45359, Sigma, St. Louis, MO, USA). Ultra-thin sections of 60-80 nm thickness were stained with alcoholic uranyl acetate and lead citrate for appropriate time intervals. The grids were then examined with Transmission Electron Microscope HT7700.

RNA-seq library construction and data analysis
Total RNA was isolated from tissues using TRIzol Reagent (Invitrogen) and further purified using RNeasy Mini Kit (Qiagen) and ran on Agilent Bioanalyzer 2100 to asses sample integrity. The RNA-seq library preparation from 1 μg of total RNA was performed with TruSeq RNA Sample Prep Kit (Illumina) according to the manufacturer's instructions, and 2x150 paired-end sequencing performed on Illumina Hiseq 2500. For RNA-Seq data Analysis, reads from RNA-seq experiments were aligned to zebrafish genome build danRer10 using STAR [52]. Duplicate reads were removed with SAM tools [53]. The mRNA expression of a gene was quantified by FPKM (Fragments Per Kilobase of transcript per Million mapped reads) based on RefSeq gene annotation using Cuffdiff [54]. Differential expressed genes were identified by a Fold change of ≥ 1.5 and FPKM of ≥ 3 for at least one of the two conditions based on R package edgeR. The GO enrichment and KEGG pathway analysis were performed using the R package cluster Profiler. The GO terms with FDR < 0.05 were consider as significant. RT-PCR assays were performed as described [18].

Chromatin immunoprecipitation (ChIP)-sequencing and data analysis
Chromatin was prepared from fresh heart tissues following the manual instruction of ChIP-IT High Sensitive Kit (Cat.No. 53040, Active Motif. Tokyo, Japan). ChIP assays were carried out using previously described protocols [55]. For each sample, fixed tissues were lysed to prepare nuclear extracts. After chromatin shearing by sonication, lysates were incubated overnight at 4°C with protein A Dynabeads (Invitrogen) coupled with 3-5 μg of antibody specific for H3K9Ac or H3K9Me3. Recovered DNA was analyzed by quantitative PCR (ChIP-qPCR) and sequencing (ChIP-Seq). ChIP-seq libraries were generated by ThruPLEX-FD Prep Kit (Product No. R40012, Rubicon, Ann Arbor, USA), and sequenced by Illumina Hiseq 2500 using the paired-end module and with 150 bp reads on each end. Sequence reads of 150bp were obtained, mapped to the zebrafish genome (danRer10), and processed as described previously.
All ChIP-Seq data sets were trimmed and aligned to zebrafish genome build danRer10 genome using Bowtie2 [56]. Duplicate reads were removed with SAM tools [53]. The Model-based Analysis for ChIP Sequencing v2.0 (MACS2) [57] was used to identify regions of ChIP-Seq enrichment over background in an unbiased manner. For histone modification H3K9Ac and H3K9Me3 ChIP-Seq, we modified the parameters to facilitate accurate detection of broad peaks (--broad --broad-cutoff 1E-3 -p 1E-5). The top 10000 peaks were used for known motif enrichment analysis, which were performed using HOMER (v4.9) [25] with default parameters. Other analysis was performed using deep tools [58]. To identify regions that were differentially enriched, the HOMER software was used [25]. Super enhancers (SEs) were identified using ROSE (https://bitbucket.org/young_computation/rose). Closely spaced peaks (except those within 2 kb of TSS) within a range of 12.5 kb were merged together, followed by the measurement of H3K9Ac signals. These merged peaks were ranked by H3K9Ac signal and then classified into SEs or (typical enhancers) TEs. Both SEs and TEs were assigned to the nearest refSeq genes.

Partition of subcellular fraction of blastema tissues
The ventricular apex of 60 adult fish were amputated and reared for six days. The newly regenerated heart tissues were removed into 200ul of cold PBS, and minced. The homogenate was serially centrifuged to collect large cell components at 2,000g for 15 minutes (Fraction II). The supernatant was further centrifuged at 16,000g for one hour to collect extrudants and EV pellets. For cell culture, the pellets were redissolved in 200ul PBS and homogenized. For biomolecule identification, the collected supernatant (fraction II) or pellets (fraction III) were separately dissolved in Trizol for RNA extraction and protein isolation according to the manufacturer's manual. The yield RNAs and proteins were used for transcriptome analyses and protein mass-spectrometer. The other 60 unampumtated fish were used as control.

Cell culture
Zebrafish embryonic fibroblast cellline (PAC2) were cultured in L-15 medium supplemented with 12% FBS at 32 C. To test the effects of the tissue blastema extract on PAC2 transcription, the collected subcellular fractions were separately added into the PAC2 cultures. After 36 hours, the treated cells were washed, and suspended in Trizol for RNA extraction and RNA sequencing as described above. To test Krt5-BMP signaling networks, immunohistochemistry (IHC) were conducted after Krt5-6His, recombinant human Noggin (C-Fc), and Liothyronine were added into PAC2 cell cultures at concentrations of 300ng/ml, 200ng/ml, 200 μM, respectively. To identify the interaction proteins, the Krt5-His and/or Noggin-Fc peptides-transfected cells were lysed and then precipitated by Ni-NTA Sefinose ™ resin (Sangon Biotech) and Protein A/G respectively. The precipitated proteins were identified by Sangon Biotech Co. (Shanghai) through Mass spectrometry analyses (TripleTOF 5600+) and ProteinPilot (V4.5) data analyses.

ATAC-seq library construction and data analysis
The ATAC-seq libraries were prepared as previously described [59,60] with minor modifications. Briefly, fresh cells were lysed in lysis buffer (10 mM Tris-HCl at pH 7.4, 10 mM NaCl, 3 mM MgCl2 and 0.1% IGEPAL CA-630 [Sigma-Aldrich]) for 10 min on ice to prepare the nuclei. Immediately after lysis, nuclei were spun at 500 x g for 5 min to remove the supernatant. Nuclei were then incubated with the Tn5 transposome (Vazyme Biotech, China). Following the tagmentation at 37 °C for 30 min, the products were purified with QIAquick PCR Purification Kit (Qiagen, Germany), PCR was performed as described in a previous study [59]. After the amplification, libraries were purified with DNA clean beads (Vazyme Biotech, China) to remove contaminating primer dimmers and quantitated using Qubit 2.0 Fluorimeter (Invitrogen, USA). All libraries were sequenced on the Illumina HiSeq X Ten with 150 bp paired-end reads (Novogene Biotech, China).
All sequencing data were mapped onto the danRer10 zebrafish genome assembly using Bowtie2 [56] with command line parameter -X 2000. For all data files, duplicate reads were removed with SAMtools [53], only non-duplicate reads kept in the BAM format were used for the subsequent analysis. MACS2 [57] was used for peak scanning with the parameters --shift -100 --extsize 200. The BAM files were converted to the BigWig format using the bamCoverage script in deepTools [58] with RPKM normalization. For purposes of assessing the alteration of chromatin accessibility on a genome-wide scale, deep-Tools was used to compute the average scores (RPKM normalized values) for transposase hypersensitive sites (THSs) in each sample. It was defined as hyper-accessible regions associated with increased ATAC-seq signal if the regions showed fold change ≥ 1.5 in the experimental group compared to the control group. On the contrary, the hyper-accessible regions associated with decreased ATAC-seq signal was defined if the regions showed fold change ≥ 1.5 in the control group compared to the experimental group. HOMER [25] was used for peak annotation and motif searching around the peak summits.

Image quantification and statistical analyses
IHC-stained tissue sections were imaged on a confocal microscope (Leica TCS SP8). Signal intensity was quantified using ImageJ software. All experiments were carried out in triplicate and repeated at least twice. Statistical analyses were performed using SPSS 17.0 software. Significant differences for all quantitative data were considered when * P ≤ 0.05, ** P ≤ 0.01, *** P ≤0.001, and **** P ≤ 0.0001. The antibodies and primers used are listed in Supplemental table 9.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Morphological analysis of heart regenerative process and blood clot replacement by Hematoxylin and eosin (H&E) staining. Figure S2: Immunohistochemical detection of eight common epigenetic marks during zebrafish heart regeneration. Figure S3: Subcellular localizations of H3K9Ac and H3K9Me3 within the transforming cells at 6dpa and 9dpa. Inlet showed magnification. Figure S4: Antagonist pretreatments altered H3K9Me3 and H3K9Ac depositions and Flk1-GFP distribution. Figure S5: ChIP-seq identification of H3K9Ac-and H3K9Me3-specific poised enhancers. Figure S6: Enrichments of H3K9Ac and H3K9Me3 marks at the nearest genes. Figure S7: ChIP-qPCR reevaluation of the enrichment of the two histone modification marks (H3K9Me3 and H3K9Ac) at the promoters of the target genes. Figure S8: Top 10 enriched TF motifs for H3K9Ac and H3K9Me3 ChIPseq. Figure S9: Localization of the Krt5 peptide in the Krt5 and/or noggin transfected PAC2 cells. Figure S10: Bafilomycin A1 (BA) pretreatment reduced tissue degradation and enhanced Flk-1 indicated angiogenesis at 6dpa. Table S1: Functional comparison of identified genes between blastema tissues and subcellular fractions-induced PAC2 cells. Table S2: H3K9Ac and H3K9Me3 reciprocal depositions at chromatin. Table S3: H3K9Ac and H3K9Me3 enrichments and the nearest genes' transcription. Table S4: Corresponding motifs at the predicted H3K9Ac-/H3K9Me3-specific enhancers. Table S5: Prediction of TF motifs associated with H3K9Ac-/H3K9Me3-specific enhancer at three core blastema genes. Table S6: Key blastema effectors identified by functional annotation clustering of the subcellular fraction proteins. According to human intermediate filament database (http://www.interfil.org/proteins.php), other intermediate filaments (IFs), such as, vimentin (vim), desmin (desma), and synemin (synm) were identified in the partitioned cell components. Table S7: ATAC-seq analyses of chromatin accessibility and known motif enrichment in PAC2 cells. Table S8: Mass spectrometry identifications of affinity purified 6-His-Krt5 and/or Noggin-Fc protein complex.