Gene and metabolite expression dependence on body mass index in human myocardium

We hypothesized that body mass index (BMI) dependent changes in myocardial gene expression and energy-related metabolites underlie the biphasic association between BMI and mortality (the obesity paradox) in cardiac surgery. We performed transcriptome profiling and measured a panel of 144 metabolites in 53 and 55, respectively, myocardial biopsies from a cohort of sixty-six adult patients undergoing coronary artery bypass grafting (registration: NCT02908009). The initial analysis identified 239 transcripts with biphasic BMI dependence. 120 displayed u-shape and 119 n-shape expression patterns. The identified local minima or maxima peaked at BMI 28–29. Based on these results and to best fit the WHO classification, we grouped the patients into three groups: BMI < 25, 25 ≤ BMI ≤ 32, and BMI > 32. The analysis indicated that protein translation-related pathways were downregulated in 25 ≤ BMI ≤ 32 compared with BMI < 25 patients. Muscle contraction transcripts were upregulated in 25 ≤ BMI ≤ 32 patients, and cholesterol synthesis and innate immunity transcripts were upregulated in the BMI > 32 group. Transcripts involved in translation, muscle contraction and lipid metabolism also formed distinct correlation networks with biphasic dependence on BMI. Metabolite analysis identified acylcarnitines and ribose-5-phosphate increasing in the BMI > 32 group and α-ketoglutarate increasing in the BMI < 25 group. Molecular differences in the myocardium mirror the biphasic relationship between BMI and mortality.

Elevated Body Mass Index (BMI) is an important risk factor for heart failure and cardiovascular death 1 . However, recent studies have reported a biphasic u-shaped relationship between increasing BMI and mortality in clinical settings characterized by acute metabolic stress such as cardiac surgery 2,3 , acute coronary syndromes 4 , heart failure 5 , and in patients requiring dialysis 6 . Here people with BMI between 25 and 35 have paradoxically better survival than those with low or normal BMI, or very high BMI (> 35) 2 These observations may be attributable to reverse epidemiology where people who are underweight or who have severe obesity have worse outcomes attributable to frailty or sarcopenia, or to unmeasured confounding such as fitness, or the presence or absence of metabolic syndrome. Both severe obesity and frailty lead to mitochondrial dysfunction 7 and dysregulated bioenergetics 8 . Imbalance of NADH/NAD+ ratio in sarcopenic underweight or normal patients can lead to mitochondrial dysfunction and consequently to higher levels of ROS production and low levels of chronic inflammation that affect protein and mitochondria turnover 9 . In obese patients, nutrient overload can overwhelm the TCA cycle, also leading to ROS production and inflammation 10 . We hypothesized that changes in expression of mitochondrial genes and metabolites may contribute to the biphasic association between BMI and adverse events following cardiac surgery. We tested this in human myocardial biopsies subjected to untargeted next generation sequencing and targeted metabolomics acquired as part of an ongoing observational study. The aim was to identify genes and energy metabolites whose expression in the biopsies show a biphasic response to increasing BMI and to evaluate how these differences could translate into differences in clinical outcomes.

