Dynamic transcriptional responses to injury of regenerative and non-regenerative cardiomyocytes revealed by single-nucleus RNA sequencing

Experimental Animals

All animal work described in this manuscript has been approved and conducted under the oversight of the UT Southwestern Institutional Animal Care and Use Committee. Timed-pregnant ICR/CD-1 mice (Charles River Laboratories) were used to deliver pups for surgical procedures on postnatal day 1 or 8. Timed-pregnant C57BL/6 mice (Charles River Laboratories) were used to deliver pups for systemic AAV9 injection on P4 and surgical procedures on P8. Neonatal Sprague-Dawley rats (Envigo) were used to isolate neonatal rat ventricular cardiomyocytes (NRVMs). Animals were housed in a 12 h light/dark cycle in a temperature-controlled room in the Animal Research Center of UT Southwestern, with ad libitum access to water and food. Sex was not determined for neonatal pups.

Cell Lines

All cells were cultured at 37°C with 5% CO₂. NRVMs were maintained in Dulbecco’s modified Eagle’s medium (DMEM) (Sigma, D5796)/199 medium (Gibco, 11150–059) (3:1) with 3% fetal bovine serum (FBS) (Gemini Bio Products, 100–106), and 1% penicillin-streptomycin (Sigma, P0781). Adeno-X 293 cells (Clontech, 632271) were maintained in DMEM with 10% FBS, and 1% penicillin-streptomycin.

Neonatal MI

Neonatal mice were anesthetized by hypothermia, and neonatal MI was performed as described previously (Mahmoud et al., 2014; Porrello et al., 2011; Porrello et al., 2013). Mice were euthanized at 1 d (24 h) or 3 d (72 h) post-surgery. Each heart was dissected in ice-cold PBS and a transverse cut on the ligation plane (for MI hearts) or similar level (for sham hearts) was made. Heart tissue below the ligation plane was collected for nuclei isolation.

Histology

Postnatal mouse hearts were fixed in 4% paraformaldehyde (Electron Microscopy Sciences, 15710) in PBS (Sigma, D8537), embedded in paraffin, and sectioned at 5 μm intervals. For trichrome staining, samples were sectioned at different levels below the suture (for MI hearts) or comparable levels (for sham hearts), and stained using Masson’s trichrome staining method.

Cardiomyocyte Nuclei Isolation

Cardiomyocyte nuclei isolation was performed as previously described (Quaife-Ryan et al., 2017) with modifications. All chemicals used were purchased from Sigma Aldrich, unless otherwise indicated. For each nuclear sample, 8–12 tissue samples collected from the same litter were pooled together and minced with a razor blade. Minced samples were resuspended in 5 mL Lysis Buffer (320 mM sucrose, 10 mM Tris-HCl, 5 mM CaCl₂, 5 mM Magnesium Acetate, 2 mM EDTA, 0.5 mM EGTA, 1 mM DTT, 1x complete protease inhibitor (Roche), 200U/mL RNase OUT Recombinant Ribonuclease Inhibitor (Invitrogen, 10777–019), pH=8) and homogenized with a dounce grinder for 15 strokes. The lysate was sequentially filtered through 70 μm and 40 μm cell strainers (Falcon, 352350 and 352340) and centrifuged at 1,000 × g for 5 min at 4 °C to pellet nuclei. The nuclear pellet was subsequently resuspended in 2 mL Sucrose Buffer (1 M sucrose, 10 mM Tris-HCl, 5 mM magnesium acetate, 1mM DTT, 1x complete protease inhibitor, 200 U/mL RNase OUT Recombinant Ribonuclease Inhibitor, pH=8) and the suspension was cushioned on top of 4 mL Sucrose Buffer, followed by centrifugation at 1,000 × g for 5 min at 4 °C to pellet nuclei. The n uclear pellet was subsequently washed once in 1 mL of Nuclei Storage Buffer (NSB) (440 mM sucrose, 10 mM Tris-HCl, 70 mM KCl, 10 mM MgCl₂, 1.5 mM spermine, 1x complete protease inhibitor, 200 U/mL RNase OUT Recombinant Ribonuclease Inhibitor, pH=7.2). The washed nuclei were resuspended in 1 mL NSB and stained with anti-PCM1 antibody (Sigma, HPA023374, 1:200) for 45 min at 4 °C. After primary antibody staining, nuclei were centrifuged at 1,000 × g for 5 min at 4 °C, washed twice with 1 mL NSB, and stained with secondary antibody conjugated with Alexa Fluor 647 (Invitrogen, A21244, 1:500) for 30 min at 4 °C. Subsequently, nuclei were washed twice with 1 mL NSB and resuspended in 2% BSA/PBS before fluorescence-activated cell sorting (FACS). For FACS, PCM1⁺ nuclei were sorted using a BD FACSARIA cell sorter in the UT Southwestern Flow Cytometry Facility. Nuclei without primary antibody staining were used as a negative gating control. After sorting, PCM1+ nuclei were pelleted at 1,500 × g for 15 min at 4 °C, resuspended in 2% BSA, and counted using a Countess II Automated Cell Counter (Thermo Scientific).

