Spatiotemporal single-cell analysis reveals critical roles of mechano-sensing genes at the border zone in remodeling after myocardial infarction

The left ventricular remodeling after myocardial infarction (MI) causes heart failure, but its underlying mechanisms remain largely unknown. Here, we performed an integrative analysis of spatial transcriptomics and single cell RNA-seq in murine MI model and found that mechanical stress-response genes are expressed at the border zone and play a critical role in left ventricular remodeling after MI. Single-cardiomyocyte RNA-seq of the heart after MI divided cardiomyocytes into 8 clusters. We identied the cell cluster which was specically derived from the ischemic area in an early stage after MI and expressed mechano-sensing genes such as Csrp3, Ankrd1 and Ankrd2. Spatial transcriptomic analysis demonstrated that cardiomyocytes with this gene expression prole were specically localized at the border zone in post-MI day 1. AAV9-mediated gene silencing of Csrp3 exacerbated cardiac dilatation and dysfunction after MI. Collectively, our datasets and ndings not only provide an insight into spatiotemporal molecular pathogenesis of myocardial infarction but also highlight mechano-sensing genes at the border zone as an adaptive regulator of left ventricular remodeling.


Introduction
Despite decades of intensive research and therapeutic development, ischemic heart diseases are still a leading cause of death worldwide. Myocardial infarction (MI) causes various complications such as arrhythmia, valvular disease and heart failure 1 . In particular, the number of patients with heart failure has recently been increasing and becoming a more serious problem in many countries 2,3 . MI induces changes of left ventricular (LV) size, shape and function which were known as LV remodeling, leading to the development of heart failure [4][5][6] . Therefore, it is an urgent need to clarify molecular mechanisms of LV remodeling after MI to inhibit the development of heart failure and improve the prognosis of patients with MI 7-10 . However, cellular and molecular behaviors change spatially and temporally after MI, it is di cult to precisely understand the molecular mechanisms of LV remodeling.
The infarcted heart is usually roughly divided into three zones such as an infarct zone (IZ), an ischemic border zone (BZ; de ned as the sectors immediately adjacent to IZ) and a remote zone (RZ) 11 . Although previous studies showed transcriptional differences between myocardium proximal and distal to the IZ [11][12][13] , the precise characteristics of cardiomyocytes in each region have not been well explored. Furthermore, previous studies examined gene expressions using mRNA isolated from bulk samples and did not separate cardiomyocytes and non-cardiomyocytes. Here, we identi ed and characterized the distinctive transcriptional properties in each region of infarcted hearts by integrative analysis of singlecardiomyocyte RNA-seq and spatial transcriptomics. We unveiled a novel molecular adaptation of cardiomyocytes in the BZ characterized by transcriptional activation of mechano-sensing genes such as Csrp3 (also known as MLP, muscle LIM protein), Ankrd1 and Ankrd2.