Results
Study cohort and group definition. The ObCard study (NCT02908009) screened 558 participants between August 2017 and May 2019, out of which 221 were approached. Sixty-eight consecutive patients were recruited to the study, with one protocol violation recorded and one participant withdrew from the study (Fig. 1A). The distribution of BMI in the sample is shown in Fig. 1B.
Out of 66 collected samples, we analyzed the transcriptome in 53 and metabolites in 55 myocardial biopsies. (Fig. 1A and Table S1). Metabolomics data for one sample was not included in the analysis because we were not able to detect majority of the metabolites. After data processing, 14,967 transcripts were identified as proteincoding, 2102 were long non-coding RNAs (lncRNA), 662 were pseudogenes, and the remaining 239 transcripts included unknown transcripts, micro-RNAs, rRNA and other non-coding RNAs ( Fig. 2A). All samples were very similar in transcript biotype composition (Fig. 2B).
To identify transcripts and metabolites that display biphasic u-or n-shape BMI relationships, we used the Two Lines method developed by Simonsohn 11 . The algorithm does not impose an assumption of form on the data, which overcomes a limitation of quadratic regression known to lead to false positives. The Two Lines method tests two regression lines and detects a breakpoint between them. The analysis identified 239 transcripts with a biphasic relationship, out of which 107 had their local minimum or maximum at BMI 28 or 29 ( Fig. 3A and Table S2). None of the analyzed metabolites displayed a biphasic BMI dependence.
Given the limitations of BMI as a measure of obesity and the ad hoc nature of the WHO BMI obesity classifications we adopted a data-driven approach to group participants. Based on the analysis of transcripts, we defined the group that included local minima or maxima of gene expression as BMI 28/29 ± 3, range 25-32 (Fig. 3B). The other two groups were BMI < 25 and BMI > 32. The three patient groups were well matched with respect to baseline demographics, clinical characteristics, medications, and operative procedures (Table 1) apart from post-surgery worst multiple organ dysfunction score (MODS, higher in 25 ≤ BMI ≤ 32, and BMI > 32) and  www.nature.com/scientificreports/ non-red blood cells transfusion (lowest in BMI > 32) that were different between groups. Levels of missing data were low (Table 1). The analysis cohort included 53 participants (80.3%) with complete transcriptome data and 55 participants (83.3%) with complete metabolomics data (Table S1). Baseline characteristics in the complete case cohort were comparable to the analyzed cohort (Table S3).
Transcriptomics and metabolomics data analysis. Hierarchical clustering and principal component analyses were performed on log-transformed and quantile normalized transcriptomics data with batch effect removed. Neither analysis revealed any clear differences between groups when all transcriptomics data were included (Fig. 3C, D). However, only 15.3% of sample variability could be explained indicating that the majority of the heterogeneity is not captured by the first two components (Fig. 3D). Clustering of the metabolomics data also did not show any clear differences between groups (Fig. 3E, F). Instead, it identified samples 34 and 58 as clear outliers. These samples were removed from further analysis. Differential expression (DE) analysis identified only one significantly different transcript (AC006059.2, logFC = − 3.91, adjusted p value = 0.008) between BMI > 32 and 25 ≤ BMI ≤ 32 groups and also across all groups (adjusted p value = 0.018). AC006059.2 encodes an uncharacterized protein within the hypoxia-induced gene 1 domain.
To determine differences between groups of transcripts, we performed a differential pathway analysis using competitive gene set testing as described by Wu and Smyth 12 to identify groups of transcripts annotated to a specific pathway (Reactome and Gene Ontology) that behave differently between the BMI groups. Out of 63 identified pathways (adjusted p value < 0.05), 29 were related to translation and mRNA processing (over 300 transcripts, Fig. 4A and Table S4). These pathways were downregulated in both the 25 ≤ BMI ≤ 32 and BMI > 32 groups as compared with the BMI < 25 group. Transcripts involved in striated muscle contraction (29 transcripts; Fig. 4A and Table S4) behaved in an opposite manner. They were expressed at higher levels in both 25 ≤ BMI ≤ 32 and BMI > 32 groups as compared with the BMI < 25 group. In addition, cholesterol biosynthesis (24 transcripts), and innate immunity pathways like FCGR activation (58 transcripts) or creation of C2 and C4 activators (56 transcripts) were upregulated in the BMI > 32 group as compared with the 25 ≤ BMI ≤ 32 ( Fig. 4A and Table S4). Ribosomal proteins L3, L4 and L6, signal peptidase complex catalytic subunit (SEC11C), SMG5 and proline-rich nuclear receptor coactivator 2 (PNRC2) were the most variable between the groups (p value < 0.05, Table S5) and among translation-related transcripts. These transcripts also displayed a biphasic u-shape (or n-shape, SMG5) BMI relationship (Fig. 4B). Other transcripts that had membership in significant pathways  www.nature.com/scientificreports/ and displaying biphasic BMI dependence are plotted in Figure S1. These include CARS1 (translation), UPF3A, PPP2R1A (nonsense-mediated decay), NHP2, THUMPD1 (rRNA processing), IGKV2-28, FCN3 (creation of C4 and C2 activators), OAS3 (interferon signaling), KYAT3, RSAD2 (metabolism of amino acids) and TNNT2 (striated muscle contraction).

