Transcriptomic Landscape of Circulating Extracellular Vesicles in Heart Transplant Ischemia–Reperfusion

1. Introduction

Heart transplantation (HTx) remains the ultimate therapeutic option for patients with end-stage heart failure. Although outcomes have improved due to advances in immunosuppressive care and organ preservation methods, the short- and long-term survival of the allograft is negatively affected by factors such as ischemia–reperfusion injury (IRI). This inevitable complication after HTx is induced by periods of cold and warm ischemia during preservation and surgical procedure, followed by reperfusion and concomitant reoxygenation after implantation, altogether eventually causing exacerbation of cellular damage and primary graft dysfunction (PGD) [1,2]. Therefore, studies and strategies for the prevention of this phenomenon are of high clinical significance.

Extracellular vesicles (EVs) are a heterogeneous group of nanosized membranous vesicles secreted by all cell types. EVs play a role in cell-to-cell communication in various physiological processes as they carry, e.g., proteins, carbohydrates, lipids, messenger RNAs (mRNAs), and microRNAs that can interact with target cells [3]. In the regulation of inflammation, they have now been recognized as potential biomarkers and immune mediators in organ transplantation [4].

Ischemia–reperfusion-induced EV release may worsen myocardial injury by promoting inflammation [5]. However, their role in IRI after cardiac transplantation remains poorly understood. We aimed to characterize the transcriptomic landscape of plasma-derived EVs during IRI in human HTx using RNA sequencing. We compared gene expression changes in EVs derived from the plasma of heart transplant recipients before and at 12 h after reperfusion. Furthermore, we correlated clinical data with the plasma EV transcriptome to identify EV transcripts that are strongly linked to IRI-related myocardial injury and the development of PGD.

2. Materials and Methods

2.1. Study Design

This study was designed to characterize the recipient plasma EV transcriptome during ischemia–reperfusion in human HTx. In addition, we integrated the clinical data of the recipients to identify significant correlations between the EV transcriptome and occurrences of myocardial injury and PGD at 12 h after reperfusion. The patient cohort was initially designed for a clinical trial on donor simvastatin treatment [6,7] and was conducted as a prospective, double-blind, randomized single-center study at the Helsinki University Hospital between 2010 and 2016. The study was approved by the Ethical Committee Board (permission number 358/13/03/02/2009), and it was conducted in accordance with the principles of the Declaration of Helsinki. The blood samples used in this study were collected at two timepoints: immediately before reperfusion (0 h) and at 12 h (12 h) after reperfusion. Using next-generation sequencing technology, we profiled the whole EV transcriptome in an unbiased approach to thoroughly understand expression changes within the EV transcriptomic landscape. The classification of PGD grades was made within 24 h after the transplantation, in accordance with the ISHLT guidelines [8]. For transcriptomic profiling and subsequent correlation analyses with clinical data, we initially used paired plasma samples from 41 recipients. Following stringent sequencing data control, the finally analyzed dataset consisted of paired plasma samples from 38 recipients.

2.2. Sample Collection and Processing

Blood samples from 41 heart transplant recipients were collected into lithium heparin BD Vacutainer blood tubes (Becton Dickinson, Franklin Lakes, NJ, USA) immediately before reperfusion (0 h) and at 12 h after reperfusion (12 h) during the HTx procedure. Blood was centrifuged at 1600× g for 10 min at RT. Plasma was transferred to cryotubes (Nunc Thermo Scientific, Waltham, MA, USA) and stored at −80 °C. Upon analysis, 900 μL of thawed plasma was diluted at a 1:1 ratio with 0.22 μm filtered PBS and filtered with 0.8 μm Millex AA syringe filters (Millipore Sigma, Burlington, MA, USA). To isolate EV, diluted plasma was ultracentrifuged using SW 55 Ti rotor and Polypropylene Thickwall tubes (Beckmann-Coulter, Brea, CA, USA) at 40,000 rpm at +4 °C for 90 min. EV pellets were washed with 2 mL of PBS and ultracentrifuged, resuspended in 100 μL of filtered PBS, and stored in DNALobind tubes (Eppendorf, Hamburg, Germany). To eliminate heparin residues, heparinase I (New England Biolabs, Ipswich, MA, USA) was added to the EV resuspension (0.1 U/μL) and incubated at 37 °C for 1 h. RNA was isolated from EVs by using a Qiagen miRNeasy micro kit (Qiagen, Hilden, Germany). The quality and concentration of EV RNA were measured using the Bioanalyzer Pico kit (Agilent, Santa Clara, CA, USA).

2.3. Transmission Electron Microscopy and Nanoparticle Tracking Analysis

To examine the presence and morphology of isolated EV, representative samples were prepared for transmission electron microscopy (0 h, n = 6 and 12 h, n = 5) [9]. To analyze the concentration and size distribution of isolated EV, 20 μL of the representative eluted EV samples (0 h, n = 5, and 12 h, n = 8) were diluted in Dulbecco’s PBS and analyzed on NanoSight LM14C (Nanosight, Salisbury, UK) and Zetaview PMX-120 model (ParticleMetrix, Ammersee, Germany) equipped with a blue laser (404 nm and 488 nm, respectively) and a high sensitivity CMOS camera. EV size and concentration data were obtained with the NanoSight instrument by recording 90 s of 3 videos, and EV size and concentration data analyzed with Zetaview were captured by scanning 11 positions with 30 frames and 12.41 μS/cm sensed conductivity. To determine the concentration and the size of isolated EV, data obtained were analyzed using NTA software 3.0 and ZetaView software (Version 8.05.12 SP2).