To isolate 2n and 4n cardiomyocyte nuclei, we followed the same cardiomyocyte nuclear isolation protocol with a modification of adding Hoechst (Invitrogen, H3570, 1:2,000) during the secondary antibody incubation to stain DNA. Nuclei without PCM1 antibody staining, or without Hoechst staining were used as negative controls for gating (Figure S4C–S4F).

Single Nucleus RNA-seq Library Preparation and Sequencing

Single nucleus RNA-seq libraries were generated using Single Cell 3’ Reagent Kits v2 (10xGenomics) according to the manufacturer’s protocol. Briefly, cells, reagents (10xGenomics, PN-120237), gel beads (10xGenomics, 220104), and partitioning oil (10xGenomics, 220088) were loaded into the Chromium Single Cell A Chip (10xGenomics, 2000019), targeting 2000∼5000 nuclei to be recovered per channel. Chip was loaded to the Chromium Controller (10xGenomics) in the Next Generation Sequencing Core of Children’s Research Institute at UT Southwestern to generate Gel Bead-In-Emulsions (GEMs). After reverse transcription using GEM-RT incubation protocol, GEMs were broken with recovery reagent (10xGenomics, 220016), cDNA was purified with DynaBeads MyOne Silane beads (Thermo Scientific, 37002D), amplified with PCR for 14 cycles, and further purified with SPRIselect reagent (Beckman Coulter, B23318). Quality control of recovered PCR product was performed on a Bioanalyzer (Agilent) with D5000 high sensitivity screentapes (Agilent, 5067–5592), and cDNA yield was measured by Qubit DNA high sensitivity assay (Invitrogen, Q32854). After fragmentation, end-repair, A-tailing and size selection purification, cDNA was ligated with adaptor and subsequently purified with SPRIselect reagents. Libraries were amplified using PCR with sample indexes from Chromium i7 Multiplex kit (10xGenomics, PN-120262), for a total of ∼14 sample index cycles, as estimated from original cDNA yield using Qubit assay. Final PCR products were subjected to double-sided size selection with SPRIselect reagents. Library quality control was performed on the Agilent Bioanalyzer with D1000 screentapes (Agilent, 5067–5582), and library yield was quantified by Qubit DNA high sensitivity assay. Sequencing was performed on an Illumina Nextseq 500 system operated by the Next Generation Sequencing Core of Children’s Research Institute at UT Southwestern using the 150bp high output sequencing kit (Illumina), with the following pair-end sequencing settings: Read1 – 26 cycles, i7 index - 8 cycles, Read2 – 124 cycles.

Pre-processing of snRNA-seq Data