Left ventricular ejection fraction
Left main stem disease-n (%) 2 (13%) 9 (24%) 2 (17%) 0.761 0  (Table 2). α-ketoglutarate also has membership in one of the significant pathways, i.e., selenoamino acid metabolism (Fig. 4A). Two Lines analysis of α-ketoglutarate expression did not identify significant u or n-shape BMI dependence (Fig. 4B). However, metabolite-transcripts association analysis found its expression patterns similar to some of the transcripts with annotations in translation (RPL3, RPL4, RPL14, RSL24D1) and RNA processing (PNRC1 and RPPH1, Fig. 4C). Long and middle chain acylcarnitines expressed at highest levels in the BMI > 32 group as compared with both the 25 ≤ BMI ≤ 32 and BMI < 25 groups. There was no significant difference between the BMI < 25 and 25 ≤ BMI ≤ 32 groups. Two line analysis indicated that acylcarnitines significantly increased in patients with BMI greater than 25 ( Figure S2 and Table S2). Four out of eight acylcarnitines (C16-OH, C18-OH, C18:1-OH, C18:2-OH) were hydroxylated. Their unhydroxylated counterparts (apart from carnitine C18) did not show any significant BMI dependence (Table S2). Ribose-5-phosphate was also significantly upregulated in the BMI > 32 group as compared with the BMI < 25 group and its expression increased significantly with BMI ( Figure S2 and Table S2).

Extent of coronary disease
Next, we examined whether groups of genes without any form of pathway information were differential, and whether such groups show biphasic BMI relationship. We explored correlations between expression patterns of 429 transcripts that were most variable between the three BMI groups (p value < 0.05, Table S5) using weighted gene correlation networks analysis. The analysis grouped the transcripts into twelve networks. Eigengene values of five showed either u or n shape BMI relationship (Turquoise, Purple, Grey, Green and Yellow in Fig. 5A) and eigengene values of one network significantly increased with BMI across the whole range (Black in Fig. 5A). The networks' eigengene values for each of the six networks were different between the three BMI groups as indicated by one-way ANOVA (Fig. 5B).
To identify the biological function of the six significant networks, we submitted the specific transcripts to Reactome pathway enrichment tool 13 . Risosomal protein transcripts from the Turquoise network significantly (adjusted p value < 0.05) enriched translation and RNA processing pathways; transcripts from the Purple network (Hydroxymethylglutaryl-CoA synthase, Insulin-induced gene 1, Glutathione hydrolase 1, Tumor necrosis factor ligand superfamily member 10 and Low-density lipoprotein receptor) enriched lipids and steroid metabolism, biological oxidations, and apoptotic pathways. Transcripts from the Yellow network enriched muscle contraction pathways with tropomyosin β-chain, desmin, myosin regulatory light chain (MYL12A), myosin-binding protein C and Ca 2+ /calmodulin-dependent protein kinase type II subunit gamma (CAMK2G). Interleukin-2 receptor alpha chain from the Green network enriched 'RUNX1 and FOXP3 control the development of T reg ' . Cytochrome P450 (CYP4FA) present in the black network enriched 'Cytochrome P450 substrates' (Fig. 5C). The Grey network included ten transcripts that were not clustered in any other network and, as expected, did not significantly enrich any pathway.