Results
Single-cardiomyocyte RNA-seq identi es spatiotemporally distinct cell clusters after MI We made a murine MI model by ligating the left anterior descending artery. We separately isolated cardiomyocytes from ischemic region (IR, area including infarct zone and 2 mm of its lateral margin) and remote region (RR, area other than IR) and conducted single-cell RNA sequencing (scRNA-seq) and spatial transcriptomic analysis by Visium (10X Genomics) (Fig. 1a). We sequenced 775 cardiomyocytes from 7 murine hearts and analyzed 761 cardiomyocytes (sham, 88 cardiomyocytes; IR in post MI (pMI) day 1, 169 cardiomyocytes; RR in pMI day 1, 237 cardiomyocytes; IR in pMI day 7, 151 cardiomyocytes; RR in pMI day 7, 130 cardiomyocytes) in which >2500 genes were detected. Using clustering analysis for scRNA-seq data, cardiomyocytes derived from ve different conditions (sham, pMI day 1 IR, pMI day 1 RR, pMI day 7 IR and pMI day 7 RR) were classi ed into eight cell clusters (Fig. 1b). We statistically analyzed the proportion of each cluster and found that there were some clusters only appeared in the speci c conditions (Fig. 1c). For example, cluster 6 accounts for the majority of sham cardiomyocytes and highly expressed genes associated with mitochondrial functions and cardiomyocyte metabolism ( Fig. 1c-e). By contrast, cluster 0 and cluster 1, mainly derived from cardiomyocytes of pMI day 1 and day 7, expressed genes associated with collagen production and in ammatory responses (in ammatory cytokine production), respectively ( Fig. 1c-e). Although the proportions of many clusters were very similar between IR and RR, only cluster 7 was speci cally derived from the cardiomyocytes of IR at pMI day 1 ( Fig. 1c). Differential gene expression analysis showed that one of the key features in the cluster 7 was high expression of genes encoding for mechano-sensing proteins such as Csrp3, Ankrd1 and Ankrd2 (Fig.  1d, e). The cluster 7 cardiomyocytes occupied about 30% of all cardiomyocytes in IR of pMI day1, but they were scarcely observed in other four conditions (Fig. 1c).
Integrative analysis of spatial transcriptomics and scRNA-seq identi es spatiotemporally-regulated cell clusters after MI To know the speci c positional information of these cardiomyocyte clusters seen in our scRNA-seq analysis, we performed spatial transcriptomic analysis using Visium (10X Genomics). A mean of 1,974.4 spots (1,153-3,154 spots) from seven murine hearts (sham, n=1; day1, n=2; day7, n=2; day14, n=2) were sequenced and a median of 3,627 genes per spot (interquartile range 2,073-5,425) were detected. Fig. 2a showed both hematoxylin and eosin (H&E) staining of the heart at each time point after MI surgery and spot region clusters revealed by clustering analysis of Visium data (Fig. 2b). Although the size of infarct regions was slightly different among mice, the characteristic of each cluster's distribution was very similar and reproducible (Extended Data Fig. 1a). We also performed differential gene expression and gene ontology enrichment analysis for each cluster (Fig. 2c, d). Spot regions of cluster Green, characteristic for genes involved in mitochondrial functional pathway, occupied the majority of sham heart and RZ of infarcted hearts. Cluster Brown, characterized by expression of collagens and extracellular matrix genes, was speci cally localized in IZ (Fig. 2a, c, d). Spot regions of cluster Red were mainly seen in BZ of the infarcted hearts (Fig. 2a). Almost all marker genes of cluster Red were shared with those of cluster 7 in single-cardiomyocyte RNA-seq analysis (Fig. 1d, e, 2c, d). Genes encoding mechano-sensing proteins such as Csrp3 and Ankrd1 were highly expressed in cluster Red, which were uniquely localized in BZ of the infarcted hearts in acute phase of MI (Fig. 2e). Spatially-regulated expression of these genes in cardiomyocytes was validated by RNA in situ hybridization (Extended Data Fig. 1b-e). Subsequently, we determined cell types composing each spot region by integrative analysis of scRNA-seq data from both cardiomyocytes and non-cardiomyocytes (Extended Data Fig. 2) and Visium data. Cardiomyocytes occupied the majority of the spot regions of cluster Green and Red, whereas cardiac broblasts accounted a large part of cluster Brown (Fig. 2f).
Weighted gene co-expression network analysis of spatial transcriptomics reveals mechano-sensing genes as a distinct key factor expressed in BZ after MI To clarify the transcriptional network of BZ, we performed weighted gene co-expression network analysis (WGCNA) on the Visium data. WGCNA identi ed nine gene modules (Fig. 3a, Extended Data Fig. 3) and each module contained 40-954 genes. Among nine modules, expressions of Module 1(M1) genes encoding collagens and extracellular matrix proteins were up-regulated selectively in IZ especially in pMI day 7 (Extended Data Fig. 3). The M2 genes were widely seen in the sham hearts and RZ of the infarcted hearts (Extended Data Fig. 3). By contrast, the M3 genes were speci cally expressed in IZ of pMI day 1 hearts (Extended Data Fig. 3). Expression patterns of M8 genes were unique and up-regulated only in BZ of the pMI day1 hearts (Fig. 3b). WGCNA of M8 identi ed mechano-sensing genes such as Csrp3 and Ankrd1 as its central genes (Fig. 3c, d). M8 genes were highly overlapped with the marker genes of cluster Red (Extended Data Fig. 4). Furthermore, single-cardiomyocyte RNA-seq analysis revealed that the cardiomyocytes with lower Csrp3 expression had high expression levels of p53 signaling associated genes such as Cdkn1a, Bax and Fas (Fig. 3e). p53 signaling in cardiomyocytes was shown to induce heart failure both after MI and after pressure overload [14][15][16] , suggesting that Csrp3 in cardiomyocytes adaptively functions against the induction of failing cardiomyocytes.
Csrp3 expressed in cardiomyocytes of BZ adaptively prevents adverse LV remodeling after MI To clarify the roles of the mechano-sensing genes speci cally upregulated in BZ of the pMI day 1 heart, we generated knockdown AAV9 vectors to deliver shRNA against Csrp3, which was identi ed as one of the key genes in M8 (Fig. 3c). We injected the knock-down vectors (Extended Data Fig. 5a) into mice through their tail veins two weeks before MI surgery and performed transthoracic echocardiography to analyse their cardiac function (Fig. 4a). The delivered shRNA vectors signi cantly reduced expressions of Csrp3 at mRNA and protein levels in the heart (Fig. 4b, c). Compared with controls, the AAV9-shCsrp3 injected mice showed more enlarged left ventricles and reduced LV contraction after MI ( Fig. 4d−f). Although there was no statistical signi cance, the AAV9-shCsrp3 injected mice also tended to die after MI as compared with control mice (Extended Data Fig. 5b).

