The epicardium is a mesothelial tissue layer that envelops the heart. Cardiac injury activates dynamic gene expression programs in epicardial tissue, which in zebrafish enables subsequent regeneration through paracrine and vascularizing effects. To identify tissue regeneration enhancer elements (TREEs) that control injury-induced epicardial gene expression during heart regeneration, we profiled transcriptomes and chromatin accessibility in epicardial cells purified from regenerating zebrafish hearts. We identified hundreds of candidate TREEs, which are defined by increased chromatin accessibility of non-coding elements near genes with increased expression during regeneration. Several of these candidate TREEs were incorporated into stable transgenic lines, with five out of six elements directing injury-induced epicardial expression but not ontogenetic epicardial expression in larval hearts. Whereas two independent TREEs linked to the gene gnai3 showed similar functional features of gene regulation in transgenic lines, two independent ncam1a-linked TREEs directed distinct spatiotemporal domains of epicardial gene expression. Thus, multiple TREEs linked to a regeneration gene can possess either matching or complementary regulatory controls. Our study provides a new resource and principles for understanding the regulation of epicardial genetic programs during heart regeneration.
This article has an associated ‘The people behind the papers’ interview.
The zebrafish heart is capable of complete or near-complete regeneration after injury, based on proliferation of spared cardiomyocytes (CMs) (Poss et al., 2002). The pro-regenerative environment provided by non-muscle cells, such as the epicardium, endocardium, vasculature and immune cells, contributes to this potential (Cao and Poss, 2018; Gemberling et al., 2015; Gonzalez-Rosa et al., 2012; Hui et al., 2017; Karra et al., 2018; Kikuchi et al., 2011a; Lepilina et al., 2006; Masters and Riley, 2014; Wang et al., 2013). For example, genetic ablation of the epicardium, a thin mesothelial layer that envelops all vertebrate hearts, blocks heart muscle regeneration and coronary angiogenesis in zebrafish (Wang et al., 2015). Key developmentally potent genes, such as retinaldehyde dehydrogenase 2 (raldh2), T-box transcription factor 18 (tbx18), fibronectin 1 (fn1) and neuregulin 1 (nrg1), are induced in epicardial tissue upon cardiac injury, first organ-wide and then resolving to the site of trauma (Fig. 1A), in a phenomenon known as ‘epicardial activation’ (Gemberling et al., 2015; Lepilina et al., 2006; Wang et al., 2013). Understanding the gene expression responses to injury that define epicardial activation can illuminate defining aspects of heart regeneration (Cao and Poss, 2018).
Enhancers are a class of cis-regulatory elements that help orchestrate gene expression during animal development and in response to environmental changes (Pennacchio et al., 2013). Kang et al. first reported a short non-coding DNA sequence upstream of the gene leptin b (lepb) that can direct expression in zebrafish hearts and fins upon injury and during regeneration, referring to these context-preferential sequences as tissue regeneration enhancer elements (TREEs). TREEs can be identified by comparing profiles of chromatin structure or decorations from uninjured and regenerating tissues. From these profiles, sequences near genes that increase RNA levels during regeneration, and in which increased marks of activated enhancers are evident in the regeneration contexts, represent candidate TREEs. Generation of transgenic animals and/or mutant animals is essential to validate TREEs, which have been described and validated in many contexts, including zebrafish hearts, zebrafish and killifish fins, and Drosophila imaginal discs (Begeman et al., 2020; Goldman et al., 2017; Harris et al., 2016; Kang et al., 2016; Pfefferli and Jazwinska, 2017; Wang et al., 2020). The discovery of specific regulatory sequences that underlie regeneration programs can reveal candidate upstream and downstream factors in regeneration, while also providing tools with which to manipulate regeneration (Chen et al., 2020; Kang et al., 2016; Sugimoto et al., 2017; van Duijvenboden et al., 2019).
Initial studies have made it clear that different cell populations engage distinct compendia of TREEs during their respective regenerative responses (Goldman et al., 2017; Lee et al., 2020; Thompson et al., 2020; Vizcaya-Molina et al., 2018). Here, to elucidate candidate TREEs responsible for epicardial gene expression responses, we profiled the chromatin accessibility of epicardial cells during heart regeneration in zebrafish by ATAC-seq (Buenrostro et al., 2013), validating several candidate TREEs using stable transgenic reporter lines. Our study provides a resource of gene regulatory changes in epicardial cells during heart regeneration and reveals new concepts in TREE-based gene regulation.
ATAC-seq analysis of epicardial chromatin structure during heart regeneration
To identify candidate TREEs that direct epicardial gene expression, we profiled transcriptomes and whole-genome chromatin accessibility by bulk RNA-seq and ATAC-seq from purified epicardial cells. We postulated that these datasets would reveal areas of active regulation, including enhancer elements linked to genes involved in epicardial responses to injury and regeneration. To isolate epicardial cells, we used an EGFP reporter driven by the regulatory sequences of tcf21. Although epicardial cells are a heterogeneous population and tcf21 conceivably does not label the entire population, it is the best available pan-epicardial marker that labels both quiescent and injury-responding epicardial cells in zebrafish (Cao et al., 2016; Kikuchi et al., 2011a; Weinberger et al., 2020). To elicit a strong epicardial injury response, we performed ventricular resection injuries on tcf21:nucEGFP animals and collected ventricles at 3 and 7 days post-amputation (dpa) together with uninjured ventricles (Ctrl) and isolated EGFP+ cells at >95% purity (Fig. 1A,B). These two time-points represent a stage of organ-wide epicardial activation (3 dpa) and injury site-restricted activation (7 dpa, Fig. 1A) (Lepilina et al., 2006). From four biological replicates, we identified 315,000 open chromatin sites on average in each experimental group (Table S1). In 3 dpa samples, 16,676 sites displayed increased chromatin accessibility compared with samples from uninjured hearts, whereas 5588 had decreased accessibility [Fig. 1C and Table S2, false discovery rate (FDR) <0.05]. By 7 dpa, chromatin accessibility appeared to largely normalize: 2583 sites had significantly increased chromatin accessibility and 169 sites had decreased accessibility (Fig. 1C). Clustering of these differential sites across samples derived three clusters with distinct changes in chromatin accessibility. The first cluster displayed increased accessibility during regeneration in both 3 and 7 dpa samples (Fig. 1D-F, increased), while the second cluster showed reduced accessibility (Fig. 1D-F, decreased). The last cluster has mixed trends (increased or decreased) in 3 and 7 dpa samples, with the 7 dpa sample more similar to the control (Fig. 1D-F, mixed). We next analyzed peak distribution by genomic region. In all three experimental groups, about 25% of the peaks resided in promoters, ∼15% in introns, ∼5% in exons and ∼50% are intergenic. The rest are within the untranslated regions (UTRs) and immediate downstream regions (Fig. 1G). The increased regions at 3 dpa (versus Ctrl) mostly resided in the intergenic (∼55%) and intronic regions (∼25%) (Fig. 1H). Together, these results suggest that there is substantial chromatin remodeling in epicardial cells after injury, and this remodeling is more extreme at 3 dpa than at 7 dpa.
To more closely examine these data for active regulatory elements, we integrated our ATAC-seq dataset with the published histone H3K27Ac (H3 acetylation at lysine 27) signature captured from regenerating zebrafish ventricles (Kang et al., 2016). This dataset comprises profiles of two biological replicates of ventricles regenerating after partial genetic ablation of CMs (regenerating) and uninjured ventricles (control). We examined trends of differential ATAC peaks, finding that regions with either increased or decreased accessibility during regeneration bear H3K27Ac signatures in whole-ventricle samples (Fig. 2A). However, only those regions with increased accessibility correlate well with a signature of increased H3K27Ac marks in regenerating samples (Fig. 2B, regenerating/control ratio >1). By contrast, those regions decreasing do not show changes in the H3K27Ac signature (Fig. 2B, bottom, ratio=1). This analysis implicates regions of DNA with increased chromatin accessibility in epicardial cells during heart regeneration as candidate TREEs. As an example, Fig. 2C shows the genomic region that contains Wilms tumor 1 transcription factor b (wt1b), a key epicardial transcription factor induced by injury in zebrafish (Kikuchi et al., 2011a). Four peaks upstream (−13 kb, −19 kb, −20 kb and −24 kb) of the transcription start site (TSS) of wt1b demonstrated increased accessibility at both 3 and 7 dpa. Three of these peaks have increased H3K27Ac marks in the regenerating ventricle dataset (peaks 1-3). Similarly, ATAC-seq peaks with increased accessibility are observed in the genomic regions that contain aldehyde dehydrogenase 1 family, member A2 [aldh1a2, also known as raldh2; a retinoic acid (RA)-synthesizing enzyme expressed in the epicardium upon heart injury (Fig. 2D) (Kikuchi et al., 2011b; Lepilina et al., 2006)], fibronectin 1a [fn1a; an epicardially expressed gene that is required for zebrafish heart regeneration (Wang et al., 2013)] and other key epicardial transcription factors: tcf21 and tbx18 (Fig. S1) (Kikuchi et al., 2011a).
The epicardium serves as a cell source and signaling hub for heart regeneration (Cao and Poss, 2018). Paracrine signals from the epicardium that support CM proliferation include the growth factor neuregulin 1 (Nrg1) and follistatin-like 1 (Fstl1) (de Bakker et al., 2021; Gemberling et al., 2015; Wei et al., 2015). We analyzed the genomic regions that included nrg1, fstl1a or fstl1b to help understand their regulation. Our analysis indicates the presence of at least two distinct transcripts of nrg1 (nrg1-202 and nrg1-205). Two promoters and three putative enhancers are identified for nrg1 (Fig. S1E). Further motif analysis of these five regions indicates the presence of numerous binding sites for TF activator protein 1 (AP-1) subunits, retinoid X receptors (RXRA and RXRG) and retinoic acid receptors (RARA, RARB and RARG) (Fig. S2 and Table S3). This result suggests that the AP-1 complex and RA signaling may regulate nrg1 expression during heart regeneration. For the follistatin-like factors, we identified increased transcript levels and several ATAC-seq peaks with increased accessibility for both genes during regeneration (Fig. S3A,B). In situ hybridization indicated injury-induced expression in presumed epicardial cells for each gene (Fig. S3A,B).
Motif and pathway enrichment analysis of regions with context-specific accessibility changes
Recurrent consensus motifs that gain active enhancer marks are likely to contain binding sites for transcription factors (TFs). To identify candidate transcriptional regulators active in epicardial cells during heart regeneration, we assayed for enriched nucleotide motifs within regions with differential accessibility at 3 dpa (versus Ctrl) using HOMER (Heinz et al., 2010). Among the top hits of regions with increased accessibility are binding motifs of the AP-1 complex subunits, such as Atf3 (activating transcription factor 3), JunB and Fos (Fig. 2E), which are present in about half of the analyzed regions. These subunits are highly expressed in the epicardium (Table S4), and in situ hybridization results indicated expression of junba in presumed epicardial cells both before and after heart injury (Fig. 2F). The AP-1 complex was recently implicated in the control of CM gene expression during zebrafish heart regeneration (Beisaw et al., 2020). In agreement with our finding, a recent preprint also reported the increased presence of the AP-1 complex subunit binding motifs in ATAC-seq peaks preferentially accessible in the injured epicardium (Weinberger et al., 2021 preprint). However, AP-1 motifs are common within regulatory sequences (Umer et al., 2019 preprint), and thus a requirement for these motifs might not reflect specificity for regeneration-related gene expression. Other enriched motifs highlight pathways and TFs known to regulate epicardium development and/or regeneration are Tcf21 (Hu et al., 2020), Runx1 (runt-related transcription factor 1) (Koth et al., 2020), the Hippo/Yap pathway (TEA domain family members, TEADs) (Xiao et al., 2018), C/EBPb (CCAAT/enhancer-binding protein b) (Huang et al., 2012), the TGFβ pathway (Smad2/3/4) (Chablais and Jazwinska, 2012) and the Hedgehog pathway (GLI family zinc finger 2, Gli2) (Choi et al., 2013; Sugimoto et al., 2017; Wang et al., 2015). Additional implicated TFs that have not been connected to epicardial functions include NFATC1 (nuclear factor of activated T cells 1), Nrf2 (nuclear factor-E2-related factor 2), Stat3, FoxO1 (Forkhead box O1), FoxO3, FoxO6 and ZBTB7A (zinc finger and BTB domain containing 7A) (Fig. 2E). Interestingly, the top hits of regions with decreased accessibility are binding motifs of Tcf21 and WT1, which are signature TFs of the epicardium. This may suggest a transition in cell state, which may warrant further investigation. Other enriched top hits include motifs belonging to Gata6 (Kolander et al., 2014), ERG (ETS Transcription Factor ERG), Meis1 (Crespillo et al., 2021; Huang et al., 2012) and Foxo3 (Fig. 2E).
For an overview of biological functions of genes linked to dynamic chromatin regions at 3 dpa, we performed Gene Ontology (GO) enrichment analysis. Compared with epicardial cell samples from uninjured hearts, we found several enriched pathways, including the TGFβ signaling pathway and FoxO signaling pathway, that match the motif analysis results (Table S5). The FoxO pathway regulates many cellular physiological processes, such as apoptosis, cell cycle, metabolism and oxidative stress resistance, by acting downstream of growth factors, insulin, glucose, TGFβ and other stimulators (Lu and Huang, 2011; Nakae et al., 2001; Wang et al., 2014). Other enriched pathways include cellular senescence and adherens junctions. It was recently reported that p53 induces senescence in mouse epicardial cells upon heart injury, and induced cellular senescence promotes neonatal heart regeneration (Feng et al., 2019; Sarig et al., 2019). The abovementioned transcription factor Nrf2 and Stat3 are key regulators of cellular senescence (Yan et al., 2021). The enriched processes of blood vessel development and angiogenesis are consistent with the role of the epicardium in supporting revascularization during regeneration (Marin-Juez et al., 2019; Wang et al., 2015) (Fig. S4 and Table S5). Biological processes such as mesenchymal cell development and differentiation, and stem cell development and differentiation implicate the progenitor feature of the epicardial cells. Other enriched processes include heart morphogenesis, apoptotic signaling pathway, cell migration involved in heart development and, unexpectedly, axon development-related processes (Table S5). Enriched molecular functions further emphasize ligand-receptor activities, extracellular matrix binding and SMAD binding (Table S5).
Direct tests of candidate regulatory elements for enhancer activity during regeneration
To prioritize a list of candidates for functional validation of enhancer activity, we first examined the top 20 distal regions with the highest fold increases in accessibility during heart regeneration (Fig. 3A). These peaks were assigned to nearest genes, which included fn1a (Wang et al., 2013), neural cell adhesion molecule 1a [ncam1a; a gene involved in axon development (Siles et al., 2018); listed by the enrichment annotation results in Fig. S4 and Table S5], repulsive guidance molecule BMP co-receptor b [rgmb; encoding a TGFβ superfamily signaling component that participates in neuronal development (Liu et al., 2016; Samad et al., 2005) (Table S5)] and guanine nucleotide binding protein, alpha inhibiting activity polypeptide 3 (gnai3; encoding a G protein). We next combined the bulk RNA-seq and ATAC-seq datasets, and identified 3026 ATAC-peaks with their related 1008 genes and RNA levels, with both features increased at 3 dpa (Fig. 3B, red dots; Fig. 3C; Table S6). Eight of the top 20 distal regions with increased accessibility during regeneration are on the list. These eight regions are assigned to genes fn1a, rgmb, gnai3, Rho GTPase activating protein 4a (arhgap4a), procollagen-lysine 2-oxoglutarate 5-dioxygenase (plod2), phosphatidylinositol-specific phospholipase C X domain containing 3 (plcxd3), triple QxxK/R motif containing (triqk) or family with sequence similarity 98 member B (fam98b) (Fig. 3B, Figs S1 and S3). Last, we looked for conserved regions by comparing our dataset with the published zebrafish Conserved Non-genic Elements (zCNEs) database, which contains conserved regions from fish to human (Hiller et al., 2013). We identified 887 ATAC-seq peaks that contain conserved regions and gain accessibility at 3 dpa (Fig. 3D and Table S7), including a top differential peak assigned to ncam1a (Fig. 3F, enhancer 2 or E2). 565 decreased ATAC-seq peaks also contain conserved regions. Motif analysis of these differentially regulated regions yielded a comparable result (Fig. 3E) to that of all differential peaks (Fig. 2E), with the AP-1 motifs being the most positive regulators, and Tcf21 and WT1 motifs being the most negative regulators. This result may suggest a conserved epicardial regeneration program. As an example of conserved ATAC-seq peaks, Fig. 3F shows four emerging ATAC-seq regions residing in intronic regions of ncam1a, which increase accessibility at both 3 and 7 dpa. ncam1a-E2 (+181 kb), which is among the top 20 regions (Fig. 3A), contains conserved sequences and displays strong enrichment with histone H3K27Ac marks in samples of whole regenerating ventricles. ncam1a-E4 (+110 kb) also has a significantly enriched histone H3H27Ac signature. To select candidate enhancers for functional evaluations, we performed in situ hybridization for ncam1a, rgmb, gnai3, plod2, triqk, plcxd3 and arhgap4a, observing induced transcript expression for all genes apart from arhgap4a in presumed epicardial cells upon injury (Fig. S3 and the following text). Regions of the first three genes received the highest priority for functional tests in transgenic lines, given that they have not yet been implicated in epicardial biology.
Distinct ncam1a-linked TREEs direct different domains of epicardial gene expression during heart regeneration
The human homolog NCAM1 encodes a cell adhesion protein of the immunoglobulin superfamily that regulates cell-cell and cell-matrix interactions during development and differentiation processes (Duncan et al., 2021; Siles et al., 2018). We searched for ncam1a in the Zebrafish Regeneration Database that includes published transcriptome datasets of heart, fin or spinal cord regeneration (http://zfregeneration.org) (Nieto-Arellano and Sanchez-Iranzo, 2019), finding that ncam1a RNA levels increased during regeneration of all three tissues (Fig. S5). To test the efficacy of these candidate enhancers, we subcloned each regulatory region upstream of a c-fos minimal promoter and EGFP cassette (Fig. 4A). Without a regulatory sequence, this cassette has minimal expression in embryos and adult hearts, even after cardiac injury (Goldman et al., 2017). Multiple stable lines were established. We identified two stable lines for ncam1a-E2:EGFP. In 3 dpf larvae, line 1 directs strong EGFP expression in the eye, brain, spinal cord and notochord, resembling ncam1a expression in the embryo [Fig. 4B (adapted, with permission, from ZFIN) (Ruzicka et al., 2019; Thisse and Thisse, 2004; Thisse et al., 2008), Fig. 4C]. Line 2 has dimmer expression, but EGFP signals are still visible in these tissues (Fig. 4C). In addition, whole-mount images of dissected 6 dpf (days post-fertilization) hearts indicate EGFP expression primarily in the outflow tract and atrioventricular valves, but not in tcf21+ epicardial cells (Fig. 4D). After resection of the adult ventricle, we observed ncam1a expression in the ventricular surface, which is enriched around the injury site at 3 dpa, suggesting an epicardial expression pattern (Fig. 4E). Similarly, images of both ncam1a-E2:EGFP lines demonstrated injury-induced EGFP expression around the wound at 3 and 7 dpa in tcf21:H2A-mCherry+ cells (Fig. 4F-K). In addition, we observed mCherry+EGFP+ cells aligned in parallel and these EGFP+ cells expressed the mural cell marker pdgfrb (Ando et al., 2021) (Fig. 4G, arrowheads; Fig. 4J, arrows), suggesting ncam1a-E2 activity in tcf21+ perivascular cells. By 14 dpa, EGFP expression is much weaker but broadly distributed across the entire ventricular surface and primarily in the perivascular cells (Fig. 4F,Ge,f). EGFP was almost undetectable at 30 dpa (Fig. 4F). Although line 2 has lower embryonic expression than line 1, its adult heart expression pattern is consistent with line 1 (Fig. 4H,I,K; data not shown). To test whether ncam1a-E2 is active specifically in ncam1a-expressing cells, we performed hybridization chain reaction (HCR) staining (Choi et al., 2018) of heart sections carrying the enhancer reporters. As shown in Fig. 4L, ncam1a-E2:EGFP is expressed in ncam1a+ cells upon injury confirming its specific activity linked to the target gene. We also noticed ncam1a expression in EGFP− cells, suggesting that ncam1a-E2 only contributes partially to the injury-induced ncam1a expression in the heart. No definitive ncam1a expression was detected in the uninjured ventricles.
We next characterized ncam1a-E4. Three ncam1a-E4:EGFP lines displayed very weak whole-body larval EGFP expression without clear tissue specificity (Fig. 5A,B). With an anti-EGFP antibody staining in adult hearts, we consistently observed a small population of EGFP+; tcf21:H2A-mCherry+ cells on the ventricular surface after injury (Fig. 5C). Enhancers may reside upstream or downstream of the TSS of the regulated gene and may function together as a cluster to exert additive and synergistic actions (Choi et al., 2021; Hnisz et al., 2013). We asked whether inserting multiple enhancers in series would enhance the activity, by placing a second E4 element after the poly-A signal of the ncam1a-E4:EGFP construct (Fig. 5D). One stable line was identified (ncam1a-E4E4:EGFP) with expression in eye, tail and heart muscles but not larval tcf21+ epicardial cells (Fig. 5D,E). Prominent cardiac EGFP in the adults was observed only upon heart injury. At 3 dpa, strong EGFP expression in tcf21:H2A-mCherry+ cells was located at the injury site, whereas slightly weaker expression was also observed in the entire epicardium distal to the injury site (Fig. 5F,H), distinguishing its activity from the ncam1a-E2:EGFP lines. EGFP expression is reduced after 3 dpa, restricted to the injury site at 7 dpa, and undetectable at 30 dpa (Fig. 5F,G). Unlike ncam1a-E2, no definitive perivascular cell expression was observed for either ncam1a-E4E4:EGFP or ncam1a-E4:EGFP lines. HCR staining further confirmed ncam1a-E4E4:EGFP expression in ncam1a+ cells both in the wound region and the distal ventricular regions (Fig. 5I,J). Following the same strategy of stacking enhancers, we asked whether combining two ncam1a-E2 sequences would enhance its activity (Fig. S6A). Compared with ncam1a-E2:EGFP, ncam1a-E2E2:EGFP lines had relatively similar EGFP expression in the embryos (Fig. S6B,C). These ncam1a-E2E2:EGFP lines demonstrated apparently identical expression patterns and the same trend of expression levels as the ncam1a-E2:EGFP reporters (peaked at 3 dpa, decreased through 30 dpa; Fig. S6D-H). Although larval expression was variable between lines of each enhancer construct, adult epicardial expression was remarkably consistent across lines. These results suggest that ncam1a-E2 and ncam1a-E4 are epicardial TREEs that regulate ncam1a expression during heart regeneration. Interestingly, their distinct expression dynamics indicate that individual enhancers may combine to generate an overall expression domain. This is comparable with the overlapping but distinct activities of multiple developmental enhancers linked to one gene (Dickel et al., 2018; Dunipace et al., 2019; Osterwalder et al., 2018), a regulatory strategy that ostensibly applies to regeneration contexts.
A rgmb-linked TREE directs injury-induced epicardial gene expression
Rgmb, a TGFβ superfamily signaling component, has been reported to bind to BMP2 or BMP4 as a co-receptor, to pattern the developing nervous system or to inhibit renal cyst development (Liu et al., 2016; Samad et al., 2005). Although no function of rgmb in heart regeneration has been reported, bmp2b overexpression enhances CM proliferation after cardiac injury in zebrafish (Wu et al., 2016). rgmb RNA levels increased at 3 dpa [log2FC=1.74, Padj (adjusted P-value)=7.42E-12] and 7 dpa (log2FC=1.38, Padj=1.92E-06) in our datasets (Table S4). The Zebrafish Regeneration Database recorded increases in rgmb RNA levels in heart or spinal cord regeneration, consistent with a pro-regenerative function in different tissues (Fig. S7). Two putative enhancers (rgmb-E1, −68 kb; rgmb-E2, −71 kb) upstream of the TSS were identified from our analyses, with rgmb-E1 showing enriched H3K27Ac marks in whole-ventricle samples (Fig. 6A). Multiple stable lines were established for each candidate enhancer. We found that rgmb-E1 sequences directed strong EGFP expression in F1 embryos, including heart, eye and notochord expression (Fig. 6B), which resembled the published in situ hybridization findings (Fig. 6C, adapted, with permission, from ZFIN) (Ruzicka et al., 2019; Thisse and Thisse, 2004; Thisse et al., 2008). Cardiac expression in the 6 dpf larvae was restricted to the myocardium but not the epicardium (Fig. 6D). In adult hearts, injury-induced rgmb transcripts were detected on the ventricular surface (Fig. 6E). For the rgmb-E1 reporter lines, EGFP expression in adult epicardial cells (tcf21:H2A-mCherry+) was induced after a heart injury, both locally and distal to the injury site (Fig. 6F-I), including expression in pdgfrb+ perivascular cells (Fig. 6H), suggesting that rgmb-E1 is an epicardial TREE. We noticed patchy EGFP expression that was often adjacent to vessels (Fig. 6F-H). This may suggest pro-angiogenic or related functions of rgmb. HCR staining confirmed that rgmb-E1 activity is restricted to rgmb+ cells (Fig. 6J). However, numerous rgmb+EGFP− cells were observed in the epicardium, indicating that rgmb-E1 only contributes partially to the gene activity. By contrast to rgmb-E1, three stable lines of rgmb-E2:EGFP were identified without consistent EGFP expression in larval or adult hearts (Fig. 6K,L). These results indicate that, whereas rgmb-E1 is sufficient to direct injury-induced gene expression, rgmb-E2 on its own is inadequate.
Distinct gnai3-linked TREEs direct similar injury-induced epicardial gene expression
Gnai3 is a G protein that binds to G protein-coupled receptors (GPCRs) to regulate various transmembrane signaling pathways (Syrovatkina et al., 2016). The Zebrafish Regeneration Database reports increased RNA levels of gnai3 during heart, fin and spinal cord regeneration, suggesting functions in multiple regeneration contexts (Fig. S8) (Nieto-Arellano and Sanchez-Iranzo, 2019). In our datasets, gnai3 RNA was increased at 3 dpa (log2FC=2.08, Padj=6.43E-11) and 7 dpa (log2FC=0.85, Padj=0.023) (Table S4). We identified two putative enhancers within intron 1 (gnai3-E1, +8.8 kb) and intron 4 (gnai3-E2, +21 kb), both of which have H3K27Ac marks in the whole-ventricle profile (Fig. 7A). Three stable lines for gnai3-E1 and four lines for gnai3-E2 were established. Each gnai3-E1:EGFP line displayed whole-body larval EGFP expression without clear specificity (Fig. 7B). No EGFP expression was noticed in larval epicardium by 6 dpf (Fig. 7C). In adult hearts, injury-induced gnai3 transcripts were detected on the ventricular surface (Fig. 7D). From each of three gnai3-E1 lines, we consistently observed a small population of EGFP+; tcf21:H2A-mCherry+ cells on the ventricular surface after injury (Fig. 7E,F; 3 dpa, based on anti-EGFP antibody staining). HCR staining confirmed EGFP expression in gnai3-expressing cells (Fig. 7G).
gnai3-E2:EGFP is expressed in the notochord and regions of skin consistently in several lines, in addition to the weak whole-body expression (Fig. 8A). Similar to gnai3-E1:EGFP, no embryonic epicardial expression was seen for gnai3-E2:EGFP lines by 6 dpf (Fig. 8B). In adults, images of lines 1 and 4 of gnai3-E2:EGFP demonstrated injury-induced expression in epicardial cells at 3 dpa (Fig. 8C,D,F). These expression patterns are diminished by 7 or 14 dpa (Fig. 8C), mimicking the reduced RNA levels of gnai3 at 7 dpa (Fig. 7A). The injury-induced EGFP expression of line 3 was relatively weak but restricted to tcf21+ cells (based on anti-EGFP antibody staining, Fig. 8E). HCR staining confirmed that EGFP expression is restricted to but does not cover all gnai3-expressing cells (Fig. 8G). These results suggest that gnai3-E1 and gnai3-E2 are epicardial TREEs, and their similar expression dynamics indicate how the expression pattern of one gene during regeneration may receive similar regulatory instructions from two distinct TREEs. Although genes are commonly regulated by multiple enhancers with redundant activities in developmental contexts (reviewed by Kvon et al., 2021), our finding here extends this concept to tissue regeneration.
The profiling we present can guide identification of factors that bind and are upstream of candidate and validated regulatory sequences, as well as target genes linked to these sequences, to help elucidate the epicardial injury response during heat regeneration. Huang et al. previously defined epicardial enhancers by analyzing the evolutionarily conserved regions linked to epicardial genes, including Raldh2 and Wt1, and found enhancers that can direct expression both in developing hearts and in response to injury (Huang et al., 2012). By contrast, the epicardial TREEs we identified direct injury-induced but not developmental expression in the epicardium, suggesting they are customized to the epicardial regeneration machinery. Notably, a recent preprint used a similar approach to identify epicardial enhancers, reporting distinct regulatory programs during epicardial development and regeneration (Weinberger et al., 2021 preprint). Enhancers linked to genes loxa, ppfibp1a, col12a1a and mdka were found to be sufficient to direct gene expression in the embryonic epicardium, although whether they have enhancer activity during regeneration was not addressed.
In our study, we noticed variations between stable lines employing the same enhancer, likely explained by different genome insertion sites and copy numbers among lines. Establishing and analyzing multiple lines for each candidate enhancer is crucial to define the activity patterns with this transgenic strategy. We found that multiple TREEs linked to a regeneration-regulated gene can possess either matching (gnai3-E1 and gnai3-E2) or partially overlapping/complementary (ncam1a-E2 and ncam1a-E4) regulatory controls, suggesting additive, redundant or synergetic effects of enhancers during regeneration. In this regard, deleting one of the enhancers for each gene may not fully abolish the injury-induced gene expression, which could explain why deletions of single TREEs typically can cause minor or no effect on gene expression (Thompson et al., 2020). We expect that the epicardial TREEs validated here, as well as putative TREEs implicated by the profiles, can be employed to generate transgenic lines that enable injury-induced expression of gene cassettes, for example, to direct targeted expression of candidate pro-regenerative factors (Kang et al., 2016).
MATERIALS AND METHODS
Zebrafish maintenance and procedures
Adult zebrafish of the Ekkwill and Ekkwill/AB strains were maintained as described previously (Poss et al., 2002; Wang et al., 2011). Briefly, the water temperature was maintained at 28°C, and fish were kept on a 14/10 light/dark cycle at a density of 5-10 fish per liter. Animals of both sexes were used for adult (9-week-old to 12-month-old) experiments. Heart resection injury was carried out as described previously (Poss et al., 2002). No animal was excluded from analyses unless they died after the injury. Published transgenic lines used in this study were Tg(tcf21:H2A-mCherry)pd252 (Cao et al., 2017) and Tg(tcf21:nucEGFP)pd41 (Kikuchi et al., 2011a). All transgenic strains were analyzed as hemizygotes. Animal procedures were approved by Animal Care and Use Committees at Duke University and Weill Cornell Medical College.
Generation of transgenic reporters
Putative enhancer regions were amplified from genomic DNA using primers listed in Table S8 and inserted upstream of the 95 bp minimal mouse c-fos promoter directing EGFP (Fivaz et al., 2000). The entire enhancer-fos-EGFP-SV40 poly A cassette is flanked by two I-Sce meganuclease restriction sites that facilitate transgenesis (Babaryka et al., 2009). These constructs were injected into one-cell-stage wild-type embryos using standard meganuclease transgenesis techniques (Babaryka et al., 2009). To isolate stable lines, larvae were examined for EGFP expression or genotyped for EGFP insertions at 1-5 dpf. Twenty-one stable lines were established in this study (listed in Table S9).
RNA-sequencing and analysis
Partial ventricular resection was performed with tcf21:nucEGFP animals as described previously (Poss et al., 2002). Ventricles were collected at 3 and 7 dpa together with uninjured clutchmates. Thirty to 50 ventricles per group from fish of mixed sexes were dissociated using Liberase DH (Roche) and EGFP+ epicardial cells were isolated by FACS as described previously (Cao et al., 2016). After FACS sorting, we plated the isolated cells in dishes or on a coverslip and observed >95% purity of EGFP+ cells. We collected 110,000-170,000 cells from each group of two biological replicates. Total RNA was purified using a Qiagen RNeasy Plus Micro Kit. The library was constructed using SMARTer Ultra Low Input RNA Kit for Sequencing and sequencing was performed using an Illumina HiSeq 2000, with 82-94 million 50 bp single-end reads obtained for each library. Sequences were aligned to the zebrafish genome (danRer10) using TopHat2 (Kim et al., 2013). Transcript levels were quantified using HTseq (Anders et al., 2015). Differential expression analysis was carried out using the Bioconductor package DESeq2 (Love et al., 2014).
ATAC sequencing and analysis
Partial ventricular resection was performed with tcf21:nucEGFP animals. Ventricles were collected at 3 and 7 dpa together with uninjured clutchmates. Four biological replicates were performed for each group. Ventricle dissociation and epicardial cells isolation were done as described above. ATAC-seq libraries were made from 50,000-198,000 epicardial cells per sample as described previously (Buenrostro et al., 2013). Sequencing was performed using an Illumina HiSeq 2000 with 89-167 million 50-nt single reads obtained for each cell library.
ATAC-seq reads were trimmed for adaptors before aligning to the DanRer10 genome using Bowtie with two mismatches allowed and mapping up to four sites (Langmead et al., 2009). Reads mapped to mitochondrial DNA were removed, and alignments were filtered for PCR artifacts. Peaks were called using MACS v2.1.4 with P<0.01 and shift −37 bp, and extend 73 bp (Zhang et al., 2008). Normalized UCSC browser tracks were generated by conversion of bam format alignments to bp resolution bigWig files and scaling by the total number of mapped reads. The Bioconductor package DiffBind V3.0.15 was used for signal normalization and differential analysis (Ross-Innes et al., 2012). A fold change greater than two and P<0.05 were used to filter the significant differential peaks. Genomic distribution of peaks was analyzed using the Bioconductor package ChIPpeakAnno v3.27.2 (Zhu et al., 2010). Packages ChIPpeakAnno (Zhu et al., 2010) and complexHeatmap v2.6.2 (Gu et al., 2016) were used to annotate the peaks and generate heatmaps. The differential peaks were annotated by nearest gene start site to the center of peaks. ATAC-Seq peaks were paired to RNA-Seq differential expression data by annotated nearest gene symbols. Genome-wide motif analysis was performed using HOMER v4.9.1 (Heinz et al., 2010). The motif search of nrg1-related sequences was performed using motifmatchr (v 1.12.0) (Schep, 2021) with merged motifs (distance 10 bp) stored in package enhancerHomologSearch (v 1.0.0) (Ou et al., 2021). The dandelionPlot was created by trackViewer v 1.27.15 (Ou and Zhu, 2019). The conserved open-chromatin regions were identified by comparing our dataset with the published zebrafish conserved non-genic elements (zCNEs) database for at least 1 bp overlapped regions (Hiller et al., 2013). Ventricle Chip-seq (H3K27Ac) dataset was downloaded from Gene Expression Omnibus (GEO) under accession number GSE75894.
Histology and microscopy
Freshly collected hearts or larvae were fixed with 4% paraformaldehyde (PFA) for 1.5-2 h at room temperature or overnight at 4°C. Fixed hearts were either mounted with Fluoromount G (Southern Biotechnology, 0100-01) or PBS between two coverslips (allowing imaging of both ventricular surfaces) or applied to 10 µm cryosections. Immunostaining of whole-mounted hearts was carried out as described previously (Cao et al., 2017; Wang et al., 2015). The primary antibody used in this study was mouse anti-EGFP (ThermoFisher, A11120, 1:200). The secondary antibody used in this study was Alexa Fluor 488 goat anti-mouse (ThermoFisher, A11029, 1:200). In situ hybridization was performed on 10 µm cryosections of PFA-fixed hearts using digoxigenin-labeled cRNA probes as described previously (Poss et al., 2002). Probes were cloned from 2 dpf zebrafish cDNA using the primers listed in Table S8. Hybridization chain reaction (HCR) was performed using synthesized probes from Molecular Instruments following the published protocols (Choi et al., 2018). Images of fluorescent transgenes in live embryos were captured using a Zeiss Axiozoom V16 microscope. Whole-mount or sectioned heart tissues were imaged using a Zeiss LSM 800 confocal microscope or a Leica DMi8 compound microscope. Whole-mount in situ hybridization images of zebrafish embryos were retrieved from the Zebrafish Information Network (ZFIN) (Ruzicka et al., 2019; Thisse et al., 2008).
We thank the Zebrafish Core Facility staff at Duke University and Weill Cornell for fish care, and M. Harrison and J. Kang for comments on the manuscript.
Conceptualization: K.D.P., J.C.; Methodology: G.E.C., J.C.; Investigation: Y.C., Y.X., J.J.B., J.O., L.S., A.S., T.C.; Resources: G.E.C., K.D.P., J.C.; Writing - original draft: Y.C.; Writing - review & editing: K.D.P., J.C.; Supervision: K.D.P., J.C.; Project administration: K.D.P., J.C.; Funding acquisition: Y.X., K.D.P., J.C.
This work was supported by a Rudin Foundation Fellowship to Y.X., by National Institutes of Health (NIH) grants (R01HL131319 and R35HL150713) and an American Heart Association Merit Award to K.D.P., and by an American Heart Association Career Development Award (18CDA34110108), a Weill Cornell Start-up fund and a National Institutes of Health grant (R01HL155607) to J.C. Deposited in PMC for release after 12 months.
The RNA-seq and ATAC-seq datasets have been deposited at GEO under accession number GSE89444M.
Peer review history
The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.200133
The authors declare no competing or financial interests.