Discussion
Main findings. The obesity paradox refers to a u-shaped biphasic relationship between BMI and mortality in people exposed to acute metabolic stress including cardiac surgery 2 . In the current analysis we show that groups of myocardial transcripts involved in translation and muscle contraction show biphasic BMI dependence; specifically the expression of genes involved in muscle contraction were greatest, and the expression of genes associated with translation were lowest in the 25 ≤ BMI ≤ 32 group. In addition, BMI > 32 was associated with upregulated innate immunity and cholesterol synthesis. These results were consistent across both competitive gene set analyses, and weighted gene correlation networks analysis. Contrary to our original hypothesis,  www.nature.com/scientificreports/ metabolites involved in energy production did not demonstrate biphasic BMI dependence. Instead, acylcarnitines and ribose-5-phosphate increased together with BMI. Our findings are summarized in Fig. 6.
Clinical Importance. These findings provide new evidence that BMI is associated with molecular differences in human myocardium that mirror the clinical biphasic relationship between BMI and mortality in the obesity paradox. These findings do not prove causality between obesity and improved outcomes, neither do they disprove unmeasured confounding or reverse epidemiology as explanations of the obesity paradox. Rather they support the hypothesis that the obesity paradox is attributable to an underlying molecular mechanisms, and suggest that further explorations of the therapeutic potential of these findings are warranted.
A novel finding of this study is that the expression of translation-related genes showed biphasic BMI dependence. However, it is not clear why this would provide an advantage to the 25 ≤ BMI ≤ 32 participants. Sarcopenia and frailty commonly observed in people who are under-and normal weight, or who have severe obesity, are more often associated with reductions in protein synthesis. Importantly, levels of protein synthesis were not measured in this analysis, so the results do not allow us to speculate on whether this was influenced by the observed global differences in translation. Changes in protein translation may have multiple other effects on cell metabolism. For example changes in ribosomal proteins influence processes like regulation of gene expression, cell cycle control, apoptosis, cell migration in angiogenesis, and intimal thickening (reviewed in 14 ). Dysregulation of ribosomal proteins expression is also a characteristic of cancer progression 15 . In the current study, changes in the expression of ribosomal proteins in 25 ≤ BMI ≤ 32 patients appear to regulate the expression of specific transcripts. Out of 104 transcripts in the Turquoise network (Fig. 5), at least thirteen have a clear function in translation-related processes. These results warrant further investigation of translation-related and associated transcripts and proteins to determine whether they could explain the apparent survival benefits of overweight patients in cardiac surgery.
The expression of genes involved in muscle contraction were increased 1.05-1.6-fold in the 25 ≤ BMI ≤ 32 group compared with BMI < 25 and BMI > 32 groups. These included tropomyosin β-chain, desmin, myosin regulatory light chain (MYL12A), myosin-binding protein C and Ca 2+ /calmodulin-dependent protein kinase type II subunit gamma (CAMK2G). Evidence indicates that changes in the expression of proteins responsible for muscle contraction lead to cardiomyopathies. Mutations in myosin regulatory light chains associate with disruptions of contractility and hypertrophy. Some of the mutations disrupt phosphorylation sites of the myosin regulatory light chains that is proposed to be mediated by Ca 2+ /calmodulin-dependent protein kinase and essential for muscle contraction 16 . Mutations in myosin-binding protein C also leads to hypertrophy 17 and decreased desmin expression levels associated with the progression of heart failure 18 . Imaging studies demonstrate reduced contractility in the atria of people with severe obesity, with remodeling and improved contractility following weight loss interventions 19 . These observations support the further evaluation of pre-surgery weight loss as an organ protection strategy in people with severe obesity.
The muscle contraction transcripts were linked with 26 other transcripts that included tRNA methyltransferase O (TRMO) and SMG5 nonsense mediated mRNA decay factor that are involved in RNA processing. These transcripts can potentially form a regulatory network controlling expression of muscle contraction. A question appears whether there is a functional relationship between expressions of ribosomal and muscle contraction genes.  www.nature.com/scientificreports/ The metabolite analysis identified several long and middle-length acylcarnitines that were significantly increased in the BMI > 32 group. The longest-chain acylcarnitines were hydroxylated (C16-OH, C18-OH, C18:1-OH and C18:2-OH). This type of modification of acylcarnitines was recently identified as a biomarker of mitochondrial myopathy 20 . Their unhydroxylated counterparts did not change significantly between the groups, however, the significance of this finding can be biased by the fact that increasing BMI positively correlates with circulating lipid levels. In this study, we did not measure plasma metabolites to adjust for that. Also ratios of hydroxylated/unhydroxylated acylcarnitines were not statistically different. Apart from acylcarnitines, we did not detect any significant difference in expression of glycolytic and oxidative phosphorylation metabolites. The reason can be the fact that there was no significant difference between the groups in the incidence of diabetes and heart failure between the groups, which have clear associations with BMI and where the role of mitochondria is clearly established 21 . Alternatively, it is possible that the differences become apparent after surgery. Overall, the results do not allow us to speculate on the role of mitochondria in the obesity paradox. Nevertheless, the levels of α-ketoglutarate, involved in branched amino acids catabolism in mitochondria, were increased in the BMI < 25 group and its expression patterns were similar to ribosomal proteins. The significance of this finding needs further investigation.
The study has several limitations. First, the sample size was small. As a result we could not detect significant differences in single gene expression. Second, the samples collected during surgery also potentially differed in the percentage of cardiomyocytes and could include fragments of blood vessels or fat tissue. Due to small biopsy sizes (30-100 mg), it was impossible to assess the cell-type heterogeneity of the samples. Third, we recruited a very low risk cohort to the analysis; rates of acute kidney injury, lung injury, and low cardiac output were very low. This, along with the small sample size reduced our ability to demonstrate differences in clinical outcomes between the three groups. Finally we did not collect information on glucose levels, past diseases or exercise, which would further add to the interpretation of the omics data.
In summary, our results indicate that the expression of genes involved in translation and related processes are downregulated in the myocardium of 25B ≤ MI ≤ 32 patients. In contrast genes involved in muscle contraction are overexpressed. To establish whether the two groups of genes are functionally related requires further research. Future research should also assess whether interventions targeting these processes may have translational potential for clinical myocardial protection. www.nature.com/scientificreports/ Disease: Improving Global Outcomes (KDIGO) Stage 3 AKI 22 or requiring inotropes, ventilation or intra-aortic balloon pump) were excluded. Emergency or salvage procedures were also excluded.