Discussion
In the present study, we comprehensively investigated the spatiotemporal gene expression of the heart after MI and identi ed several gene programs which were precisely regulated in a spatiotemporal manner. It has been well known that LV remodeling after MI including expansion of the BZ causes LV dysfunction and dilatation, resulting in heart failure 17,18 . However, the cellular and molecular mechanisms underlying remodeling processes in the BZ remain unclear. Previous studies examined expressions of RNA or proteins using myocardial tissues around the BZ dissected from infarcted hearts and reported that expressions of genes associated with angiogenesis, in ammation, apoptotic response and B-type natriuretic peptide (BNP) were highly up-regulated in the BZ 11,[19][20][21] . However, it was di cult to precisely de ne spatiotemporal molecular behaviors after MI at the single-cell level due to the ambiguous spatial localization of the BZ and the analysis using bulk tissues. In this study, we combined spatial transcriptomics and single-cell RNA-seq to obtain spatiotemporal expression pro les after MI at the single-cell level, and identi ed the mechano-sensing gene Csrp3, expressed in the BZ at an early phase after MI, as an adaptive regulator of LV remodeling. Our integrative datasets provide a publicly available resource to further deepen the understanding of pathogenesis of myocardial infarction.
Our spatial transcriptomic analysis clearly revealed that expressions of mechano-sensing genes were upregulated in an acute phase after MI in the BZ. In a recent study using echocardiography for a porcine MI model, the end-diastolic mean wall stress and myocardial stiffness of the BZ lineally increased throughout 28 days after MI 22 , suggesting that cardiomyocytes activate gene expression in response to mechanical stresses during the very acute phase of MI. Csrp3 (MLP) has been reported to be a mechanostress sensor localized at the Z disc 23,24 and its expression in cardiomyocytes is upregulated by mechanical stretch 25 as well as pressure overload on the heart 16,26 . It has been reported that Csrp3 induces expression of cardiac hypertrophic genes such as Nppb, Acta1 and Myh7 in response to pressure overload and phenylephrine stimulation through the activation of calcineurin/NFAT pathway 23,27,28 .
Rcan1, whose expression is correlated with the activity of calcineurin/NFAT pathway 29 , was also upregulated in the cardiomyocytes with high expression levels of Csrp3 in our study (Fig. 3e). By contrast, DNA damage and p53 signaling pathway, an inducer of heart failure [14][15][16] , was activated in cardiomyocytes with low expression levels of Csrp3 (Fig. 3e), suggesting that activation of Csrp3 in the BZ after MI inhibits DNA damage and p53 signaling activation to prevent the induction of failing cardiomyocytes. In skeletal muscle, Csrp3 has also been reported to regulate autophagy and protect against muscular dystrophy 30,31 . Since autophagy is one of the important cardioprotective mechanisms during post-infarction cardiac remodeling [32][33][34][35] , Csrp3 may also protect against adverse cardiac remodeling through the regulation of autophagy. In this study, mice with reduced levels of Csrp3 showed more LV dilatation and dysfunction, indicating exacerbation of LV remodeling after MI. These results suggest that the upregulation of Csrp3 as an acute response to mechanical stress in the BZ soon after MI is important to inhibit LV remodeling and may be a novel target of heart failure.

