Abstract
The exact mechanisms that underlie the pathological processes of myocardial ischemia in humans are unclear. Cardiopulmonary bypass with cardioplegic arrest allows the authors to examine the whole transcriptional profile of human left ventricular myocardium at baseline and after exposure to cold cardioplegia-induced ischemia as a human ischemia model.
The authors obtained biopsies from 45 patients undergoing aortic valve replacement surgery at baseline and after an average of 79 min of cold cardioplegic arrest. Samples were RNA sequenced and analyzed with the Partek® Genomics Suite (Partek Inc., St. Louis, MO) for differential expression. Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA) and Biobase ExPlain (Biobase GmbH, Wolfenbuettel, Germany) systems were used for functional and pathway analyses.
Of the 4,098 genes with a mean expression value greater than 5, 90% were down-regulated and 9.1% were up-regulated. Of those, 1,241 were significantly differentially expressed. Gene ontology analysis revealed significant down-regulation in immune inflammatory response and complement activation categories and highly consistent was the down-regulation of intelectin 1, proteoglycan, and secretory leukocyte peptidase inhibitor. Up-regulated genes of interest were FBJ murine osteosarcoma viral oncogene homolog and the hemoglobin genes hemoglobin α1 (HBA1) and hemoglobin β. In addition, analysis of transcription factor–binding sites revealed interesting targets in factors regulating reactive oxygen species production, apoptosis, immunity, cytokine production, and inflammatory response.
The authors have shown that the human left ventricle exhibits significant changes in gene expression in response to cold cardioplegia-induced ischemia during cardiopulmonary bypass, which provides great insight into the pathophysiology of ventricular ischemia, and thus, may help guide efforts to reduce myocardial damage during surgery.
Abstract
The authors have performed a comprehensive genomic analysis of human heart exposed to cold cardioplegia-induced myocardial ischemia providing novel information regarding gene regulation during ischemia. The results further our understanding of the pathological response of the human heart to stress and identify novel interventional targets for future therapies.
Supplemental Digital Content is available in the text.
Myocardial ischemia remains a major source of perioperative morbidity and mortality
The mechanisms involved in the pathological process of myocardial ischemia at the level of gene expression in the human heart remain elusive
The authors have performed a comprehensive genomic analysis of human heart exposed to cold cardioplegia-induced myocardial ischemia providing novel information regarding gene regulation during ischemia
The results further our understanding of the pathological response of the human heart to stress and identify novel interventional targets for future therapies
MYOCARDIAL ischemia occurs when there is a mismatch between coronary oxygen delivery and metabolic requirements of the myocardium, which may be clinically manifested during angina, coronary angioplasty, or cardiopulmonary bypass (CPB). Myocardial ischemia may lead to a spectrum of myocardial stunning, hibernating myocardium, and ultimately cell death if the ischemic insult is severe. In the human heart, irreversible damage begins after approximately 20 to 40 min of oxygen deprivation.1 Observed molecular and cellular changes of myocardial ischemia are characteristic of an inflammatory response,2 but the exact mechanisms that underlie this pathological process are unclear and may not be fully emulated by animal models of ischemia or infarction.
Thus, we felt it valuable to investigate a human ischemia model. During cardiac surgery, CPB with aortic cross-clamping (AoXC) and cardioplegic arrest is associated with excellent clinical outcomes and suitable operative conditions. However, despite the use of cardioprotective strategies, AoXC during CPB is accompanied by a variable yet obligate ischemic period lasting 1 to 3 h, resulting in hypoxia, metabolic substrate depletion, reperfusion injury, apoptosis, and necrosis.3 Cardiac-specific biomarkers of ischemia and infarction, including troponin, are increased even after routine coronary artery bypass graft surgery and correlate with the duration of ischemia from AoXC.4–6
This process of CPB provides us with the ability to examine the transcriptional profile before and after an expected, consistent, and reproducible human ischemic event although induced by cold cardioplegic arrest and not coronary occlusion. In addition, the absence of reperfusion in this time period allows us to examine the transcriptomic response to intermittent ischemia, without having to account for the perturbations of reperfusion injury. Although various animal models have been used to examine the effects of ischemia on cardiac function, no human data exist that examine the early transcriptomic response to a left ventricular (LV) ischemic insult. We therefore characterized the effect of cold cardioplegia-induced acute ischemia on the transcriptional profile of the LV by performing whole-transcriptome next-generation RNA sequencing (RNA-seq) in patients undergoing cardiac surgery by sampling human LV tissue before and after the obligate ischemia during AoXC. We hypothesized that the cold cardioplegia-induced ischemic injury will dramatically alter transcription in the human myocardium, and that we would identify genes and pathways, which will identify interventional targets for pharmacological therapy.
Materials and Methods
Patients presenting for elective aortic valve replacement surgery by a single surgeon provided written informed consent with institutional review board (Boston, Massachusetts) approval.
Human LV Samples
We obtained punch biopsies (approximately 3 to 5 µg total RNA content) from the site of a routinely placed surgical vent in the anterolateral apical LV wall of patients undergoing elective aortic valve replacement surgery with CPB. After 79 min (interquartile range, 65 to 104) of AoXC and intermittent cold blood cardioplegia for myocardial protection every 20 min,7 a second biopsy was obtained in the same manner (fig. 1).
Sample acquisition timeline. The baseline sample is acquired within minutes of the first dose of cardioplegia and application of the aortic cross-clamp. The postischemia sample is obtained shortly before the end of aortic cross-clamping. The red arrows show the timed application of cardioplegia every 20 min.
Sample acquisition timeline. The baseline sample is acquired within minutes of the first dose of cardioplegia and application of the aortic cross-clamp. The postischemia sample is obtained shortly before the end of aortic cross-clamping. The red arrows show the timed application of cardioplegia every 20 min.
RNA Sequencing
Tissue samples were immediately placed in RNAlater® (Ambion; ThermoFisher Scientific, Waltham, MA) and after 48 h at 4°C were stored at −80°C until RNA extraction. Total RNA was isolated with TRIzol (Invitrogen, Carlsbad, CA), and RNA quality was assessed using the Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA) with no samples being excluded for poor quality. Library preparation and sequencing has been described previously (see Methods, Supplemental Digital Content 1, https://links.lww.com/ALN/B130), but briefly, ribosomal RNA was removed by performing 1 to 2 washings of RNA annealed to poly-T oligos beads (Invitrogen; Life Technologies, Grand Island, NY). RNAs were reverse transcribed using random hexamers (Invitrogen). Double-stranded DNA synthesis was performed using Pol I and RNA-ase H. Short fragments were purified with QiaQuick polymerase chain reaction (PCR) extraction kit (Qiagen, Hilden, Germany) and resolved with elution buffer for end reparation and poly(A) addition and then ligated with sequencing adaptors for cluster generation and sequencing on the Illumina HiSeq 2000 (Illumina, San Diego, CA). As samples were analyzed at different times, different read lengths were used, initially using single-end reads and then transitioning to paired-end reads ranging from 36 to 100 base pairs (see table 1, Supplemental Digital Content 1, https://links.lww.com/ALN/B130, which shows the description of samples at the two sequencing centers).
Digital Droplet PCR Validation
Validation of expression levels was performed in six patients, on both the baseline and the postischemia samples, using the Bio-Rad QX-100 digital droplet PCR (ddPCR) system (Pleasanton, CA). ddPCR assays were designed to target genes (MGST1, intelectin 1 [ITLN1], FBJ murine osteosarcoma viral oncogene homolog [FOS], hemoglobin β [HBB]), and normalized to glyceraldehyde 3-phosphate dehydrogenase (Gapdh). Absolute concentration of complementary DNA was determined using standard ddPCR protocols. The ratio of target gene to Gapdh concentrations provides normalized gene expression.
Technical Replicates
We performed technical replication with two samples (baseline and postischemia) at a single sequencing center using two different library creation methods.
Analysis
Raw reads produced by the Illumina sequencer imaging files were filtered to remove reads containing adaptor sequences, containing more than 5% unknown nucleotides, or having more than 50% of reads with base quality scores less than 5.
Using Tophat v2.0.58,9 and Bowtie2,10 reads were aligned to the Homo Sapiens reference genome (UCSC hg19). For paired-end reads, mate inner distance was set at 165 base pairs and mate standard deviation (SD) at 37 base pairs. For other parameters, we used the default settings in Tophat and Bowtie2. Final read alignments having more than two mismatches or more than two gaps or more than two edits distances were discarded. The anchor length was 8 without any mismatches being tolerated. The minimum and maximum intron lengths were 70 and 500,000, respectively. The maximum insertion and deletion lengths were 3.
Analysis was performed on aligned BAM files using Partek® Genomics Suite 6.6 (Partek Inc., St. Louis, MO). Partek® Genomics Suite estimates the most likely relative expression levels of known isoforms and reads that map to the genome using an expectation/maximization algorithm, similar to that described by Xing et al.11 To identify genes with differing relative abundances of isoforms across samples, similar to Marioni et al.,12 fold change was calculated as the normalized ratio of mapped reads between baseline and postischemia.
Multivariate testing was performed with a mixed-model ANOVA with the computation of P values based on least-squares means, adjusted for other factors. To build a parsimonious ANOVA model that would be generalizable for other studies, we included only clinical and demographic variables considered to have the greatest importance based on the literature and prior published studies. These variables were age (in years), sex, presence of diabetes mellitus (defined as an HgbA1C ≥6.5 or documented in the chart on oral or parenteral treatment), presence of coronary artery disease, LV ejection fraction, time of sample (baseline vs. postischemia), center where RNA-seq was performed, and minutes of AoXC. Subsequently, these variables were then examined for sources of variation by plotting the mean signal-to-noise ratio (F ratio) of all sequenced genes. Based on their low F ratio, age, diabetes mellitus, LV ejection fraction, coronary artery disease, or sex were not included in the final model. Their exclusion also did not change the significance of the genes included in the final model. With the outcome of absolute gene expression measured in reads mapping to the genome per kilobase of transcript per million reads sequenced (RPKM), we then adjusted for time of sample, center where RNA-seq was performed, and duration of AoXC. Significance of differential expression was defined as the expression difference between frequently expressed transcripts (>5 reads per million; RPKM) between postischemia and baseline at a false discovery rate (FDR) of P value less than 0.05.13
Analysis of Transcription Factor–binding Sites
To identify which transcription factors might be important in regulating gene expression, we searched for reported and putative transcription factor–binding sites (TFBSs) in the promoter region (−500 to +100 from transcription start site) of differentially expressed genes using Biobases ExPlain MATCH program (Biobase GmbH, Wolfenbuettel, Germany).14,15 We assumed that coexpressed genes are coordinately regulated by a limited number of transcription factors. Overrepresented bindings sites, in our case the top 100 up-regulated and down-regulated genes, were compared with an unaltered background set, which in our analysis included 1,000 genes with a nonsignificant fold change of 1.0 ± 0.15. TFBSs are predicted from positional weight matrices from the TRANSFAC® (TRANScription FACtor; Biobase GmbH, Wolfenbuettel, Germany) database.16 Enriched matrices are shown as the ratio of the average number of putative binding sites per 1,000 base pairs for the query and the background sets (yes/no ratio). Values greater than 1 indicate overrepresentation of the motif in the query set. For validation, we performed all analyses three times using a background set of 500 genes picked randomly from all genes with a nonsignificant fold change of 1.0 and no overlap between sets. Overrepresentation was determined using a one-tailed Fisher exact test. We controlled for multiple testing using an FDR.
Functional Analysis of Differentially Expressed Genes
We used gene ontology (GO), which provides an extensive ontological description of biological processes, molecular functions, and cellular components, to test for enrichment of categories in sets of genes or proteins. GO ANOVA was performed using the same covariate model as for differential expression described above. Probe sets were excluded with a maximum signal less than three and an FDR P value greater than 0.05, which reflects the differential expression and disruption of GO functional categories instead of individual genes. Analysis was restricted to functional groups between 3 and 50 genes. To compare individual genes to all other genes in a functional group, we used a Z test to calculate the disruption score. A low disruption score implies a change in the pattern of gene expression within the functional group, which may or may not be accompanied by differential expression.
A GO enrichment analysis was performed using Fisher exact test restricting analyses to functional groups with more than two genes. The enrichment score was calculated using a chi-square test comparing the proportion of the gene list in a group to the proportion of background in the group, with a value of 3 or greater corresponding to significant overexpression. Furthermore, a gene set enrichment analysis, which compares the frequency of genes detected in the experiment versus random chance, was performed restricting analysis to functional groups between 15 and 500 genes and using 10,000 permutations.17
The lists of differentially expressed transcripts were uploaded to Ingenuity Pathway Analysis 9.0 (Ingenuity Systems, Redwood City, CA) and BioBase ExPlain system (Biobase GmbH) to identify regulatory networks operant in human ischemic myocardium. In both applications, the P value is calculated by determining the probability that a disease or biological function assigned to the data is random. To show the connectivity between proteins, and classify them according to the therapeutic target, biomarker, or molecular mechanism, differentially expressed genes and transcripts are categorized by Medical Subject headings (MeSHG) terms and proprietary ontology. The significance of the association is measured by the ratio of the number of proteins from the dataset that map to the category divided by the total number of proteins that map to a canonical pathway.
Results
Forty-five, predominantly male subjects with a median age of 72 yr were enrolled (table 1). The patients represent the average of our patient population, who often present with diabetes mellitus (36%), coronary artery disease (46%), and preserved LV systolic function (median ejection fraction 60%). The majority of patients presented with severe aortic stenosis (91%) with minimal aortic regurgitation (64%), mild LV septal thickening, and normal LV dimensions. AoCx time varied depending on whether patients had concomitant coronary artery bypass graft surgery (46%) in addition to their aortic valve replacement. Although the median postoperative day 1 creatine kinase MB fraction of 23.0 µg/l (interquartile range, 19.2 to 27.9 µg/l) was well below the 10 times the 99th percentile of acute myocardial infarction by conventional cardiac biomarker methods (66 µg/l in our laboratory), two patients exceeded this value.18
Gene Expression
Median RNA integrity score was 7.0 (interquartile range, 6.3 to 7.5). After quality control measures, we examined 35,470 transcripts from 21,308 genes. When we compared the postischemia sample to the baseline sample, 16,821 genes (78.9%) were down-regulated, and 4,460 genes (20.9%) were up-regulated. When we restricted to genes with a mean RPKM of greater than 5 across all samples (N = 4,098), 3,724 (90.0%) were down-regulated and 374 (9.1%) were up-regulated (fig. 2).
Volcano plot. Fold change on X-axis plotted against P value on Y-axis. The horizontal line marks the unadjusted P = 0.05 and the red vertical lines mark a fold change of ±1.5. All dots in the upper left and right quadrant are considered significant before adjusting for multiple comparisons.
Volcano plot. Fold change on X-axis plotted against P value on Y-axis. The horizontal line marks the unadjusted P = 0.05 and the red vertical lines mark a fold change of ±1.5. All dots in the upper left and right quadrant are considered significant before adjusting for multiple comparisons.
Of those genes with an RPKM greater than 5, we identified 1,241 genes that were significantly differentially expressed (FDR P < 0.05) between postischemia and baseline after multivariate modeling. The top 16 genes along with their baseline and postischemia expression changes with a mean RPKM greater than 5, fold change greater than 2, and an FDR P value less than 0.05 are shown in table 2.
Shown Are Top Up-regulated and Down-regulated Genes between Postischemia and Baseline for the Ventricle at Baseline and Postischemia and for the Atrium at Baseline

In all analyses, we applied an FDR (Benjamini & Hochberg algorithm) to present the significant P value after multiple testing. Alternatively, we could have used more stringent multiple testing procedures such as the Bonferroni correction, which can be somewhat conservative if there are a large number of tests and/or the test statistics are positively correlated, as is often the case in genetics where many genes are inherited together and have a strong relationship to each other (strong linkage disequilibrium). Had we used a strict Bonferroni correction, only one gene (hemoglobin α1 [HBA1]) in table 2 would have retained significance.
Functional Classification
GO Analysis.
To understand their collective biological function, we examined differentially expressed genes for enrichment of GO categories. The majority of gene sets with significant enrichment in a specific GO category were down-regulated, and were categorized in cellular components and biological processes, with few identifying with molecular functions (see table 2, Supplemental Digital Content 1, https://links.lww.com/ALN/B130, which shows the up-regulated and down-regulated GO annotation). The top five down-regulated GO terms with an FDR significance less than 4.7 × 10−9 were extracellular component, complement activation, protein activation cascade, immune effector process, and humoral immune response. Highly consistent across almost all down-regulated GO groups identified was the presence of complement genes (e.g., C1, C3, C6, C7).
However, significantly up-regulated gene sets were all categorized as biological processes and were associated with GO terms of regulation of catalytic activity and molecular function (see table 2, Supplemental Digital Content 1, https://links.lww.com/ALN/B130). In total, of the GO groups that met an FDR threshold of P value less than 0.01, 55 were down-regulated and only 4 were up-regulated. Furthermore, the significant P values were much lower in the down-regulated groups as well.
Upon direct comparison between the up-regulated and down-regulated gene sets, the highly significant down-regulated GO categories were associated with the immune and inflammatory response and complement activation. Conversely, the only significantly up-regulated GO category was vasculature development (see tables 3 and 4 in Supplemental Digital Content 1, https://links.lww.com/ALN/B130, which show the GO annotation of the differential expression and the gene set enrichment analysis, respectively, and see figs. 1 and 2 in Supplemental Digital Content 1, https://links.lww.com/ALN/B130, which show a graphical representation of the functional analysis of the GO terms).
GO Disruption.
The GO analysis revealed a substantial number of functional categories with low disruption scores, implying changes in the pattern of gene expression within the functional group and potential disruption of normal operation of the group. The most pronounced disruption occurred in the molecular function group encompassing the categories of antioxidant, peroxidase, and oxidoreductase activity containing the hemoglobin genes HBA1 and HBB (fig. 3). Although the expression of the majority of genes in this category is very similar between baseline and postischemia (23.3 vs. 20.1), HBB and HBA1 are significantly up-regulated compared with all other genes in this category (disruption P < 0.001).
Gene ontology group disruption (A) oxidoreductase activity and (B) oxygen transport. Shown are the two groups with the most pronounced disruption of a functional category in the gene ontology ANOVA, implying changes in the pattern of gene expression within the functional group and potential disruption of normal operation of the group. The hemoglobin genes hemoglobin α1 (HBA1) and hemoglobin β (HBB) disrupt the molecular function group oxidoreductase activity and oxygen transport. Although the expression of the majority of genes in these categories is very similar between baseline and postischemia, HBB and HBA1 are significantly up-regulated compared with all other genes in this category (disruption P < 0.001).
Gene ontology group disruption (A) oxidoreductase activity and (B) oxygen transport. Shown are the two groups with the most pronounced disruption of a functional category in the gene ontology ANOVA, implying changes in the pattern of gene expression within the functional group and potential disruption of normal operation of the group. The hemoglobin genes hemoglobin α1 (HBA1) and hemoglobin β (HBB) disrupt the molecular function group oxidoreductase activity and oxygen transport. Although the expression of the majority of genes in these categories is very similar between baseline and postischemia, HBB and HBA1 are significantly up-regulated compared with all other genes in this category (disruption P < 0.001).
Because the nature of the experiment was to examine a timed insult, a GO category of particular interest was “response to stress,” with stressors most commonly associated with oxygen, such as anoxia, hypoxia, or hyperoxia (fig. 4). Oxidative stress was also one of the most significant toxicity lists in Ingenuity Pathway Analysis (P = 4.4 × 10−5; fig. 5).
Gene ontology enrichment. Forest Plot of postischemia compared with baseline in response to stress category. Green bars show down-regulation, red bars show up-regulation with dark colors showing percent of significantly down-regulated or up-regulated genes (false discovery rate <0.05).
Gene ontology enrichment. Forest Plot of postischemia compared with baseline in response to stress category. Green bars show down-regulation, red bars show up-regulation with dark colors showing percent of significantly down-regulated or up-regulated genes (false discovery rate <0.05).
Oxidative stress toxicity list in Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA). Shown are groups with percent down-regulation (green) and up-regulation (red) matching the dataset. The numbers represent the number of genes in each category and the orange line represents the negative logarithmic P value. PPARα = peroxisome proliferator–activated receptor α.
Oxidative stress toxicity list in Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA). Shown are groups with percent down-regulation (green) and up-regulation (red) matching the dataset. The numbers represent the number of genes in each category and the orange line represents the negative logarithmic P value. PPARα = peroxisome proliferator–activated receptor α.
TFBS Analysis
To identify which transcription factors are involved in the regulation of gene expression in differentially expressed genes, we scanned the promoter region for overrepresentation of TFBSs. Similar to the GO processes, TFBSs are overrepresented in down-regulated genes; we identified 10 TFBSs that were significantly overrepresented in the promoter region of up-regulated genes and 40 TFBSs that were overrepresented in down-regulated genes (table 3).
Closer examination of the promoter region of the most significantly down-regulated genes (fig. 6) shows multiple matrices matching to similar composition of TFBS although their location in relevance to the transcription start site can vary widely between different promoters. The transcription factors, such as RORA, STAT1, STAT5a, and CEBP, act in regulation of reactive oxygen species production,19 apoptosis,20 immunity,21 and inflammatory response.22
TRANSFAC® (TRANScription FACtor) database positional weight matrices in promoters. Shown are maps of predicted binding sites of selected TRANSFAC® (Biobase GmbH, Wolfenbuettel, Germany) positional weight matrices in promoters of genes down-regulated in response to ischemia. TSS = transcription start site.
TRANSFAC® (TRANScription FACtor) database positional weight matrices in promoters. Shown are maps of predicted binding sites of selected TRANSFAC® (Biobase GmbH, Wolfenbuettel, Germany) positional weight matrices in promoters of genes down-regulated in response to ischemia. TSS = transcription start site.
ddPCR Validation and Technical Replication.
All four genes validated with ddPCR (MGST1, ITLN1, FOS, and HBB) showed a similar expression profile to their messenger RNA RNA-seq transcript when normalized to Gapdh. There was a high degree of correlation between the expression values obtained from RNA-seq and the normalized concentration from the ddPCR validation (r2 = 0.83 for baseline values, r2 = 0.76 for postischemia, and r2 = 0.85 for the fold changes). Normalized concentration and log2 fold changes between postischemia and baseline were ITLN1 (mean baseline/post 3,251/14, SD baseline/post 5,068/10, fold change −7.8), MGST1 (mean baseline/post 173/96, SD baseline/post 158/48, fold change −0.9), FOS (mean baseline/post 285/722, SD baseline/post 318/1,062, fold change 1.3), and HBB (mean baseline/post 2,672/3,254, SD baseline/post 2,981/2,579, fold change 0.3). Figure 7 clearly shows the high baseline expression and very low postischemia expression for ITLN1.
Digital droplet polymerase chain reaction validation. Fluorescence amplitude plots for carboxyfluorescein (FAM)–labeled droplets in a representative (A) baseline and (B) postischemia sample for the intelectin 1 (ITLN1) gene. Each dot represents a single droplet or event. The event number is displayed on the horizontal axis, whereas the fluorescence amplitude is shown on the vertical axis. The blue dots represent high fluorescence intensity (positive droplets), whereas the black dots represent low intensity (negative droplets).
Digital droplet polymerase chain reaction validation. Fluorescence amplitude plots for carboxyfluorescein (FAM)–labeled droplets in a representative (A) baseline and (B) postischemia sample for the intelectin 1 (ITLN1) gene. Each dot represents a single droplet or event. The event number is displayed on the horizontal axis, whereas the fluorescence amplitude is shown on the vertical axis. The blue dots represent high fluorescence intensity (positive droplets), whereas the black dots represent low intensity (negative droplets).
The results of our technical replication show that for the genes for which we performed the ddPCR, the fold changes show a similar magnitude and the same directionality, and the coefficient of variation is 79%. Likewise, when we compare all genes in both samples, the average coefficient of variation is 58% (baseline) and 51% (postischemia), and when we restrict to those with an RPKM 5 or greater, the coefficient of variation is 49 and 42%, respectively.
Discussion
Our comprehensive analysis of whole-genome RNA-seq data on human LV myocardium identifies the genes and pathways that are significantly and strongly regulated in cold cardioplegia-induced myocardial ischemia. Our results comprehensively delineate gene regulation at the transcriptional level in patients exposed to CPB and suggest that messenger RNA isoform regulation plays a critical role in the cardiac response to cold cardioplegia-induced ischemia. However, of the identified pathways that are activated in response to this myocardial injury, it is still not clear what part of the transcriptional response to ischemia is a protective mechanism versus a reflection of genetic activation of maladaptive injury response.
Hallmarks of the transcriptional response to ischemia included changes in genes and pathways associated with inflammation and the oxidative response to stress, which themselves can mediate myocardial injury. Nevertheless, the most profoundly differentially expressed genes have yet to be directly associated with human ischemic myocardial injury.
Down-regulation
The three most down-regulated genes, ITLN1, proteoglycan 4, and secretory leukocyte peptidase inhibitor, are significantly down-regulated by more than three-fold (P < 0.01).
Omentin/intelectin is a novel adipokine discovered in 2005, codified by two genes, with omentin 1 (ITLN1) as the primary circulating form. This protein stimulates insulin-mediated glucose transport in human adipocytes and is expressed in epicardial and omental adipose tissue. Its role in cardiovascular disease has been well described. Specifically, the actions of omentin are responsible for increasing endothelium-dependent vasodilation (rat ex vivo),23 decreasing vascular inflammation (human in vitro),24 and decreasing extracellular migration and angiogenesis (human in vitro).25 Hence, omentin appears to be a “protective adipokine,” with antiinflammatory properties as well as reducing endothelial dysfunction.26
Regulation of omentin in an acute setting has been studied in mice, showing that omentin promotes endothelial cell function and the revascularization process under conditions of ischemic stress, likely through adenosine monophosphate–activated protein kinase and Akt.27
The highly significant and consistent down-regulation of ITLN1 in myocardial tissue after exposure to cold cardioplegia-induced ischemia suggests that ITLN1 also plays an important role in myocardial injury. In addition, expression in nonischemic right atrial tissue is very similar to baseline values in ventricular myocardium (table 2) further supporting that ITLN1 expression in nonischemic cardiac tissue is high. As with ITLN1, secretory leukocyte peptidase inhibitor is also expressed in omental adipose tissue.28 This gene encodes a secreted inhibitor, which protects epithelial tissues from serine proteases. Serine protease inhibitors, such as the withdrawn Aprotinin, have a long clinical history in cardiac surgery and are attractive agents in cardiac surgery to protect the myocardium from ischemia–reperfusion injury.29,30 In summary, this suggests that ITLN1 and secretory leukocyte peptidase inhibitor, two genes without prior association with ischemic injury in human myocardium, are the targets of interest in modulation of ischemia and warrant further investigation.
Proteoglycan 4 produces a protein synthesized by chondrocytes at the surface of articular cartilage and synovial lining cells although expression in human and mouse cardiac tissue has been described.31 Similar to ITLN1, proteoglycan 4 expression at baseline is also comparable with atrial expression levels. The reason for profound down-regulation remains to be elucidated because the only link to cardiac disease is that a homozygous deletion in proteoglycan 4 causes camptodactyly-arthropathy-coxa vara-pericarditis syndrome.
Up-regulation
In general, up-regulation was not as profound in fold change or in number of genes affected as down-regulation. Nevertheless, there are several very intriguing targets, namely FOS and the hemoglobin genes HBA1 and HBB (fold change >2, P < 0.01).
The gene FOS is part of the FOS gene family, whose proteins can dimerize with proteins of the JUN family, thereby forming the transcription factor complex activator protein 1.32 These proteins are regulators of cell proliferation, differentiation, transformation, and even apoptotic cell death.33 Specifically in cardiac tissue in rats, up-regulation of c-Fos by tumor necrosis factor-α has been shown to increase the transcription of transforming growth factor β1, a key cytokine that promotes postmyocardial injury cardiac fibrosis.34 This up-regulation of FOS goes along with what we observed in our patients and is consistent with an activation of inflammatory mediators after CPB that ultimately can result in postischemia myocardial dysfunction.
We also show for the first time the profound up-regulation of HBA1 and HBB in the LV in response to cold cardioplegia-induced ischemia. The gene expression atlas shows an increased expression of HBA1 and HBB in heart tissue as compared with mean expression level for all samples,35 yet to our knowledge the response to ischemia in human LV myocardium has not been examined. Arab et al.36 showed a rapid increase in hemoglobin expression in right ventricular myocardium in infants and children with congenital heart disease in response to ischemia, with similar findings in isolated cardiomyocytes during simulated ischemia–reperfusion injury. Likewise, studies showed increased hemoglobin messenger RNA expression in neurons in embryonic and adult brains in the setting of transient and permanent cerebral artery occlusion in rats.37 Both findings are indicative of nonhematologic origins of hemoglobin gene transcription. A possible mechanism is the myocardium’s attempt to minimize nitrosative stress by sequestration of excess nitric oxide generated during ischemia, which is supported by our GO analysis showing the response to nitrosative stress as the top down-regulated category (fig. 4). This would imply a protective mechanism.
Functional Analysis
Our GO and network pathway analyses identified a number of pathways that are related to oxidative stress, defined as the imbalance between the production of reactive oxygen species and the detoxification of reactive intermediates. Although severe oxidative stress can trigger apoptosis and necrosis and is involved in cardiovascular diseases such as atherosclerosis, reactive oxygen species are not necessarily damaging as they are an integral part of the immune systems response to pathogens and part of normal cell signaling.
Complement activation is a key mediator of ischemia-induced myocardial tissue damage. Complement inhibitors have been applied successfully in animal models in reducing tissue damage from ischemia–reperfusion injury, yet anticomplement therapy in human myocardial ischemia–reperfusion injury has yielded inconsistent results.38 A direct link exists between the activation of the complement system and the regulation of the immune system and past evidence points to an increase in both in response to ischemia–reperfusion.39 However, the significant down-regulation of the humoral immune response concomitantly with both the alternative and classical complement pathway in our cohort would suggest that ischemia, in the absence of reperfusion, may not lead to an immediate increase in these pathways as previously described. This is supported by very low complement deposition (C3d) in brain-dead donor hearts that had ischemic injury without reperfusion,40 by the lag in C3 deposition compared with the loss of plasma membrane integrity during early myocardial ischemia,41 and by the initial decrease in serum C3 levels in pediatric patients undergoing cardiac surgery.42
Limitations
The focus of this study was the transcriptional response of the human LV myocardium to cold cardioplegia-induced ischemia. We used CPB in patients undergoing cardiac surgery as the model, as exposing normal human myocardium to ischemia in vivo is not feasible. In addition, our model of cold cardioplegia-induced myocardial ischemia provides important contributions to the field as many inconsistencies in prior studies using human myocardium stems from the fact that the most experiments are inherently uncontrolled using tissue that becomes available only as a result of death or cardiac transplantation. Our prospective study has several features that are novel: (1) every sample is collected from the exact same physical location on the LV apex by the same surgeon at the same time point in the operation, reducing intersample variability; (2) the time between sample collection and aliquoting into RNA preservative is approximately 1 min, resulting in extremely high-quality RNA; (3) the tissue sample we collect is from LV tissue (as opposed to atrial or right ventricular tissue used in prior studies) in nonfailing hearts (unlike LV assist device cores from patients with end-stage heart failure). Although the hearts do have concentric hypertrophy from long-standing aortic stenosis, the surgical procedures are all elective, and risk factors are optimally medically managed in preparation of the nonurgent operative procedure and adjusted for with multivariable modeling.
At the same time, our experimental design of ischemia induced by cold blood cardioplegia is not exact equivalent of in vivo ischemia, and our results might underestimate the true extent of ischemic injury. The down-regulation observed in the majority of our samples could also be attributed to hypothermia-induced suppression of transcription, and further work animal studies and in vitro experiments will be necessary to discern biologically relevant differences. In addition, lack of validation at the protein level, secondary to the small amount of tissue obtained, is also a limitation.
Conclusion
We have shown that the human LV exhibits significant yet variable changes in gene expression in response to cold cardioplegia-induced ischemia during CPB. Investigation of this process on a molecular level by examining whole-genome transcriptional activity provides great insight into the pathophysiology of ventricular ischemia and, thus, may help guide efforts to reduce myocardial damage during surgery as well as to identify interventional targets for pharmacological therapy and biomarkers of injury in the future.
Acknowledgments
The authors thank the outstanding contributory efforts of the Coronary Artery Bypass Graft Genomics research staff: James Gosnell, R.N.; Kujtim Bodinaku, M.D.; Svetlana Gorbatov, M.P.H. (all part of Brigham and Women’s Hospital, Department of Anesthesiology, Perioperative and Pain Medicine, Harvard Medical School, Boston, Massachusetts). The authors also thank all study subjects who participated in the TRANSCRIBE Program.
Supported by the National Institutes of Health (Bethesda, Maryland; grant no. R01HL118266 to Dr. Muehlschlegel), the American Heart Association (Dallas, Texas), The Harvard Catalyst (Boston, Massachusetts), The Watkins Cardiovascular Discovery Award (Brigham and Women’s Hospital, Boston, Massachusetts), the Scholars in Clinical Science Program (Harvard Medical School, Boston, Massachusetts), and a Mercedes Concepcion Faculty Development Fellowship and Scholar in Translational Anesthesia Research grant from the Department of Anesthesiology, Perioperative and Pain Medicine (Brigham and Women’s Hospital, Harvard Medical School, Boston, Massachusetts).
Competing Interests
The authors declare no competing interests.