Revealing tissue-specific metabolic crosstalk after a myocardial infarction

Myocardial infarction (MI) promotes a range of systemic effects, many of which are unknown. Here, we investigated the alterations associated with MI progression in heart and other metabolically active tissues (liver, skeletal muscle, and adipose) in a mouse model of MI (induced by ligating the left ascending coronary artery) and sham-operated mice. We performed a genome-wide transcriptomic analysis on tissue samples obtained 6- and 24-hours post MI or sham operation. By generating tissue-specific biological networks, we observed: (1) dysregulation in multiple biological processes (including immune system, mitochondrial dysfunction, fatty-acid beta-oxidation, and RNA and protein processing) across multiple tissues post MI; and (2) tissue-specific dysregulation in biological processes in liver and heart post MI. Finally, we validated our findings in independent MI cohorts, using both bulk and single-cell transcriptomic data. Overall, our integrative analysis highlighted both common and specific biological responses to MI across a range of metabolically active tissues.

signaling pathway, and several metabolic pathways (glycogen, inositol phosphate and purine) 248 (Table S6). 249 250 Mitochondrial dysfunction in the heart leads to disturbance of energy (ATP) production 251 (Kiyuna et al., 2018;Palaniyandi et al., 2010) and, in the presence of oxygen, to 252 accumulation of reactive oxygen species (ROS), which can cause oxidative stress. Vdac1, a 253 key gene for regulation of mitochondria function and one of the central genes in the heart-254 specific cluster (see above), is significantly downregulated in MI (Camara et al., 2017). Vdac1 255 is located in the outer mitochondrial membrane and is involved directly in cardioprotection 256 (Schwertz et al., 2007) within the cGMP/PKG pathway ( Figure S3A). In the same pathway, 257 we also observed down-regulation of the reporter metabolite hydrogen peroxide (Table S5), a 258 ROS that is related to cardioprotection (Schwertz et al., 2007;Yada et al., 2006). We also 259 observed downregulation of Pdha1, which is known to have a substantial role in both the 260 HIF-1 signaling pathway and the pyruvate metabolism pathway that converts pyruvate to 261 acetyl-CoA in the mitochondria ( Figure S3B). Acetyl-CoA is used in the TCA cycle to 262 produce NADH and FADH2, which are both needed for ATP production and were 263 downregulated in our reporter metabolite analysis of the heart. Our findings are thus 264 consistent with dysfunctional mitochondria and ATP production in the heart in response to an 265 MI. Pdha1 has been also been linked to the heart sensitivity during to ischemic stress, where 266 its deficiency can compromise AMP-activated protein kinase activation (Sun et al., 2016). 267

268
In skeletal muscle and adipose tissue, we found that central genes in their respective tissue-269 specific clusters related to fatty acid metabolism and lipid metabolism were significantly 270 altered (Table S6, Figure 5). In liver-specific cluster, we found that their central genes were 271 related to fatty-acid beta oxidation (Cyp4a31, Cyp4a32) and glutathione metabolism (Gstm3) 272 (Table S6, Figure 5A). Alterations of fatty acid beta-oxidation and glutathione metabolism 273 have previously been reported in non-alcoholic fatty liver disease, a known risk factor of 274 CVD (Alexander et al., 2019;Mardinoglu et al., 2017a). Moreover, in liver, we also found 275 that retinol metabolism was uniquely related to genes in the liver-specific cluster, mainly 276 driven by four significantly differentially expressed central genes of the clusters, i.e. 277 Cyp26a1, Cyp4a31, Cyp4a32, and Hsd17b6 (Table S6). A previous study showed that 278 mortality from CVD in older individuals was accompanied by impaired liver ability to store 279 retinol (Lima et al., 2018). 280 281

Multi-tissue modeling reveals key metabolic pathways affected post MI 282
To investigate the metabolic responses to MI in and across tissues in the mice, we constructed 283 a multi-tissue genome-scale metabolic model. The model consisted of five tissue-specific 284 genome scale metabolic models, namely heart, liver, skeletal muscle, adipose, and small 285 intestine. The small intestine model (for which we do not have transcriptomic data) was added 286 to include ingestion and conversion of dietary nutrients into chylomicrons, which are directly 287 secreted into blood and transport lipids to other tissues . The final 288 mouse multi-tissue model included 19,859 reactions, 13,284 metabolites, 7,116 genes and 41 289 compartments. We predicted the metabolic fluxes in mice 24 h after an MI or sham operation 290 by integrating the dietary input, tissue-specific resting energy expenditure and transcriptomics 291 data. 292