Animal models
All the animal experiments were approved by the University of Tokyo Ethics Committee for Animal Experiments and strictly adhered to the guidelines for animal experiments of the University of Tokyo (Approved Number P14-106). All wild type C57BL/6 male mice were purchased from CLEA JAPAN.
All mice were housed in separate cages at a maximum of 6 mice per cage in a speci c-pathogen-free, temperature-controlled vivarium under a 12-h light/dark cycle with ad libitum access to food and water.

Operation of myocardial infarction model and echocardiography
Myocardial infarction (MI) was induced as previously described 36 . Age of mice operated and used for experiments were 9 week old male mice. Mice were anaesthetized by inhalation of 2% iso urane. MI was induced by ligation of the left anterior descending artery. Mice that failed to develop MI or died within one day after the operation were excluded from the analysis. Transthoracic echocardiography was performed on conscious mice, using a Vevo 2100 imaging system (Visualsonics, Inc.). To minimize variation in the data, cardiac function was assessed only when the heart rate was 600-700 beats per minute. M-mode echocardiographic images were obtained from a longitudinal view to measure the size and function of the left ventricle.
Single-cardiomyocyte RNA-seq analysis of mouse samples Cardiomyocytes were isolated using Langendorff perfusion from the left ventricle as previously described 16  Single-cardiomyocyte cDNA libraries were generated using the Smart-seq2 protocol 37 and the efficiency of reverse transcription was assessed by examining the cycle threshold (Ct) values of control genes (Tnnt2) from quantitative real-time polymerase chain reaction (qRT-PCR) using a CFX96 Real-Time PCR Detection System (Bio-Rad), and the distribution of cDNA fragment lengths was assessed using LabChip GX (Perkin Elmer) and/or TapeStation 2200 (Agilent Technologies). The following primer sets were used for qRT-PCR: Tnnt2 mRNA forward: TCCTGGCAGA GAGGAGGAAG; Tnnt2 mRNA reverse: TGCAGGTCGA ACTTCTCAGC. A Ct value of 25 was set as the threshold. The sequencing libraries were subjected to paired-end 51-bp RNA sequencing on a HiSeq 2500 (Illumina) in rapid mode.
For the isolation and collection of non-cardiomyocytes, hearts were minced and enzymatically dissociated using 2 mg/mL type 2 collagenase (Worthington), 1 mg/mL dispase (Roche), and 20 U/mL DNase I (Roche), with 5 cycles of digestion for a total 40 min at 37 °C. After using a 40-μm cell strainer (Greiner) to remove the cardiomyocytes, 5,000 cells were prepared to a concentration of 1000 cells/µL and loaded into the Chromium Controller and single-cell cDNA libraries were generated using the Chromium 3′ v3 chemistry kit (10x Genomics) according to the manufacturer's instruction. Libraries were sequenced on a NovaSeq 6000 System (Illumina) using a NovaSeq S4 Reagent Kit (200 cycles,