2.4. RNA Sequencing

EV RNA samples from 17 recipients were sequenced on Illumina HiSeq 2500 and samples from 21 recipients were sequenced on Illumina NovaSeq 6000 (Illumina, San Diego, CA, USA). Pair-ended cDNA libraries were prepared using SMARTseqv4 Ultra Low Input RNA kit (Takara Bio Inc., Kusatsu, Shiga, Japan). One hundred picograms of amplified cDNA were tagmented and indexed using the Nextera XT DNA library prep kit (Illumina). Finally, cDNA libraries were sequenced using a sequencing v4 kit (Illumina) to obtain at least 15 million reads on the HiSeq and 25 million reads on the NovaSeq platform.

2.5. Bioinformatics

RNAseq raw data were processed with the Trimmomatic software (version 0.39) [10] and were aligned to the Human Genome (GRCh38) with the STAR software (version 2.7.10a) [11]. Gene expression read counts were calculated by using the Rsubread software (version 2.16.0) [12]. Sequencing data were normalized with the TMM normalization method prior to the statistical analysis to avoid differences between sequencing platforms and batches. Raw counts, normalized CPM, and RPKM values are given in . Outliers were identified by generating a principal component analysis (PCA) plot from counts per million (CPM) values for each gene in the sample. Analyses such as PCA and unsupervised hierarchical clustering were performed with the ropls (version 1.34.0) [13] and factoextra R packages (version 1.0.7). Differential gene expression testing was performed using edgeR (version 4.0.1) [14]. Protein-coding genes presenting absolute log-change greater than 1.5 and q-value (FDR-adjusted p-value) less than 0.05 were considered to be significant. In addition, the R package CombiROC (version 0.3.4) and pROC (version 1.18.5) were used to calculate the power of discrimination for 1317 differentially expressed genes (DEGs). Heatmaps of dDEGs were visualized by the pheatmap R package (version 1.0.12) [15]. Functional annotations of DEGs were investigated with the clusterProfiler (version 4.4.4) using its functions enrichGO and enrichKEGG. Enriched terms with a q-value (FDR-adjusted p-value) of less than 0.05 were considered significant. Outputs of enrichment analyses were visualized using the dotplot function in the clusterProfiler R package [16]. We summarized and reduced the redundant enriched gene ontology terms by using the Revigo-based [17] R package RRVGO (version 1.10.0) [18]. Summarized terms were visualized using the scatterplot function included in the same package. To examine the protein network of the genes, we used the STRING analysis tool (https://string-db.org/, version 12.0, last accessed on 2 October 2023) [19]. Furthermore, we used Web-based cell-specific enrichment analysis (WebCSEA, https://bioinfo.uth.edu/webcsea/, last accessed on 6 October 2023) [20] to elucidate the cell types that may be associated with the significant 1317 protein-coding genes found at 12 h post-reperfusion when compared to 0 h. This tool enables the enrichment of systematic gene sets for human tissue cell-type expression signatures. The most significant cell types were identified and visualized based on the −log10 fold change combined with the p-value. For correlation analysis of high-sensitivity plasma troponin I level (TnI), PGD, and AR with EV transcriptome, we used the Wilcox rank sum test. The final set of significant genes was determined with a p-value cutoff of less than 0.005.

3. Results

3.1. Recipient Characteristics

Plasma EV RNA was extracted from 41 heart transplant recipients, from samples taken immediately before and at 12 h after reperfusion and then subjected to next-generation RNA sequencing. Three recipients were excluded due to the low RNA quality of their samples. Of the remaining 38 recipients, the median age was 60 (range 27–67), and 9 (23.7%) of them were female. The most common cause for HTx was dilated cardiomyopathy (46.2%). A total of 26.2% of patients were bridged to transplantation with mechanical circulatory support . A CONSORT flowchart describes the study flow (Figure 1).

Figure 1. CONSORT-type workflow chart of the study. Plasma EV transcriptomes before reperfusion and at 12 h after reperfusion of the heart transplant were compared with each other. Only the EV transcriptome of plasma samples taken at 12 h after reperfusion was correlated with the clinical data of the heart transplant recipient. Figure created with BioRender.com (accessed on 14 November 2023).

3.2. EV Characterization

Plasma EVs showed a typical cup-shaped morphology in the electron microscopic images (0 h, n = 6 and 12 h, n = 5) (Figure 2A,B). Nanoparticle tracking analysis revealed that the size distribution of isolated EVs from 0 h (n = 5) and 12 h (n = 8) samples ranged within approximately 50 nm to 1355 nm and 50 nm to 1095 nm, respectively. The plasma EV concentration was similar at both timepoints (Figure 2C,D).

Figure 2. Characterization of isolated EVs. Transmission electron microscopy and nanoparticle tracking analysis were performed to confirm and examine the presence, morphology, size distribution, and concentration of recipient plasma EVs. (A) Transmission electron microscopy of plasma EVs at 0 h and (B) at 12 h after reperfusion of the heart. Arrows indicate the typical cup-shaped morphology of EVs. (C) Concentration (particles/mL) and (D) size distribution data were obtained from nanoparticle measurement equipment (NanoSight and ZetaView); the dashed line in each plot was plotted by geom_smooth function in R, using the LOESS method to represent the trend of concentration values throughout the samples according to the size of the samples.

3.3. Reperfusion Altered mRNA Expression Profiles of Recipient Plasma EV

Unsupervised PCA showed distinct plasma EV transcriptomic profiles at 12 h after heart transplant reperfusion compared to 0 h before reperfusion (Figure 3A,B). Raw counts, normalized CPM, and RPKM values are given in . DEG analysis, with a cut-off of log-fold change greater than 1.5 and a q-value less than 0.05, revealed a total of 1317 (3.4%) protein-coding genes. Among these, 25.6% of the genes were found to be upregulated, while 74.4% were downregulated (Figure 3C and ). We additionally conducted ROC analysis to assess different powers of discrimination for DEGs . When applied minimum specificity of 80% and sensitivity of 80% for high-power discrimination, 27% of DEGs (357 out of 1317) at 12 h passed the cutoff criteria . These genes exhibited similar biological functions to those we observed in pathway analyses when using 1317 DEGs .

Figure 3. Differentially expressed genes at 12 h post-reperfusion compared to pre-reperfusion. (A) Unsupervised principal component analysis showed that different timepoint (blue, 0 h; red, 12 h) samples have distinct EV transcriptomic profiles. (B) Hierarchical clustering analysis confirmed clear differences in the gene expression patterns between 0 h and 12 h groups. (red, upregulated; green, downregulated) (C) The number of differentially expressed genes at 12 h after reperfusion is illustrated in a pie chart format.

3.4. Biological Functions Related to the Differentially Expressed Genes after Reperfusion

To gain a better understanding of the impact of heart transplant reperfusion injury on the plasma EV transcriptome, we investigated the 1317 DEGs by enriching biological annotations. Firstly, we performed gene ontology enrichment over-representation analysis to explore relevant biological processes. Upregulated genes at 12 h after reperfusion were mainly linked to cell division, leukocyte-mediated immunity, and metabolic processes (Figure 4A), whereas downregulated genes were related to extracellular structure, tissue development, and cell proliferation (Figure 4B). To refine and reduce possible redundant enriched gene ontology terms, we summarized those terms based on a semantic similarity score by the Relevance method [21]. The calculated scores thus yielded more precise terms related to biological processes. We found that upregulated genes were related to leukocyte-mediated immunity, acute-phase response, regulation of coagulation, lipid metabolic processes, and mitosis , while downregulated genes were linked to complement activation, cell chemotaxis, wound healing, response to oxygen levels, organ development, and collagen processes . Next, we conducted a KEGG pathway enrichment analysis to explore possible biological pathways related to the DEGs. In alignment with the findings in the gene ontology analysis, upregulated genes were associated with cytokine receptor interaction, cell cycle, natural-killer-cell-mediated cytotoxicity, retinol metabolism, antigen processing and presentation, allograft rejection, and complement and coagulation cascades (Figure 4C and ). Downregulated genes were related to cytokine interaction, cell signaling pathways (e.g., PI3K-Akt, AGE-RAGE, Rap1, Hippo, and TGF-β), focal adhesion, complement and coagulation cascades, and ECM–receptor interaction (Figure 4D and ). Lastly, we performed WebCSEA analysis to get a better insight into the possible origin of EV transcripts, based on the tissue cell-specific signatures. The results suggested that DEGs at 12 h after reperfusion showed various immune cell and organ tissue cell-type signatures in the overall organ system (Figure 4E). Since we analyzed plasma samples from heart transplant recipients, we specifically investigated whether the DEGs were related to heart tissue. The results suggested that DEGs at 12 h after reperfusion showed cell-type-specific signatures associated with fibroblasts, stromal cells, smooth muscle cells, pericytes, mesenchymal stem cells, epicardial adipocytes, endocardial cells, endothelial cells, and Schwann cells in the heart, based on the Tabula Sapiens database [22] (Figure 4F and ). Although Tabula Sapiens contains hepatocytes in the heart tissue cell signatures, we excluded hepatocytes from our results due to their irrelevance to heart histology. In addition, with detailed analysis by separating upregulated and downregulated genes, we observed heart-related immune cells were mainly upregulated and cell structure component genes were downregulated .

Figure 4. The biological relevance of differentially expressed genes found at 12 h post-reperfusion using pathway enrichment analyses. The 1317 differentially expressed genes (DEGs) were analyzed to explore their biological functions by enriching relevant terms with different databases. Focusing on biological processes, gene ontology enrichment revealed terms associated with (A) upregulated DEGs and (B) downregulated DEGs. Enriched pathways based on the KEGG database suggested relevant pathways of (C) upregulated and (D) downregulated DEGs. (E) WebCSEA results to visualize overall organ system-specific cell signatures, and (F) WebCSEA results to visualize heart tissue cell-specific signatures. Since our patient cohort consists exclusively of cardiac transplant recipients, we focused primarily on heart-related genes.

3.5. Post-Reperfusion EVs Carry More Adaptive Immunity-Related Genes but Less Extracellular Matrix Component Protein-Coding Genes

Furthermore, we hypothesized that the genes with higher log-fold changes may represent the most significant biological functions related to those genes. Therefore, we filtered the DEGs based on the log fold-change and selected 50 upregulated and 50 downregulated genes for STRING analysis, rather than focusing on genes related to the most significantly enriched GO and pathway terms to avoid biased interpretation of established gene networks. STRING analysis was performed to explore the protein–protein network of DEGs, as well as relevant biological annotations of these genes. The most upregulated 50 genes were associated with retinol metabolism, complement system, and graft-versus-host disease, sharing a gene set related to allograft rejection (Figure 5A). Downregulated genes showed well-established protein–protein networks, and the abundantly enriched biological annotations were related to extracellular matrix organization, including collagen-containing extracellular matrix, and organ development (Figure 5B).

Figure 5. Protein network analysis of the top 50 most significantly upregulated and downregulated genes at 12 h post-reperfusion. Probable protein–protein networks and biological functions of the most significant (A) 50 upregulated and (B) 50 downregulated genes were analyzed based on the STRING online database. Edges represent protein–protein associations based on curated databases, text mining, and experimentally determined interactions. Edges are connected with lines whose thickness represents the confidence of the network association.

3.6. Correlation of EV Transcriptome with Cardiomyocyte Injury

Our previous study has reported that plasma high-sensitivity troponins are highest at 12 h in recipients who had severe left ventricular PGD [23]. Therefore, we investigated all EV transcripts expressed at 12 h to identify significant correlations between the EV transcriptome and cardiomyocyte injury, defined by 12 h plasma level of high-sensitivity TnI. We found that 48 protein-coding genes were significantly correlated with myocardial injury (Figure 6A and , and they exhibited correlations in their expression with each other (Figure 6B). In addition, of the 48 genes correlated to TnI plasma levels, 8 were also found among the DEG after reperfusion, suggesting their expression patterns are significantly affected after reperfusion .

Figure 6. EV transcripts that showed significant correlations with clinical outcomes. (A) Forty-eight genes showing correlations between the EV transcriptome and the myocardial injury marker; high-sensitivity troponin I, measured at 12 h after reperfusion (p < 0.005). (B) Correlations among the expression levels of these 48 genes (p < 0.005). (C) Four genes significantly correlated with both TnI and PGD (p-value < 0.005). (D) Correlations of expression levels among the four genes. (E) Expression levels of the four genes (p < 0.005). p-values in B and D are denoted with asterisk marks, * p < 0.05, ** p < 0.005, *** p < 0.001.

3.7. Association of EV Transcriptome and Outcomes after Heart Transplantation

We next compared the plasma EV transcriptome of recipients without PGD (66%) with those who developed PGD of any grade (34%). The result revealed that 3194 protein-coding genes were related to the higher incidence of PGD when applying a p-value of less than 0.05. Generally, these genes that showed correlations with PGD (p < 0.05, n = 3194) were associated with lipid metabolism, acute-phase response, glucuronidation, TNF signaling pathway, and complement and coagulation cascades . We found 146 protein-coding genes when using a more stringent p-value cut-off (p-value < 0.005) and 30 protein-coding genes when applying a p-value of less than 0.001. Most of these 30 genes were upregulated compared to recipients without PGD . To investigate whether the gene expression of these 146 genes was affected by reperfusion, we compared the genes with 1317 DEGs at 12 h after reperfusion. We found 25 genes that were also present in the DEGs after reperfusion, and their biological functions were related to retinol metabolism, amyloidosis, bone marrow cells, plasma cells, liver, and digestive gland . However, no significant differences in donor characteristics were observed between the recipient PGD groups . Given that cardiomyocyte damage is a prominent manifestation of PGD after HTx, we cross-checked our correlation analyses, and identified four protein-coding genes (ITPKA, DDIT4L, CD19, and CYP4A11) that were consistently observed (Figure 6C and . These genes showed correlations in expression among each other (Figure 6D) and three of them showed upregulation in the recipient group with PGD after HTx (Figure 6E). Although reperfusion injury may have a greater impact on the early outcome (i.e., PGD) than on long-term outcomes such as acute rejection, we investigated the potential relationship between genes related to reperfusion injury and biopsy-proven acute rejection. We examined differences in the EV transcriptome at 12 h based on the recipients’ IV-treatment history within 30 days after HTx and biopsy-proven AR diagnosis history within the first year after HTx. Among the 38 recipients, 9 received IV treatment within 30 days after HTx, while 27 did not (2 not available). Additionally, 24 recipients developed any biopsy-proven AR within the first year after HTx, and 14 recipients remained stable. We found 100 genes that showed significantly different expression levels based on the history of IV treatment within 30 days and another 100 genes associated with biopsy-proven AR within the first year, respectively (p < 0.01). Moreover, 31 and 42 genes were found in each analysis when applying a p-value cutoff of 0.005 . However, we did not find any significantly enriched pathways and gene ontology terms using biopsy-proven AR prediction genes.

4. Discussion

In this study, we provide a detailed analysis of the transcriptomic landscape of plasma-derived EVs following IRI in human HTx. While transcriptomic signatures of IRI in experimental settings have been well established using bulk RNAseq, the whole transcriptomic profile of plasma EVs in clinical settings has been unexplored. We uncovered 1317 EV transcripts that were significantly altered during heart transplant IRI, including upregulation of genes related to immune activation and metabolism mediated by cytochrome P450 (CYP) and downregulation of genes associated with cell proliferation, extracellular matrix, smooth muscle contraction, and TGF-β signaling pathway. In addition, EV transcripts were linked to cardiomyocyte injury and PGD as 48 protein-coding EV transcripts were associated with higher levels of the cardiomyocyte injury marker TnI and 30 EV genes with the development of PGD after HTx. Importantly, four transcripts—ITPKA, DDIT4L, CD19, and CYP4A11—were closely associated with both cardiomyocyte injury and PGD, indicating that these protein-coding genes may have important roles in regulating heart transplant IRI.

A recent systematic review of bulk RNA analyses of myocardial tissue suggests that the transcriptomic hallmarks of IRI include increased mRNA expressions related to response to stress, cell proliferation, inflammatory response, and cell pathway (i.e., TNF, NF-kB, IL-17, MAPK, TLR, and NOD-like signaling pathways), while downregulated genes were associated with cytoskeletal structures, cell development, cell survival, ion channels, and cAMP signaling pathway [24]. Our findings showed that recipient plasma EV transcripts were associated with upregulated genes related to inflammation, T-cell chemotaxis, and cytokine interactions, while downregulated genes were linked to ion transport and extracellular matrix . Interestingly, we found that circulating EVs contain gene transcripts related to allograft rejection already at 12 h after reperfusion, along with transcripts associated with cell proliferation and metabolism. Moreover, WebCSEA analysis focused on the heart tissue database suggested that DEGs at 12 h after reperfusion were related to a broad range of tissue structural components, such as fibroblasts, stromal cells, smooth muscle cells, pericytes, epicardial adipocytes, and others related to heart tissue. Nevertheless, our pathway analyses and tissue-specific cell type enrichment analyses suggested that circulating EVs after reperfusion may be involved in or reflect the regulation of cellular structure, hence affecting tissue remodeling and repair in the cardiovascular system after reperfusion. Collectively, our results may help to understand how the cellular structure disruption and inflammation after IRI contribute to the development of PGD after HTx.

Upon reperfusion of the heart transplant, oxygen radicals may cause tissue damage by interacting with polyunsaturated fatty acids and by the formation of lipid peroxides and hydroperoxides. Oxygen free radicals induce the release of platelet-activating factors by endothelial cells and, thereby, aggravate neutrophil-mediated immune responses. Consequently, elevated proinflammatory signals and activation of downstream signaling pathways may ultimately result in the development of parenchymal fibrosis and decreased cardiac function. Given that, we speculate that lipid metabolism and lipid-metabolism-induced stress are key factors for understanding the mechanism of IRI. For example, polyunsaturated fatty acids such as arachidonic acid and linoleic acid are well-known to participate in the induction of endothelial dysfunction, vascular tone, calcium mobilization, and oxidative stress [24,25,26,27,28,29]. Especially, arachidonic acid is an essential fatty acid and precursor of prostaglandins, thromboxanes, and leukotrienes, which are important players in endothelial dysfunction [30]. Of the four genes that showed a moderate correlation with both cardiomyocyte injury and PGD, CYP4A11 has been shown to have a role in lipid metabolism, reactive-oxygen-species-induced lipid peroxidation, and inflammation in non-alcoholic fatty liver disease [31]. There is growing evidence suggesting that CYP enzymes play a pivotal role in the cardiovascular system [32]; for example, polymorphisms in CYP4A11 have associations with coronary artery diseases [33]. Furthermore, DDIT4L, another one of the four genes highlighted in the results, promotes cardiomyocyte cell death by inhibiting the mTOR signaling pathway under hypoxia/reoxygenation setting [34] and increases oxidized LDL-induced cytotoxicity [35], which may indicate the relation between oxidative stress and lipid metabolism. Furthermore, DDIT4L expression is reported in pathological cardiac hypertrophy, and overexpression of DDIT4L increases autophagy and causes mild systolic dysfunction [36].

Next, we focused on the calcium signaling pathway. During IRI, calcium overload can cause cardiac systolic dysfunction, because calcium ions play an essential role in cardiac action potential [37]. We found both ITPKA and CD19 may participate in the regulation of calcium transport in cardiomyocytes and the PI3K signaling pathway. Although CD19 is well known as a marker for B cells, we focused on the role of CD19 and B-cell receptor pathway cascades that participate in calcium homeostasis [38]. BCR-induced calcium transport has a connection with G-protein-coupled receptors, which are closely related to cardiovascular health [39], and enhanced BCR signaling under deficiency of CD19 is related to PI3K activity [40], and changes in myocardial B-cell population may affect myocardial growth and contractility [41]. Downregulation of calcium signaling pathway-related transcripts may result in impaired cardiac contractibility and remodeling due to dysregulation of excitation–contraction coupling [42]. Altogether, an imbalance in calcium homeostasis and activation of downstream signaling pathways may significantly affect the outcome of a heart transplant by altering the metabolism and result in mitochondrial dysfunction, inflammation, and reduced heart function.

One limitation of this study is the relatively small cohort size of 38 recipients. Additionally, validation of RNAseq data by qRT-PCR was not performed due to a limited volume of samples. Revalidation of the data in a larger cohort using both RNAseq and qPCR could further enhance our understanding of the mechanistic connection between EV and IRI and the role of EV in HTx. Moreover, we cannot differentiate transplantation-associated IRI from normal IRI, as well as alloimmune response reflected in circulating EVs after reperfusion. In addition, comparing blood samples from the coronary sinus and peripheral blood may provide enhanced specificity of EV analysis related to cardiac conditions. However, we did not consider this when initiating the original randomized control trial. Future studies using external cohorts and comparisons with samples taken from recipients before transplantation and during ischemia may enable us to differentiate the effects of ischemia-induced and reperfusion-induced injuries on the EV transcriptome.

5. Conclusions

In conclusion, we provide a comprehensive insight into the whole transcriptomic profile of plasma EVs from heart transplant recipients after reperfusion. In addition, we highlight four transcripts that EVs carry after reperfusion and suggest their clinical relevance both in cardiomyocyte injury and in PGD after HTx. Altogether, our findings suggest that EV may provide important information about heart transplant IRI. In addition to our analysis focused on protein-coding genes in plasma EVs, investigating the expression profiles of long non-coding RNAs, pseudogenes, and microRNAs may provide a deeper understanding of the EV biomarker capabilities or intricate roles of plasma EVs as gene expression modulators in IRI. Additional in-depth studies should be carried out to further elucidate the possible tissue sources of the circulating EVs and their roles in IRI and allograft damage. Ultimately, the profiling of EV transcriptome may be utilized as a novel approach for the development of non-invasive biomarkers for PGD after HTx.

References

1. Tang, Q.; Dong, C.; Sun, Q. Immune Response Associated with Ischemia and Reperfusion Injury during Organ Transplantation. Inflamm. Res. 2022, 71, 1463–1476. [Google Scholar] [CrossRef] [PubMed]
2. Fernández, A.R.; Sánchez-Tarjuelo, R.; Cravedi, P.; Ochando, J.; López-Hoyos, M. Review: Ischemia Reperfusion Injury—A Translational Perspective in Organ Transplantation. Int. J. Mol. Sci. 2020, 21, 8549. [Google Scholar] [CrossRef]
3. Valadi, H.; Ekström, K.; Bossios, A.; Sjöstrand, M.; Lee, J.J.; Lötvall, J.O. Exosome-Mediated Transfer of MRNAs and MicroRNAs Is a Novel Mechanism of Genetic Exchange between Cells. Nat. Cell Biol. 2007, 9, 654–659. [Google Scholar] [CrossRef]
4. Ravichandran, R.; Bansal, S.; Rahman, M.; Sureshbabu, A.; Sankpal, N.; Fleming, T.; Bharat, A.; Mohanakumar, T. Extracellular Vesicles Mediate Immune Responses to Tissue-Associated Self-Antigens: Role in Solid Organ Transplantations. Front. Immunol. 2022, 13, 861583. [Google Scholar] [CrossRef]
5. Ge, X.; Meng, Q.; Wei, L.; Liu, J.; Li, M.; Liang, X.; Lin, F.; Zhang, Y.; Li, Y.; Liu, Z.; et al. Myocardial Ischemia-reperfusion Induced Cardiac Extracellular Vesicles Harbour Proinflammatory Features and Aggravate Heart Injury. J. Extracell. Vesicles 2021, 10, e12072. [Google Scholar] [CrossRef] [PubMed]
6. Tuuminen, R.; Syrjälä, S.; Krebs, R.; Keränen, M.A.I.; Koli, K.; Abo-Ramadan, U.; Neuvonen, P.J.; Tikkanen, J.M.; Nykänen, A.I.; Lemström, K.B. Donor Simvastatin Treatment Abolishes Rat Cardiac Allograft Ischemia/Reperfusion Injury and Chronic Rejection Through Microvascular Protection. Circulation 2011, 124, 1138–1150. [Google Scholar] [CrossRef] [PubMed]
7. Nykänen, A.I.; Holmström, E.J.; Tuuminen, R.; Krebs, R.; Dhaygude, K.; Kankainen, M.; Jokinen, J.J.; Lommi, J.; Helanterä, I.; Räisänen-Sokolowski, A.; et al. Donor Simvastatin Treatment in Heart Transplantation. Circulation 2019, 140, 627–640. [Google Scholar] [CrossRef]
8. Kobashigawa, J.; Zuckermann, A.; Macdonald, P.; Leprince, P.; Esmailian, F.; Luu, M.; Mancini, D.; Patel, J.; Razi, R.; Reichenspurner, H.; et al. Report from a Consensus Conference on Primary Graft Dysfunction after Cardiac Transplantation. J. Heart Lung Transplant. 2014, 33, 327–340. [Google Scholar] [CrossRef] [PubMed]
9. Puhka, M.; Nordberg, M.-E.; Valkonen, S.; Rannikko, A.; Kallioniemi, O.; Siljander, P.; Af Hällström, T.M. KeepEX, a Simple Dilution Protocol for Improving Extracellular Vesicle Yields from Urine. Eur. J. Pharm. Sci. 2017, 98, 30–39. [Google Scholar] [CrossRef]
10. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A Flexible Trimmer for Illumina Sequence Data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef]
11. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast Universal RNA-Seq Aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef]
12. Liao, Y.; Smyth, G.K.; Shi, W. The Subread Aligner: Fast, Accurate and Scalable Read Mapping by Seed-and-Vote. Nucleic Acids Res. 2013, 41, e108. [Google Scholar] [CrossRef]
13. Thévenot, E.A.; Roux, A.; Xu, Y.; Ezan, E.; Junot, C. Analysis of the Human Adult Urinary Metabolome Variations with Age, Body Mass Index, and Gender by Implementing a Comprehensive Workflow for Univariate and OPLS Statistical Analyses. J. Proteome Res. 2015, 14, 3322–3335. [Google Scholar] [CrossRef] [PubMed]
14. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed]
15. Kolde, R. Pheatmap: Pretty Heatmaps. Available online: https://CRAN.R-project.org/package=pheatmap (accessed on 11 August 2022).
16. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS 2012, 16, 284–287. [Google Scholar] [CrossRef] [PubMed]
17. Supek, F.; Bošnjak, M.; Škunca, N.; Šmuc, T. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLoS ONE 2011, 6, e21800. [Google Scholar] [CrossRef] [PubMed]
18. Sayols, S. Rrvgo: A Bioconductor Package for Interpreting Lists of Gene Ontology Terms. MicroPubl. Biol. 2023, 2023. [Google Scholar] [CrossRef]
19. Szklarczyk, D.; Gable, A.L.; Lyon, D.; Junge, A.; Wyder, S.; Huerta-Cepas, J.; Simonovic, M.; Doncheva, N.T.; Morris, J.H.; Bork, P.; et al. STRING V11: Protein–Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-Wide Experimental Datasets. Nucleic Acids Res. 2019, 47, D607–D613. [Google Scholar] [CrossRef]
20. Dai, Y.; Hu, R.; Liu, A.; Cho, K.S.; Manuel, A.M.; Li, X.; Dong, X.; Jia, P.; Zhao, Z. WebCSEA: Web-Based Cell-Type-Specific Enrichment Analysis of Genes. Nucleic Acids Res. 2022, 50, W782–W790. [Google Scholar] [CrossRef] [PubMed]
21. Schlicker, A.; Domingues, F.S.; Rahnenführer, J.; Lengauer, T. A New Measure for Functional Similarity of Gene Products Based on Gene Ontology. BMC Bioinform. 2006, 7, 302. [Google Scholar] [CrossRef]
22. Jones, R.C.; Karkanias, J.; Krasnow, M.A.; Pisco, A.O.; Quake, S.R.; Salzman, J.; Yosef, N.; Bulthaup, B.; Brown, P.; Harper, W.; et al. The Tabula Sapiens: A Multiple-Organ, Single-Cell Transcriptomic Atlas of Humans. Science 2022, 376, eabl4896. [Google Scholar] [CrossRef]
23. Holmström, E.J.; Syrjälä, S.O.; Dhaygude, K.; Tuuminen, R.; Krebs, R.; Nykänen, A.; Lemström, K.B. Severe Primary Graft Dysfunction of the Heart Transplant Is Associated with Increased Plasma and Intragraft Proinflammatory Cytokine Expression. J. Heart Lung Transplant. 2023, 42, 807–818. [Google Scholar] [CrossRef]
24. Movahed, M.; Brockie, S.; Hong, J.; Fehlings, M.G. Transcriptomic Hallmarks of Ischemia-Reperfusion Injury. Cells 2021, 10, 1838. [Google Scholar] [CrossRef] [PubMed]
25. Verma, S.; Fedak, P.W.M.; Weisel, R.D.; Butany, J.; Rao, V.; Maitland, A.; Li, R.-K.; Dhillon, B.; Yau, T.M. Fundamentals of Reperfusion Injury for the Clinical Cardiologist. Circulation 2002, 105, 2332–2336. [Google Scholar] [CrossRef] [PubMed]
26. Ghosh, A.; Gao, L.; Thakur, A.; Siu, P.M.; Lai, C.W.K. Role of Free Fatty Acids in Endothelial Dysfunction. J. Biomed. Sci. 2017, 24, 50. [Google Scholar] [CrossRef] [PubMed]
27. Kim, F.; Tysseling, K.A.; Rice, J.; Pham, M.; Haji, L.; Gallis, B.M.; Baas, A.S.; Paramsothy, P.; Giachelli, C.M.; Corson, M.A.; et al. Free Fatty Acid Impairment of Nitric Oxide Production in Endothelial Cells Is Mediated by IKKβ. Arterioscler. Thromb. Vasc. Biol. 2005, 25, 989–994. [Google Scholar] [CrossRef]
28. Zhang, C.; He, M.; Ni, L.; He, K.; Su, K.; Deng, Y.; Li, Y.; Xia, H. The Role of Arachidonic Acid Metabolism in Myocardial Ischemia–Reperfusion Injury. Cell Biochem. Biophys. 2020, 78, 255–265. [Google Scholar] [CrossRef]
29. Chauhan, P.K.; Sowdhamini, R. Transcriptome Data Analysis of Primary Cardiomyopathies Reveals Perturbations in Arachidonic Acid Metabolism. Front. Cardiovasc. Med. 2023, 10, 1111019. [Google Scholar] [CrossRef]
30. Zhou, Y.; Khan, H.; Xiao, J.; Cheang, W.S. Effects of Arachidonic Acid Metabolites on Cardiovascular Health and Disease. Int. J. Mol. Sci. 2021, 22, 12029. [Google Scholar] [CrossRef]
31. Gao, H.; Cao, Y.; Xia, H.; Zhu, X.; Jin, Y. CYP4A11 Is Involved in the Development of Nonalcoholic Fatty Liver Disease via ROS-induced Lipid Peroxidation and Inflammation. Int. J. Mol. Med. 2020, 45, 1121–1129. [Google Scholar] [CrossRef] [PubMed]
32. Elbekai, R.H.; El-Kadi, A.O.S. Cytochrome P450 Enzymes: Central Players in Cardiovascular Health and Disease. Pharmacol. Ther. 2006, 112, 564–587. [Google Scholar] [CrossRef]
33. Sirotina, S.; Ponomarenko, I.; Kharchenko, A.; Bykanova, M.; Bocharova, A.; Vagaytseva, K.; Stepanov, V.; Churnosov, M.; Solodilova, M.; Polonikov, A. A Novel Polymorphism in the Promoter of the CYP4A11 Gene Is Associated with Susceptibility to Coronary Artery Disease. Dis. Markers 2018, 2018, 5812802. [Google Scholar] [CrossRef]
34. Hernández, G.; Lal, H.; Fidalgo, M.; Guerrero, A.; Zalvide, J.; Force, T.; Pombo, C.M. A Novel Cardioprotective P38-MAPK/MTOR Pathway. Exp. Cell Res. 2011, 317, 2938–2949. [Google Scholar] [CrossRef]
35. Cuaz-Pérolin, C.; Furman, C.; Larigauderie, G.; Legedz, L.; Lasselin, C.; Copin, C.; Jaye, M.; Searfoss, G.; Yu, K.T.; Duverger, N.; et al. REDD2 Gene Is Upregulated by Modified LDL or Hypoxia and Mediates Human Macrophage Cell Death. Arterioscler. Thromb. Vasc. Biol. 2004, 24, 1830–1835. [Google Scholar] [CrossRef]
36. Simonson, B.; Subramanya, V.; Chan, M.C.; Zhang, A.; Franchino, H.; Ottaviano, F.; Mishra, M.K.; Knight, A.C.; Hunt, D.; Ghiran, I.; et al. DDiT4L Promotes Autophagy and Inhibits Pathological Cardiac Hypertrophy in Response to Stress. Sci. Signal. 2017, 10, eaaf5967. [Google Scholar] [CrossRef] [PubMed]
37. Fearnley, C.J.; Roderick, H.L.; Bootman, M.D. Calcium Signaling in Cardiac Myocytes. Cold Spring Harb. Perspect. Biol. 2011, 3, a004242. [Google Scholar] [CrossRef] [PubMed]
38. Buhl, A.M.; Pleiman, C.M.; Rickert, R.C.; Cambier, J.C. Qualitative Regulation of B Cell Antigen Receptor Signaling by CD19: Selective Requirement for PI3-Kinase Activation, Inositol-1,4,5-Trisphosphate Production and Ca²⁺ Mobilization. J. Exp. Med. 1997, 186, 1897–1910. [Google Scholar] [CrossRef] [PubMed]
39. Li, Y.; Li, B.; Chen, W.-D.; Wang, Y.-D. Role of G-Protein Coupled Receptors in Cardiovascular Diseases. Front. Cardiovasc. Med. 2023, 10, 1130312. [Google Scholar] [CrossRef]
40. Xu, Y.; Fairfax, K.; Light, A.; Huntington, N.D.; Tarlinton, D.M. CD19 Differentially Regulates BCR Signalling through the Recruitment of PI3K. Autoimmunity 2014, 47, 430–437. [Google Scholar] [CrossRef] [PubMed]
41. Adamo, L.; Rocha-Resende, C.; Lin, C.-Y.; Evans, S.; Williams, J.; Dun, H.; Li, W.; Mpoy, C.; Andhey, P.S.; Rogers, B.E.; et al. Myocardial B Cells Are a Subset of Circulating Lymphocytes with Delayed Transit through the Heart. JCI Insight 2020, 5, e134700. [Google Scholar] [CrossRef]
42. Gorski, P.A.; Ceholski, D.K.; Hajjar, R.J. Altered Myocardial Calcium Cycling and Energetics in Heart Failure—A Rational Approach for Disease Treatment. Cell Metab. 2015, 21, 183–194. [Google Scholar] [CrossRef] [PubMed]