293
The modeling showed that oxygen uptake, carbon dioxide production and the oxidative 294 phosphorylation pathway in heart, adipose and skeletal muscle were decreased in MI mice, in 295 agreement with the downregulation of oxidative phosphorylation we observed in these tissues 296 (Table S7). By contrast, liver showed slightly increased oxygen uptake, which might due to 297 the slightly (not statistically significant) upregulated oxidative phosphorylation (Table S7). 298 These findings indicate that the changes in oxygen and carbon dioxide fluxes and the 299 oxidative phosphorylation pathway could serve as a positive control for predicting the 300 changes due to MI in the fluxes. 301 302 Next, we investigated the tissue-specific metabolic flux changes in the same model (Table  303 S7). We found that the pentose phosphate pathway was upregulated in heart 24 hours post MI, 304 consistent with upregulated glucose metabolism after an MI. Elevated glycolysis could allow 305 the heart to rapidly generate energy under stress conditions, and the enhanced pentose 306 phosphate pathway could increase the NADPH level, which could help maintain the level of 307 reduced glutathione in heart (Tran and Wang, 2019). We also found that adipose tissue 308 secreted more ketone bodies, including acetoacetate and butyrate, into plasma; the plasma 309 level of ketone bodies has been reported as a stress marker in acute MI (Miyamoto et al., 310 1999). Notably, relatively small metabolic changes were found in liver and skeletal muscle, 311 which is probably due to the small number of transcriptomic changes in metabolic pathways 312 in these tissues. 313 314

Single-cell analysis highlights cell-type patterns among central genes in a heart-specific 315
cluster 316 We then questioned whether and how our gene set from both central and heart-specific 317 clusters identified by bulk RNA-seq data would perform at the single cell level. We examined 318 the DEGs from most central and heart-specific cluster genes in an independent single-cell 319 dataset of heart failure in human (Wang et al., 2020). We observed that the majority of the 320 heart-specific genes clustered together in both cardiomyocyte-(CM) and non-CM-enriched 321 data during heart failure, with higher expression levels in CM-enriched data ( Figure S4). On 322 the other hand, the genes from the most central cluster showed no particular clustering 323 pattern. We also analyzed the clustered tissue-specific genes, and found three distinct groups: 324 group 1, which consisted of genes related to the TCA cycle, glycolysis and HIF-1 signaling 325 pathways; group 2, which consisted of genes with higher expression in all cell types and were 326 related to heart-specific functions and propanoate and pyruvate metabolism; and group 3,

Validating our findings with publicly available datasets 347
We validated our observations in heart tissue in two independent cohorts of bulk RNA-seq 348 data from mouse heart (Table S8). We filtered both validation cohorts to get and analyzed 349 only 24 hours post-MI data. We found that there were 2169 DEGs in heart 24 h after 350 infarction from our data were validated in at least one of the independent cohort (959 in both) 351 ( Figure 6G). We also found that 109 out of the 123 most connected genes in our heart-352 specific cluster were also significantly differentially expressed in at least one of the 353 independent cohorts (81 in both). By performing functional analysis of the validation cohorts, 354 we found that ~61% of GO BP and 84% of KEGG pathways identified in our analysis of the 355 heart were also present in at least one of the validation cohorts 24 h after infarction ( Figure  356 6H-I). In both cohorts, we observed downregulation of mitochondrial functions and fatty acid 357 metabolism processes. We also observed upregulation of processes and pathways related to 358 retinol metabolism and inflammatory response in both validation cohorts. 359 360