20027466, Illumina).
Single-cardiomyocyte RNA-seq data processing Raw sequencing reads from single-cardiomyocyte RNA sequencing libraries were trimmed to remove adapter sequences and low-quality bases using fastp-0.20.0 38 with the parameters "--cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 30 --length_required 30". The reference transcript data and the gene annotation le for mouse was downloaded from GENCODE (release 23, https://www.gencodegenes.org/) 39 . The clean reads were aligned to the mouse genome (mm10) using STAR (version 2.7.2b) 40 . The reads aligned to exon were counted using featureCounts 41 . TPM normalization was calculated with reads mapped to the nuclear genome. We removed the low-quality cardiomyocytes with less than 2500 detected genes, which were used for the downstream analysis. Bright eld histology images of Hematoxylin and Eosin stain were taken using a 10X objective on a BZ-X700 microscope (Keyence Corporation, Itasca, Illinois). Raw images were stitched together using BZ-X analyzer software (Keyence Corporation) and exported as TIFF les.
Libraries were prepared according to the Visium Spatial Gene Expression User Guide (CG000239 Rev D, 10X Genomics) and sequenced on a NovaSeq 6000 System (Illumina) using a NovaSeq S4 Reagent Kit (200 cycles, 20027466, Illumina), at a sequencing depth of approximately 250-400M read-pairs per sample. Sequencing was performed using the following read protocol: read 1, 28 cycles; i7 index read, 10 cycles; i5 index read, 10 cycles; read 2, 91 cycles.

Spatial transcriptome data processing
Raw FASTQ les and histology images were processed by sample with the Space Ranger software (ver 1.2.1, 10X Genomics), against the Cell Ranger mm10 reference genome "refdata-gex-mm10-2020-A", available at: https://cf.10xgenomics.com/supp/spatial-exp/refdata-gex-mm10-2020-A.tar.gz. Raw counts were used as the input of Seurat R package (version 4.0.0), log normalization was implemented separately for each dataset and integrated them. The integrated data were used for dimensionality reduction and cluster detection. Differentially expressed genes were detected by using FindMarkers in Seurat (log2fc.threshold > 0.25 and p_val_adj < 0.05). Top DEGs were presented according to the log2 fold-change of the average expression (avg.log2FC). Finally, GO enrichment analysis was performed using Metascape with GO biological processes 43 .
Weighted gene co-expression analysis CPM (counts per million) counts were used as the input of WGCNA analysis (version 1.69) 44 . All genes expressed at more than 700 spots were used for analysis. The soft power threshold was analyzed with the "pickSoftThreshold" function and was applied to construct a signed network and calculate module eigengene expression using the "blockwiseModules" function. Modules with <30 genes were merged to their closest larger neighboring module. To visualize the weighted co-expression networks, Cytoscape (version 3.8.0) 45 with the "prefuse force-directed layout" was used. Signed eigengene-based connectivity of a gene in a module re ected node size.

AAV9-shCsrp3 infection
The AAV vectors were prepared by VectorBuilder Inc. (https://en.vectorbuilder.com) according to established procedures 46 . Brie y, AAV vectors of serotypes 2 and 9 were generated in HEK293T cells, using triple-plasmid co-transfection for packaging. Viral stocks were obtained by CsCl2-gradient centrifugation. Titration of AAV viral particles was performed by real-time PCR quanti cation of the number of viral genomes, measured as CMV copy number. The viral preparations had a titer between 1×10 12 and 5×10 12 genome copy (GC)/mL. Viruses were administered in 100-μL saline via tail-vein injection. 3×10 11 GC doses of AAV9-eGFP or 3×10 11 GC doses of AAV9-shCsrp3 were administered to the uninjured mice 2 week before MI surgery.

RNA in situ hybridization
For RNA uorescent in situ hybridization, the RNAscope system (Advanced Cell Diagnostics) was used with a probe against murine Csrp3 and Tnnt2 mRNA as previously described 47 . Frozen sections (10 µm) were xed in PBS containing 4% paraformaldehyde for 5 min at room temperature, dehydrated by serial immersion in 50%, 70%, and 100% ethanol for 5 min at room temperature, and treated with protease for 30 min at room temperature. The probe was then hybridized for 2 h at 40 °C, followed by RNAscope ampli cation, and co-stained with DAPI to detect nuclei. Images were obtained using each 20X or 40X objective on a BZ-X700 microscope (Keyence Corporation, Itasca, Illinois).

Tissue histology
For histological analysis, mice were anaesthetized by iso urane inhalation and sacri ced by cervical dislocation. The chest was opened, and the heart was ushed with cold PBS via cardiac apical insertion of a 25-gauge needle. The right atrium was cut to allow drainage of blood from the heart, and the mice were brie y perfused with cold xative (4% paraformaldehyde in PBS) through the apex of the heart. Tissues were excised, ushed with xative, incubated in xative for 12 h at 4 °C with gentle rotation, and nally embedded in para n. Para n-embedded heart tissues were sectioned into 4-μm slices using an SM2010 R Sliding Microtome (Leica Biosystems), and sections were stained with Picrosirius Red/Fast Green dyes. Bright eld histology images were taken using a 10X objective on a BZ-X700 microscope (Keyence Corporation, Itasca, Illinois).

qRT-PCR
Expression of Csrp3 mRNA was evaluated by qRT-PCR using a CFX96 Real-Time PCR Detection System, and the relative expression levels of the target genes were normalized to the expression of an internal control gene, using the comparative Ct method. The following primer sets were used for qRT-PCR.

Declarations Data Availability
The authors declare that all data are available in the article le and Extended Data information les or available from the authors upon reasonable request. The sequencing and alignment metrics of singlecardiomyocyte RNA-seq and Visium are provided as Extended data tables. The source data underlying Figures. 1e, 2c, d, 3c, d, 4b, c, f and Extended Data Fig. 3b are provided as data source tables. Singlecardiomyocyte, single-cell RNA sequencing and spatial transcriptomic data have been deposited in GSE176092 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE176092). The secure token is "ersbeoycldebpgd").