The Cell Ranger Single-Cell Software Suit (https://support.10xgenomics.com/single-cell-geneexpression/software/pipelines/latest/what-is-cell-ranger) was used to perform sample demultiplexing, barcode processing and single-cell 3′ gene counting. The cDNA reads were aligned to the mm10/GRCm38 premRNA reference genome. Only confidently mapped reads with valid barcodes and unique molecular identifiers were used to generate the gene-barcode matrix. Further analyses for quality filtering was performed using the Seurat R package (Butler et al., 2018). Doublets were identified using Scrublet and removed from downstream analysis (Wolock et al., 2019). We further calculated the distribution of detected genes per cell and removed cells in the top 1% quantile or had fewer than 200 detected genes. We also removed cells with more than 25% of the transcripts coming from mitochondrial genes. The co-isolation of mitochondria during the nucleus isolation process is likely due their association with ER. This is consistent with reports from other groups where mitochondrial DNA was detected in single-nucleus RNA-seq (Hu et al., 2018; Lake et al., 2017) as well as ATAC-seq (which also profiles nuclei) (Rickner et al., 2019). In particular, cardiomyocytes are among cell types that have the most mitochondrial content with ∼30% of the cellular volume being filled with mitochondria (Piquereau et al., 2013). However, these mitochondrial transcripts were removed from our data before gene expression analyses, thus they do not affect our results. After quality filtering and removing unwanted cells from the dataset, we normalized the data by the total expression, multiplied by a scale factor of 10,000 and log-transformed the result. To account for variation in the number of genes and UMI detected in each cell, we normalized the scaled expression matrix using the ScaleData function and set the vars.to.regress to nUMI and nGenes.

In total, we analyzed 12,021 cardiomyocyte nuclei from control hearts and 9,716 cardiomyocyte nuclei from injured hearts. We sequenced an average of 48,896 reads with 75% mapping to the genome and a median of 925 genes per nucleus (Figure S1D).

Although the numbers and distribution of detected genes per nucleus were overall comparable across samples, we observed reduced UMI counts, possibly reflecting decreased nuclear mRNA content, and a slight increase of detected transcripts per nucleus in the injured samples, suggesting activated gene transcription after injury (Figures S1D and S1E).

Integrated Analysis of snATAC and snRNA datasets

Cardiac nuclei were isolated using the same procedure as described above. Single cell ATAC-seq libraries were generated using Chromium Single Cell ATAC library & Gel Bead Kit V1.0 (10xGenomics, 1000111) according to the manufacturer’s protocol. Briefly, nuclei were combined with ATAC enzyme to form transposition mix and were incubated for 60 min at 37°3. Transposed nuclei were then loaded to the Chromium Chip E (10xGenomics, 1000086) and run on the Chromium Controller (10xGenomics) for generation of Gel bead in emulsion (GEM). GEM was incubated and broken. GEM reaction mixture was cleaned up using Dynabeads MyOne Silane magnetic beads (Thermo Fisher Scientific, 37002D). DNA was amplified using PCR with specific primers from Chromium i7 Multiplex Kit N Set A (10xGenomics, 1000084), and purified with SPRIselect beads (Beckman Coulter, B23318). Final sequencing library size was determined using D1000 tape on the Agilent Bioanalyzer. The final library yield was measured using the Qubit DNA high sensitivity assay. Next generation sequencing was performed on an Illumina Nextseq 500 system operated by the Next Generation Sequencing Core of Children’s Research Institute at UT Southwestern using the 150bp high output sequencing kit (Illumina), with the following pair-end sequencing settings: Readl – 70 cycles, i7 index – 8 cycles, i5 index – 16 cycles, Read2 – 70 cycles.

The CellRanger ATAC v1.1.0 (https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/what-is-cell-ranger-atac) was used for primary analysis of the scATAC-seq data. Raw base call (BCL) files were converted to FASTQ files and filtered reads were aligned to the mouse mm10 reference genome. Transposase cut sites in each cell were quantified using barcoded UMIs and 10x cell barcode sequences. Count matrix for accessible chromatin peaks was generated for downstream analysis using the R package Signac (v0.1.6)

To classify cells in snATAC-seq dataset based on cell types identified in snRNA-seq, we followed the data transfer method developed in Seurat package. The detailed procedure can be found at (https://satiialab.orq/seurat/v3.1/atacseq_integration_vignette.html). In brief, LSI (Latent Semantic Indexing (Cusanovich et al., 2015) was used to learn the structure of ATAC-seq data, then FindTransferAnchors function was used to identify anchors between the scATAC-seq dataset and snRNA-seq dataset. The identified anchors were used to transfer the cell type labels learned from snRNA-seq data to snATAC-seq data with function TransferData. We filtered out cells with a prediction score < 0.4. At this threshold > 96% of cells were transferred with a cell type label from snRNA-seq dataset.

Motif enrichment analysis of differential accessible peaks

To identified differential accessible peaks between CM4 and CM1 cell clusters, we used the FindMarkers function within Seurat with parameters min.pct = 0.2, test.use = ‘LR’, and latent.vars = ‘peak_region_fragments’. Peaks that had p valued < 0.05 were designated as differential accessible peaks. We performed motif enrichment analysis using the findMotifsGenome.pl program in HOMER (Version 4.8) and searched for motif enrichment around the 200bp window of the peaks summit.

Integrated Analysis of snRNA Datasets

To account for potential batch effects, we used the IntegrateData function implemented in Seurat V3 which utilizes canonical correlation analysis that identifies a linear combination of features to construct a shared correlation structure and align the global transcriptome across datasets (Butler et al., 2018). Briefly, we performed standard preprocessing (log-normalization), and identified the top 2000 variable features for each individual dataset. We then identified integration anchors using the FindIntegrationAnchors function. We used default parameters and dimension 30 to find anchors. We then passed these anchors to the IntegrateData function to generate integrated Seurat object.

Visualization and Clustering

To visualize the data, we used Uniform Manifold Approximation and Projection (UMAP) to project cells in 2D space on the basis of the aligned canonical correlation analysis (Becht et al., 2018). Aligned canonical correlation vectors (1:30) were used to identify clusters using a shared nearest neighborhood modularity optimization algorithm (Butler et al., 2018). Using the graph-based clustering, we divided cells into 12 clusters using the FindCluster function in Seurat with resolution 0.6. We identified cardiomyocyteclusters based on expression of Myh6 and Tnnt2. We identified seven clusters expressing markers of non-cardiomyocyte populations, which likely represent contaminating non-cardiomyocyte nuclei after PCM1 sorting. These non-cardiomyocyte clusters were removed, and the clustering processes were repeated and resulted in the identification of 8 subpopulations. We merged three clusters that did not show clear molecular distinctions and were called as a single cluster when reducing the cluster resolution to 0.4. We also removed one small cluster that could not be merged based on the biological identity and was composed of too few cells (< 50) to represent a real subpopulation. The resulting datasets were divided into 5 cardiomyocyte subpopulations. Cluster markers were identified using Findmarker function with Wilcoxon rank sum test.

Analysis of 2n/4n Cardiomyocyte snRNA Datasets

We performed an additional analysis of 2n and 4n cardiomyocyte snRNA datasets. After quality filtering, we retained 3890 2n cardiomyocyte nuclei and 5540 4n cardiomyocyte nuclei. The average number of detected genes per cell in 2n and 4n datasets was 1001 and 685, respectively. To map cells in 2n and 4n datasets to the cardiomyocyte populations identified in the previous analysis, we integrated 2n, 4n, and P1-MI-d3 datasets for dimensionality reduction and subspace alignment using the first 15 canonical correlation vectors. Using the graph-based clustering, we initially divided cells into 7 clusters using the FindCluster function in Seurat with resolution 0.2. Two clusters were later identified as contaminating macrophages and endothelial cells based on their expression of C1qa and Fabp5, respectively. The identities of the remaining new clusters were assigned to the corresponding cardiomyocyte population based on 1) de novo cluster identification and 2) their integration of cells in the P1-MI-d3 dataset that had cell type labels transferred from the previous analysis. CM1, CM3, CM4, and CM5 were identified as separate clusters and their identities were assigned based on the transferred identities of P1-MI-d3 cells that were embedded in each cluster. Although the CM2 cluster was not initially identified as a separate cluster, it contains the majority of the CM2 cells from P1-MI-d3 and thus was manually selected as a separate cluster. We visualized the integrated data in UMAP (Seurat V3) and colored cells by their sample of origin (i.e. 2n, 4n, P1-MI-d3, Figure S4G) and cluster identity (Figure S4H and S4I). The proportion of each individual cluster in each sample was calculated and the Chi square test was used for enrichment analysis.