Methods
Sampling. Atrial biopsies (30-100 mg) were collected prior to cardiopulmonary bypass from the right atrium auricle from patients that fasted at least eight hours before the surgery. Samples were immediately snapfrozen in liquid nitrogen, split for RNA isolation and metabolomics analysis and stored at − 80 °C.

Outcomes. Levels of metabolites and transcripts in atrial biopsies.
RNA isolation and sequencing. RNA

Metabolomics. A panel of 144 metabolites involved in mitochondrial function and energy metabolism
were analyzed using a targeted assay on a Thermo Quantiva interfaced with a Vanquish Liquid Chromatography System as previously described in 23,24 . In brief, tissue was extracted using a modified Folch extraction into chloroform/methanol (2:1 600 μl per 50 mg of tissue, followed by 200 ul of water, 200 ul of chloroform, repeated once). For nucleotides and acyl-CoA derivatives one half of the aqueous extract was dissolved in 150 µl of 70:30 acetonitrile:water containing 20 µM deoxy-glucose 6 phosphate and 20 µM [U-13C, 15 N] glutamate. The resulting solution was vortexed, sonicated and centrifuged. Chromatography consisted of a strong mobile phase (A) was 100 mM ammonium acetate, and weak mobile phase was acetonitrile (B) and the LC column used was the ZIC-HILIC column from SeQuant (100 mm × 2.1 mm, 5 µm).
For amino acids and TCA cycle intermediates aqueous extracts were reconstituted in 50 μl of 10 mmol/l ammonium acetate in water before TCA cycle intermediates were separated using reversed-phase liquid chromatography on a C18-PFP column (150 mm × 2.1 mm, 2.0 μm; ACE). For chromatography on the UHPLC system, mobile phase A was 0.1% formic acid in water, and mobile phase B was 0.1% formic acid in acetonitrile. Data processing and statistical analysis. Transcriptomics. Sequencing data were quality-checked with Fastqc v0.11.5 25 , quantified with Salmon v1.21 26 after indexing and annotating with Gencode34 (Ensembl v100) reference genome and transcriptome files. Transcript quantities were normalized to length-scaled transcripts per million and filtered to retain only high quantities (filterByExpr function in edgeR) before downstream analysis using limma-voom model 27 with empirical Bayes moderation 28 with batch accounted for. The false discovery rate was set at 5%. Pathway enrichment in the dataset was tested with camera function (limma) 12 using Reactome 13 (protein-coding transcripts) and Gene Ontology 29,30 (lncRNA) annotations.
The biphasic BMI relationship of gene and metabolite expression was tested with the Two Line algorithm 11 Weighted gene correlation networks analysis (WGCNA) was performed with the WGCNA R package 31 . These analyses were performed on log-transformed and quantile normalized expression values with batch effect removed as they are not multivariate models.
Metabolomics. Peak area ratios of metabolites were obtained by integration within vendor software (Xcalibur QuanBrowser, Thermo Scientific, Hemel Hempstead, UK) and compared with isotopically labelled standards for quantification. Data for metabolites where no-expression values were greater than 20% were removed from the analysis. Further pre-processing included filtering features of constant value and extreme relative standard deviation or coefficient of variation, as well as missing value replacement using MetaboAnalystRv2 32 . The processed data was log-normalized and Pareto-scaled. Pairwise comparisons of sample groups were carried out using cross-validated PLSDA and t-test. Metabolites were considered differential where variable importance in projection scores (VIP) > 1 and t-test p value < 0.05.
Sample separation was visualized using principal component analysis plot of normalized transcriptome and metabolite data with R base and ggplot2 33 .
Multiomics analyses. RNA and metabolite were combined using block sparse PLSDA models using mixOmics v0.6 34 . Canonical correlation patterns and association networks derived from the model components were then used to infer relationship among genes and metabolites. Network visualization was carried out with Cytoscape 35 .