Figure 2
Integrative analysis of spatial transcriptomics and scRNA-seq identi es spatiotemporally-regulated cell clusters after MI. a, Hematoxylin and eosin (H&E) staining of whole heart at each time point after MI surgery (above). Clusters found in b. were overlaid on the H&E images above (below   Csrp3 expressed in cardiomyocytes in the BZ adaptively limits adverse LV remodeling after MI. a, Experimental design for testing the effect of AAV9-shCsrp3 vectors. b, mRNA expression levels of Csrp3 in Control and AAV9-shCsrp3 injected hearts were tested 2 weeks after injection using qRT-PCR (n = 3 vs 3).
Data are shown as mean and SD. ** P < 0.01; Signi cance was determined by unpaired t-test. c, Western blot analysis of Csrp3 using Control and AAV9-shCsrp3 injected hearts were tested 2 weeks after injection (n = 4 vs 4). Data are shown as mean and SD. * P < 0.05; Signi cance was determined by unpaired t-test.
d, Representative echocardiographic images of Control vs AAV9-shCsrp3 injected mice. e, Histochemical detection of collagen bres by Sirius Red/Fast Green dye staining in Control and AAV9-shCsrp3 injected mice at 1 week after MI surgery. f, Echocardiographic assessment of the heart of Control and AAV9-shCsrp3 injected mice after MI surgery. Data are shown as mean and SD. n = 12 (Control) vs 10 (shCsrp3). * P < 0.05, ** P < 0.01, *** P < 0.005, **** P < 0.001; signi cance was determined by two-way analysis of variance (ANOVA) with Bonferroni's multiple comparison test.