Analysis of Injury Response in CM4 Cells

CM4 cells from P1-sham-d1, P1-MI-d1 and P1-MI-d3 datasets were isolated for further cluster analysis. We recalculated the highly variable genes and performed a principal component analysis for dimensionality reduction. To visualize the data, we used UMAP to project cells in 2D space on the basis of the top 10 significant principal components, and labeled cells by their original sample IDs. We then performed pairwise comparison of transcriptomes between the three stages (P1-Sham-d1 vs P1-MI-d1, P1-Sham-d1 vs P1-MI-d3, and P1-MI-d1 vs P1-MI-d3) using the Wilcoxon rank sum test on the scaled data. Differentially expressed genes with p-value < 0.0001 in any of the comparisons were used to generate a consensus DEG list. We next plotted the z-scores of each gene’s expression in all three stages and performed unsupervised hierarchical clustering, which resulted in the identification of three groups of genes showing distinct expression dynamics (Figure 4B).

Analysis of Injury Response in CM5 Cells

To identify which cardiomyocyte population gives rise to CM5 after injury, we performed pseudotime analysis using Monocle 2 package, as described in the tutorials (http://cole-trapnell-lab.github.io/monocle-release/tutorials/) (Qiu et al., 2017; Trapnell et al., 2014). Up to 1000 cells from each major cardiomyocyte cluster in P8 hearts 1-day after MI dataset (P8-MI-d1) were isolated and used as the input expression matrix for trajectory analysis. We excluded CM4 cells from this analysis due to their low abundance. Differentially expressed genes were identified using the differentialGeneTest function with default settings. DDRTree method was used to reduce the dimension of the dataset. Cells were then ordered along pseudotime and colored by their original identities. Genes significantly changed in expression along the injury trajectory of CM5 were identified using the differentialGeneTest function with default settings. Differential gene analysis was performed by comparing CM5 in P8-MI-d3 to CM5 P8-MI-d3 (p-val <0.0001 and FC = 1.5) to identify genes that showed up-regulated and down-regulated expression in CM5 at 3-days post-MI compared to 1-day post-MI.

For comparison, we also performed pseudotime analysis using Palantir, as described in the tutorials (https://nbviewer.iupyter.org/qithub/dpeerlab/Palantir/blob/master/notebooks/Palantir_sample_notebook.ipynb) (Setty et al., 2019). Briefly, raw data of up to 400 cells from each cardiomyocyte cluster in the P8 1-day after MI dataset (P8-MI-d1) were isolated from Seurat object and used as the input expression matrix. The expression matrix was next normalized by total counts and log transformed. Principal component analysis was performed to determine metagenes which were then used to construct diffusion maps and low dimensional embedding. To determine a trajectory, the start cells were set to cells from the CM1 cluster. The computed trajectory was visualized in t-Distributed Stochastic Neighbor Embedding (t-SNE).

Cell-Cycle Analysis

Cell-cycle analysis for individual cardiomyocytes in each cluster was performed using a published approach ((Kowalczyk et al., 2015), https://satijalab.org/seurat/v2.4/cell_cycle_vignette.html). Briefly, a list of genes associated with either S or G2/M cell-cycle phase was used to assign cell-cycle scores, an index of cell-cycle activity, based on the scaled expression of these genes in each cell.

Imputation of snRNAseq Data and Analysis for Gene Expression Correlation

To correct dropout values in our datasets and to increase statistical power calculating gene expression correlation, we used Markov affinity-based graph imputation of cells (MAGIC) which is based on diffusion geometry to denoise cell count matrix and fill in missing values (van Dijk et al., 2018). For gene expression correlation analysis, MAGIC corrected gene expression matrix of CM4 cells was used to calculate expression correlation between each individual gene and Mki67 or Ccnb1 using the Spearman’s rank correlation coefficient test. Cutoffs using Spearman’s rank correlation coefficient > 0.85 and standard deviation of expression across cells greater than 0.1 were used to choose highly correlated genes Table S3. The highly correlated transcription factors were identified by intersecting with a transcription factor gene list (Zhang et al., 2012).

Upstream Regulator Analysis and Gene Regulatory Network Construction

The identified cell-cycle correlated genes in Table S3 were used as input for upstream regulator analysis. Upstream regulators were predicted using iRegulon (Verfaillie et al., 2014), to identify transcription factor (TF) binding motifs overrepresented within 10 kb genomic sequences centered around TSS for the given gene sets, with the following parameters: motif collection - 10K (9713 PWMs), putative regulatory region - 10kb centered around TSS, motif ranking database - 10 kb centered around TSS (7 species), enrichment score threshold - 3.0, ROC threshold for AUC calculation - 0.03, rank threshold - 5000, maximum FDR on motif similarity - 0.001. The TF-target gene interactions derived from this analysis were visualized using Cytoscape 3.5.1 (Shannon et al., 2003) to construct gene regulatory networks.

Cardiomyocyte in vitro Proliferation Assay

NRVMs were isolated from 1- or 2-day-old Sprague-Dawley rats with the Isolation System for Neonatal Rat/Mouse Cardiomyocytes (Cellutron, nc-6031) according to the manufacturer’s instructions. NRVMs were plated at a density of 1.5×10⁵ cells/well to gelatin-coated 12-well plates and were maintained in DMEM/199 medium (3:1, with 3% FBS), and penicillin-streptomycin for 48 h before adenoviral infection. Adenoviruses for SRF, NFYa, NFYb, NFIc, NFE2L1 and NFE2L2 were generated using the Adeno-X Adenoviral System 3 (Clontech, 632267). Packaging of adenovirus and NRVM in vitro proliferation assays were performed as previously described (Wang et al., 2019).

Cardiomyocyte in vitro Survival Assay

NRVMs were isolated and cultured as above described. NRVMs were first infected with 6–8 μL adenovirus in 1 mL DMEM/199 medium containing 3% FBS per well in a 12-well plate for 48 h and were then treated with 50 uM H2O2, or 100mM peroxynitrite (Cayman Chemical), or 0.1% Methyl methanesulfonate (MMS, Sigma) for 2 h, 1h and 1h, respectively. NRVMs treated with H₂O₂ were next washed with PBS and incubated with Calcerin AM dye (Thermo Fisher, L3224) for 20 min following the manufacturer’s protocol. Cell viability was then visualized on a Keyence BZ-X700 microscope.

Bulk RNA Sequencing in NRVMs and Data Analysis

NRVMs were collected at 48 h after infection with adenoviruses for GFP, NFYa, and NFE2L1. RNA was extracted using the RNeasy Mini Kit (QIAGEN, 74104) according to the manufacturer’s protocol. An Agilent 2100 TapeStation was used for RNA quality analysis. Stranded mRNA-seq libraries were generated using the KAPA mRNA HyperPrep Kit (Roche, KK8581) following manufacturer’s protocol. Sequencing was performed on an Illumina NextSeq 500 system using a 75-bp, high-output sequencing kit for single-end sequencing.

For data analysis, quality control of RNA-seq data was performed using FastQC Tool (Version 0.11.4). Sequencing reads were aligned to rat Rnor 6.0reference genome using HiSAT2 (Version 2.0.4) with default settings and --rna-strandness F (Kim et al., 2015). After removal of duplicate reads using Samtools v0.1.18 (Li et al., 2009), aligned reads were counted using featurecount (Version 1.6.0) per gene ID (Liao et al., 2014). Differential gene expression analysis was performed with the R package edgeR (Version 3.20.5) using the GLM approach (Robinson et al., 2010). For each comparison, genes with more than 1 CPM (Count Per Million) in at least three samples were considered as expressed and were used for calculating normalization factor. Cutoff values of absolute fold change greater than 1.5 and false discovery rate less than 0.01 were used to select differentially expressed genes between sample group comparisons. Normalized gene CPM values were used to calculate RPKM (Reads Per Kilobase per Million mapped reads) values, which were then used for heatmap plotting.

Immunohistochemistry and TUNEL Assay

Hearts were fixed in 4% paraformaldehyde in PBS (Sigma, D8537), cryopreserved in 10% sucrose/PBS and 18% sucrose/PBS solution sequentially, embedded in O.C.T. Compound (Fisher Healthcare, 23730571), and cryosectioned at 10 μm intervals. For immunofluorescent staining, cryosections were air-dried for 30 min at room temperature and fixed with 4% PFA for 20 min. Section slides were then washed twice with PBS and permeabilized with 0.3% Triton X-100 (Fisher Scientific, BP151–500) in PBS for 10 min. Sections were then blocked in 5% goat serum (Sigma, G9023)/3% BSA/0.025% Triton X-100/PBS blocking solution for 1 h and stained with the indicated primary antibodies prepared in blocking solution at 4 °C overnight using the following dilutions: cTnT (Invitrogen, MA5–12960, 1:200), Xirp2 (Thermo Fisher Scientific, PA5–57072, 1:50), pH3 (Cell Signaling Technology, 9701S, 1:200), TdTomoto (MyBioSource, MBS448092, 1:500). Sections were subsequently washed with 0.025% Triton X-100/PBS three times (5 min for each wash), and incubated with corresponding secondary antibodies conjugated to Alexa Fluor 488, 555 or 647 (Invitrogen) prepared in blocking solution at room temperature for 1.5 h. After secondary antibody incubation, sections were washed with PBS, incubated with Hoechst (Thermo Scientific, 62249, 1:2000) or DAPI (Invitrogen, D1306, 1:2000)/PBS at room temperature for 10 min, and washed twice with PBS before mounting. Images were obtained using a Zeiss LSM 800 confocal microscope. TUNEL assay was performed using Click-iT Plus TUNEL Assay for In Situ Apoptosis Detection kit (Thermo Fisher Scientific, C10619) following manufacturer’s protocol.

Generation and Injection of AAV9

AAV9 vectors expressing TdTomato, NFYa or NFE2L1 driven by the muscle-specific CK8 promoter were cloned by replacing the Cas9 coding sequence of the AAV9-CK8-Cas9 vector with the corresponding coding sequences of TdTomato, NFYa or NFE2L1, and were packaged by the Harvard Medical School/Boston Children’s Hospital Viral Core, as previously described (Amoasii et al., 2017; Martari et al., 2009). C57BL/6 mice were injected with AAV9 at P4 at a titer of 5E13 vg/kg by intraperitoneal injection.

Transthoracic Echocardiography

Cardiac function was evaluated by two-dimensional transthoracic echocardiography on conscious mice using a VisualSonics Vevo2100 imaging system as described previously (Wang et al., 2019). All measurements were performed by an experienced operator blinded to the study.

Quantification of Cardiomyocyte Size and TUNEL Positive Area

Qualification of cardiomyocyte size was performed on transverse sections using ImageJ software (https://imagej.nih.gov/ij/). Wheat Germ Agglutinin (WGA, Invitrogen, W32466) was co-stained with cTnT to mark cellular membranes of cardiomyocytes. The freehand selection tool was used to trace the cellular boundary of Xirp2⁺ and Xirp2⁻ cardiomyocytes, and the area was measured (Analysis -> Measure). Only the cardiomyocytes at the transverse view were scored. Animals (n=4) were sectioned and two sections from each animal were imaged at 4–6 non-overlapping view sights.

For quantification of the TUNEL signal intensity, heart cross-sectional area was first traced using the freehand selection tool and then measured. The TUNEL signal was measured as integrated density then normalized by the heart section area.

Quantitative Real Time PCR Analysis

Total RNA was extracted from NRVMs using Trizol and reverse transcribed using iScript Reverse Transcription Supermix (Bio-Rad) with random primers. The Quantitative Polymerase Chain Reactions (qPCR) were assembled using KAPA SYBR Fast qPCR Master Mix (KAPA, KK4605). Assays were performed using a 7900HT Fast Real-Time PCR machine (Applied Biosystems). Expression values were normalized to 18S mRNA and were represented as fold change. The following oligonucleotides were ordered from Integrated DNA Technologies to measure transcript abundance: 18S Forward: 5′- ACC GCA GCT AGG AAT AAT GGA - 3′; 18S Reverse: 5′- GCC TCA GTT CCG AAA ACC A - 3′; Ccna2 Forward: 5′- CCC GGA GCC AGA AAA CCA CTG GT - 3′; Ccna2 Reverse: 5′- GCT CAC AAG GAT GGC CCG CAT - 3′; Cdc2 Forward: 5′- TTT CGG CCT TGC CAG AGC GTT - 3′; Cdc2 Reverse: 5′- GTG GAG TAG CGA GCC GAG CC - 3′; Aurka Forward: 5′- GGG TGG TCT GTG CAT GCT CCG - 3′; Aurka Reverse: 5′- GCC TCA AAA GGA GGC ATC CCC ACT A - 3′; Aurkb Forward: 5′- ATG AGC AGC GGA CTG CCA CG- 3′; Aurkb Reverse: 5′- GTC CAG GGT GCC GCA CAT GG- 3′; Cdc20 Forward: 5′- GGC TGG GTT CCC CTG CAG ACA T - 3′; Cdc20 Reverse: 5′- TGG GCA AAG CCA TGG CCT GAG A - 3’.

Pathway and Gene Set Enrichment Analysis

Statistical analysis of gene sets was performed using the Metascape (Zhou et al., 2019). Gene Set Enrichment Analysis was used to determine the enrichment of genes that are more highly expressed either in P1 or P14 hearts as well as embryonic genes in cardiomyocytes at different stages (P1 versus P9) and in different clusters (CM4 versus CM1-CM3) (Subramanian et al., 2005). For generating the embryonic enriched gene set, the following dataset was downloaded from ENCODE portal (https://www.encodeproject.org): ENCSR597UZW, and ENCSR526SEX and analyzed using R package edgeR. Gene Set Enrichment Analysis was also performed to determine the enrichment of Hallmark pathway gene sets (Liberzon et al., 2015) in CM4 versus CM1-CM3.