Identification of driver genes in MI 361
We observed that Flnc, Lgals3, Prkaca and Pprc1 showed important response to MI. These 362 genes were 4 of 16 genes that were DEGs in at least three tissues and validated in both 363 validation cohorts (Table S9). Flnc, Lgals3 and Pprc1 were upregulated in heart, skeletal 364 muscle, and adipose, whereas Prkaca was downregulated in these three tissues. We further 365 retrieved their neighbors at each tissue specific CNs, showed their regulations from 366 differential expression results, and performed functional analysis in Table S9. 367 368 Flnc, which encodes filamin-C, was part of heart and skeletal muscle-specific CN cluster and 369 also part of gene group 2 at the single-cell level ( Figure S4). Its neighbor genes were found to 370 be significantly (FDR < 0.05) associated to several functions, including TCA cycle, pyruvate 371 metabolism, glycolysis pathway, and involved in mitochondrial functions. Specifically, they 372 were related to heart-specific processes in heart, VEGF signaling pathway in muscle, 373 carbohydrate metabolism in adipose, and to MAPK signaling pathway and muscle contraction 374 in heart and muscle.  CVD has a complex etiology and is responsible for a range of systemic effects, hindering our 395 understanding of its consequences on different tissues. Here, we took advantage of the 396 technological advances in high-throughput RNA-seq and applied integrative network analyses 397 to comprehensively explore the underlying biological effects of MI. Specifically, we 398 generated RNA-seq data from heart, liver, skeletal muscle and adipose tissue obtained from 399 mice 6 and 24 h after an MI or sham operation. We used transcriptomics data analyses 400 (differential expression, functional analysis, and reporter metabolites analysis) to determine 401 the systemic effects of the MI across multiple tissues. Moreover, we performed CN analyses 402 to pinpoint important key and tissue-specific clusters in each tissue, and identified the key 403 genes in each cluster. Finally, we used a whole-body modelling approach to identify the 404 crosstalk between tissues and reveal the global metabolic alterations, before finally validating 405 our findings with publicly available independent MI cohorts. 406 407 Based on our analyses, we observed downregulation of heart-specific functions and 408 upregulation of lipid metabolism and inflammatory response in heart, muscle, and adipose 409 tissue after an MI ( Figure 4B). Liver showed a distinct response with respect to the other 410 three tissues, including downregulation of inflammatory response. We observed that fatty acid 411 metabolism was downregulated in heart and adipose tissue, whereas fatty acid beta-oxidation 412 was upregulated and glutathione metabolism was downregulated in liver. We also observed 413 upregulation of oxidative stress in heart and skeletal muscle. We also observed 414 downregulation of mitochondrial functions in heart, muscle, and adipose tissue. Furthermore, 415 we found upregulation of retinol metabolism in heart and downregulation of retinol 416 metabolites in liver ( Figure 4B). 417

418
We hypothesized that downregulation of fatty acid metabolism from adipose tissue was due to 419 exchange of fatty acids with other tissues (liver and muscle) ( Figure 4B). We also observed 420 the flow of retinol from liver to heart during MI, consistent with previous reports (Palace et 421 al., 1999). These MI-associated alterations lead to dysfunctional mitochondria and decreased 422 energy production, especially in heart and skeletal muscle. In summary, we systematically unveiled the deregulation of biological processes and 445 pathways that resulted from MI in heart, liver, muscle, and adipose tissue by integrating 446 transcriptomic data and the use of biological networks. We also identified the key clusters and 447 central genes using generated tissue-specific CNs. In this study, we demonstrated a strategy to 448 utilize multi-tissue transcriptomic data to identify alteration of biological processes and 449 pathways to systemically explore the effect of a disease. 450 451

Limitation of the Study 452
We recognized several limitations to be noted on this research. First, only transcriptomic data 453 was analyzed in this research, hence the sensitivity might be limited especially for short 454 timepoint, e.g. 6 hours after MI. Second, we focused our analysis in this research only on 455 protein-coding genes. Third, to explore more about the shift in metabolism due to MI, longer 456 timepoints needs to be explored. This opens new opportunities for future research, including 457 analyzing the non-protein-coding gene signatures and longer timepoints.  infarction. The mice were then anesthetized with isoflurane, orally intubated, and connected 520 to a small-animal ventilator (SAR-830, Geneq, Montreal, Canada) distributing a mixture of 521 oxygen, air and 2-3% isoflurane. ECG electrodes were placed on the extremities, and cardiac 522 rhythm was monitored during surgery. An incision was made between the 4th and 5th ribs to

Echocardiography in mice 542
Echocardiographic examination, using VisualSonics VEVO 2100 system (VisualSonics Inc, 543 Ontario, Canada), which includes an integrated rail system for consistent positioning of the 544 ultrasound probe was performed 6 and 24 h after an MI to determine the size of the MI. We 545 calculated infarct size based on wall motion score index (WMSI) 24 h after myocardial 546 infarction by a 16-segments model on 3 short axis images, as 0 for normal, ½ for reduced wall 547 thickening and excursion in a segment and 1 for no wall thickening and excursion in a 548 segment. WMSI was calculated as the sum of scores divided by the total number of segments. 549 Hair removal gel was applied to isofluorane-anesthetized (1.2%) mice chest to minimize 550 resistance to ultrasonic beam transmission. The mice were then placed on a heating pad and 551 extremities were connected to an ECG. A 55 MHz linear transducer (MS550D) was used for 552 imaging. An optimal parasternal long axis (LAX) cine loop of >1000 frames/s was acquired 553 using the ECG-gated kilohertz visualization technique. Parasternal short axis cine-loops were 554 acquired at 1, 3, and 5 mm below the mitral annulus. Infarct size was calculated based on wall 555 motion score index 6 and 24 hours after myocardial infarction by a 16-segments model on 556 LAX and 3 short axis images view, as 0 for normal, ½ for reduced wall thickening and 557 excursion in a segment and 1 for no wall thickening and excursion in a segment. The data 558 were evaluated using VevoStrain™ software system (VisualSonics Inc, Ontario, Canada). The output from Kallisto, both estimated count and TPM (Trancript per kilobase million), 575 were subsequently mapped to gene using the mapping file retrieved from Ensembl BioMart 576 website, by filtering only protein coding genes and transcripts. Genes with mean expression 577 less than 1 TPM in each condition were filtered. For data exploration, we used PCA from 578 sklearn package (Pedregosa et al., 2011) in Python 3.7 and used TPM values as the input. 579 Subsequently, we performed differential gene expression analysis using DESeq2 package in 580 R. We utilized the capabilities from DESeq2 to normalize the rounded estimated count data 581 and to correct for confounding factors (such as time). To define a gene as differentially 582 expressed (DEGs), a gene has to fulfill a criterion of FDR < 5%. The results of differential 583 expression analysis were then used for functional analysis. 584 585 We checked the tissue specificity of the DEGs in each tissue with the data from Mouse Gene 586 Atlas (Su et al., 2004). For all the tissue-specific genes, we also checked their human-587 homolog genes in the human secretome database (Uhlén et al., 2019).

Co-expression network generation 598
We generated the co-expression network by generating gene-gene Spearman correlation ranks 599 within a tissue type, using spearmanr function from SciPy (Jones et al., 2001) in Python 3.7. 600 Using the same environment, we performed multiple hypothesis testing using Benjamini-601 Hochberg method from statsmodels (Perktold et al., 2017). Correlation data were filtered with 602 criterion of adjusted p-value < 5%. 603

604
The top 5% of filtered correlation results were then loaded into iGraph module (Csardi and 605 Nepusz, 2006)  to get the enriched GO BP and KEGG pathways. Criterion FDR < 0.05 were used to find the 609 significantly enriched terms. Clusters with less than 30 genes were discarded, to be able to get 610 significant functional analysis results. Since GO BP was relatively sparse, we used Revigo 611 (Supek et al., 2011) to summarize the GO BP into a higher level. Revigo was further 612 employed to build a GO BP network. Clustering coefficient was calculated based on the 613 average local clustering coefficient function within iGraph. 614 615

Multi-tissue metabolic modeling 616
We combined tissue-specific models (of heart, liver, muscle, adipose and small intestine) 617 constructed previously  in a multi-tissue model by adding an 618 additional compartment representing the plasma, which allows the exchange of metabolites 619 among different tissues. Blocked reactions that could not carry fluxes (and the unused 620 metabolites and genes linked to these reactions) were removed from the models. In addition, 621 the dietary input reactions and constraints were added to the small intestine model to simulate 622 the food intake (Table S7). Specifically, we assumed that the mice weighed 30 g and 623 consumed 4.5 g chow diet per day (15 g/100 g body weight) based on a previous study 624 (Kummitha et al., 2014). We also calculated the tissue-specific resting energy expenditures 625 and set them as mandatory metabolic constraints based on previous studies and resting energy 626 expenditure for other tissues was incorporated by including a mandatory glucose secretion 627 flux out from the system with the lower bound calculated based on ATP (Table S7)

Data and code availability 648
All raw RNA-sequencing data generated from this study can be accessed through accession 649 number GSE153485. Codes used during the analysis are available on 650 https://github.com/sysmedicine/ArifEtAll_2020_MultiTissueMI 651 Table S1 Differential Expression Analysis Results 653