The gene expression profile dataset GSE6381 was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). GSE6381 presented early gene expression profiles during intraoperative myocardial I/R in cardiac surgery under cardiopulmonary bypass. Specifically, human right ventricular samples were obtained during surgical repair of ventricular septal defect and resection of right ventricular outflow obstruction, with four replicates at each stage (immediately post-cardioplegic arrest, end-ischemia, and 5 min after reperfusion). Data were detected on the Affymetrix Human Genome U133A Array platform (GPL96). The limma package in R language was used to identify DEGs for end-ischemia vs. cardiac arrest, reperfusion vs. end-ischemia, and reperfusion vs. cardiac arrest, respectively (Ritchie et al., 2015). The thresholds for DEGs selection were set as p < 0.05 and log₂ FC (fold change) > 1. Volcano plots and heatmaps were generated to overview the gene expression profiles by using the R package.

For the DEGs from each comparison, pathway and process enrichment analyses were performed by using the Metascape (http://metascape.org/) (Zhou et al., 2019). The analyses were based on the following databases: GO Biological Processes, KEGG Pathways, Reactome Gene Sets, and Canonical Pathways. Terms with count ≥3, p-value < 0.01, and enrichment factor (observed counts/counts expected by chance) >1.5 were grouped in clusters based on the similarities. Specifically, p-values were calculated based on the accumulative hypergeometric distribution, and q-values were calculated by using the Banjamini-Hochberg procedure for multiple testing (Hochberg and Benjamini, 1990). Hierarchical clustering was performed on the enriched terms with the kappa scores, and sub-trees with a similarity >0.3 were considered a cluster (Roberts, 2008). The most statistically significant term within a cluster was used to represent the cluster. The overlaps between three gene lists were shown in a Circos plot (Krzywinski et al., 2009). Heatmaps of the top 20 and 100 enriched terms were generated by using the JTreeView (Saldanha, 2004).

To further investigate the relationships among the enriched terms, network plots of the top 20 clusters were generated, where terms with a similarity >0.3 were connected by edges. The networks were visualized by using the Cytoscape software v3.7.2 (https://cytoscape.org/), where each node represented an enriched term and was colored by its cluster ID and then by p-value (Shannon et al., 2003). In addition, the nodes were represented as pie charts, where the size of a pie was proportional to the total number of hits belonging to the specific term. The pie charts were color-coded based on the gene list identities, and the size of a slice showed the percentage of genes under the term which originated from the specific gene list.

For all DEGs, PPI network analysis was carried out by using the Metascape. The analyses were based on several databases, including BioGrid, InWeb_IM, and OmniPath (Stark et al., 2006; Turei et al., 2016; Li et al., 2017). The resultant network contained the subset of proteins with at least one physical interaction with one other member. The Molecular Complex Detection (MCODE) algorithm was applied to identify highly connected network components (Bader and Hogue, 2003). The PPI network plots were visualized by using the Cytoscape, where each node represented a hub protein and was colored by the cluster ID and gene list identity.

For the 26 key genes identified in MOCDE 1, GO enrichment analyses of biological process (BP), molecular function (MF), and cellular component (CC) were conducted by using the gProfiler (https://biit.cs.ut.ee/gprofiler/) (Raudvere et al., 2019). In addition, KEGG pathway analysis was performed to analyze significantly enriched signaling pathways. A p-value < 0.05 was applied as the threshold value.

By using the GeneMANIA (http://genemania.org/), a flexible and user-friendly web interface, PPI network and functional assays of 26 key genes were generated and visualized (Zuberi et al., 2013). Several types of relationships among proteins included co-expressions, physical interactions, shared protein domains, pathways, website predictions, co-localizations, and genetic interactions. GeneMANIA also reported weights that indicated the predictive value of each selected data set.

GSEA is a computational method used to determine whether a predefined set of genes shows significant and concordant differences between the two study groups. GSEA software v3.0 (http://www.gsea-msigdb.org/gsea/) was used to perform KEGG pathway enrichment for two comparisons: end-ischemia vs. cardiac arrest and reperfusion vs. cardiac arrest (Subramanian et al., 2005). A p-value < 0.05 was considered as statistically significant. The “c2.cp.kegg.v6.2.symbols.gmt” gene set was downloaded from the Molecular Signatures Database and used as the reference gene set.

The protocol of animal study was approved by the Ethic Committee for Animal Care and Use at Soochow University, Suzhou, Jiangsu, China. The animal experimental procedures followed the Guide for the Care and Use of Laboratory Animals (National Institutes of Health publication No. 85-23, revised in 1996). Adult male C57BL/6J mice (8 weeks old, 20–25 g) were provided by the Experimental Animal Center of Soochow University. Mice were housed in a controlled condition (temperature of 24–26°C, 12-h light/dark cycle, and free access to food and water).

We induced myocardial I/R injury in mice as previously reported (Peng et al., 2021). First, the mice were anesthetized with intraperitoneal sodium pentobarbital (50 mg/kg). Next, the mice were endotracheally intubated and received mechanical ventilation using a rodent respirator (RWD, Shenzhen, China). A left thoracotomy was performed in the left fourth intercostal space. The left anterior descending coronary artery (LAD) was temporarily ligated at ∼2 mm under the lower edge of the left auricle. Observation of tissue blanching in the left ventricular myocardium and elevation in the ST-segment of electrocardiography indicated that the myocardial I/R model was successfully established. The LAD was occluded for 30 min to induce myocardial ischemia, and the ligature was removed for reperfusion of 30 min. For mice in the sham surgery group, the same surgical procedure was performed except that the suture placed around the LAD was not tied. At the end of the reperfusion, the hearts were harvested.

Evans blue staining and 2,3,5-triphenyltetrazolium chloride (TTC) staining were used to assess the area at risk and infarct size of myocardium, as previously described (Peng et al., 2020; Song et al., 2020; Yang et al., 2021). Briefly, the LAD was re-ligated following 30 min of reperfusion, and 1 ml of 2% Evans blue dye (Sigma, St. Louis, MO) was perfused into the coronary artery. The heart was then immediately excised and quickly frozen at −80°C. The heart samples below the ligature were sliced into five transverse sections, and the slices were incubated with 1% TTC (Sigma, St. Louis, MO) at 37°C for 15 min in a dark chamber and fixed with 10% formaldehyde overnight. The non-ischemic area (blue), the myocardium area at risk (AAR, red), and infarct area (IA, pale) were measured using the Image Pro Plus software (Media Cybernetics, Silver Spring, MD, United States). The ratios of IA/AAR and AAR/left ventricular area (LV) were calculated.

Histological structure of the left ventricular myocardium samples was detected using the hematoxylin and eosin (HE) staining (Yang et al., 2021). Briefly, myocardial tissues were fixed with paraformaldehyde and embedded with paraffin. The samples were cut transversely into 5 µm-thick sections. After dewaxing with dimethylbenzene and dehydrating with graded ethanol series, the slices were stained with hematoxylin and eosin. Finally, histological damage of the myocardium was observed under a light microscope (Nikon Corporation; Tokyo, Japan).

Cell apoptosis in the left ventricular myocardium samples was assessed using the terminal deoxynucleotidyl transferase dUTP nick‐end labeling (TUNEL) analysis with the in situ Cell Death Detection kit (Roche Molecular Biochemicals, Mannheim, Germany) (Peng et al., 2020; Peng et al., 2021). The heart samples were fixed in 4% paraformaldehyde, dehydrated, embedded in paraffin, and cut into 4 μm‐thick sections. The slices were incubated with TUNEL reaction mixture and then counterstained with the 4′,6-diamino-2-phenylindole (DAPI) to show the nuclei. At least four random and nonoverlapping of views in areas surrounding the infarct area in the I/R group or the corresponding area in the sham group were observed using a fluorescence microscope (Nikon Corporation; Tokyo, Japan). The apoptosis rate was expressed as a percentage of TUNEL-positive cells to total cell nuclei.

Total RNA in the left ventricular myocardium samples was extracted using the TRIZOL reagent (Thermo Fisher Scientific, Waltham, MA). After RNA quantification and purity evaluation, reverse transcription was conducted with a cDNA Synthesis Kit (abm, Richmond, BC, Canada), in a 20 μl reaction system containing 5× All-in-one, DEPC water, and RNA samples. Then, qPCR was performed with the EvaGreen (abm, Richmond, BC, Canada) on the Roche Light Cycler R480 System (Roche, Bedford, MA). The amplification conditions were: Pre-denaturation at 93°C for 10 s, denaturation at 55°C for 15 s, and annealing at 72°C for 20 s for 40 cycles. Three replicates were used for each sample. Gene expression were analyzed by using the 2^−ΔΔCT method and normalized to β-tubulin. The sequences of primers (Sangon Biotech, Shanghai, China) are listed in Table 1.

The total protein of the left ventricular myocardium samples was extracted using the RIPA reagents (Beyotime, Shanghai, China). The protein concentrations were detected using a bicinchoninic acid protein assay kit (Beyotime, Shanghai, China). All proteins were separated with 10% sodium dodecyl sulphate-polyacrylamide gel electrophoresis, and transferred onto polyvinylidene fluoride membranes (Millipore Corp., Bedford, MA). The membranes were incubated with the primary antibodies overnight at 4°C, and then incubated with the horseradish peroxidase-conjugated secondary antibodies for 2 h at room temperature. The primary antibodies were as follows: goat anti-CCL21 (1:1000; R&D Systems, Minneapolis, MN), rabbit anti-XCR1 (1:1000; Novus Biologicals, Littleton, CO), goat anti-CXCL13 (1:1000; R&D Systems, Minneapolis, MN), goat anti-CCR7 (1:1000; Novus Biologicals, Littleton, CO), rabbit anti-GPR174 (1:1000; Biorbyt, Cambridge, United Kingdom), rabbit anti-CXCR5 (1:1000; Absin, Shanghai, China), and rabbit anti-β-tubulin (1:5000; Cell Signaling Technology, Beverly, MA). Finally, the protein bands were detected using the Tanon 5200 luminescent workstation (Tanon, Shanghai, China). The protein expression was normalized to β-tubulin as control.

Following 30 min of reperfusion, the left ventricular myocardium samples were fixed with 10% formalin and embedded in paraffin. Then, the tissues were cut transversely into 5 µm-thick sections, deparaffinized, rehydrated, and subjected to antigen retrieval. For immunofluorescence staining, the slices were incubated with primary antibodies overnight at 4°C, and then incubated with secondary antibodies for 2 h at room temperature. The following primary antibodies were used: goat anti-CCR7 (1:200; Novus Biologicals, Littleton, CO), rabbit anti-GPR174 (1:200; Biorbyt, Cambridge, United Kingdom), mouse anti-α-Actinin (1:250; Santa Cruz Technology, CA), and mouse anti-Vimentin (1:200; Abcam, CA). The nuclei were counterstained with DAPI. At least four random and nonoverlapping of views in areas surrounding the infarct area in the I/R group or the corresponding area in the sham group were observed using the fluorescence microscope. Images were analyzed using the Image Pro Plus software.

Statistical analysis was performed by using the GraphPad Prism software (version 9.0, GraphPad, San Diego, CA). Data were expressed as mean ± standard error of the mean and compared using t-test. A two-tailed value of p < 0.05 was considered statistically significant.

The expression profiles of 12,549 genes in GSE6381 of human right ventricular samples during myocardial I/R process were analyzed. The DEGs were identified in three comparisons (Figure 1A): 309 DEGs (247 upregulated and 62 downregulated) in end-ischemia vs. cardiac arrest (Supplementary Table S1), 175 DEGs (90 upregulated and 85 downregulated) in reperfusion vs. end-ischemia (Supplementary Table S2), and 217 DEGs (165 up-regulated and 52 down-regulated) in reperfusion vs. cardiac arrest (Supplementary Table S3). These DEGs showed apparently differential expression profiles (Figure 1B). The top 10 up- and downregulated DEGs for each comparison are presented in Table 2. There were 40, 34, and 51 common DEGs between two comparisons, but no gene was found as common to the three (Figure 1C; Supplementary Table S4).

The pathway and process enrichment analyses of DEGs for each comparison were conducted by using the Metascape. Figure 2A shows gene overlaps and genes that belong to the same enriched ontology term. The top 20 enriched terms are colored by p-values in the heatmap (Figure 2B), with the characteristics shown in Figure 2C. Of these, inorganic ion homeostasis, positive regulation of JNK cascade, and cytokine-mediated signaling pathway may be closely related to the myocardial I/R process. Neuroactive ligand-receptor interaction is the common KEGG pathway. The top 100 enriched terms are available in Supplementary Figure S1. Furthermore, the relationships of the significantly enriched terms are visualized in networks, colored by cluster IDs (Figure 2D), p-values (Figure 2E), and gene list identities (Figure 2F). All results of the enrichment analyses are listed in Supplementary Table S5.

The PPI networks of DEGs were performed to identify key modules with highly clustered proteins using the MCODE functionality. All interactions of the proteins are presented in Supplementary Figure S2. The PPI analyses revealed a total of 10 key modules. The MCODE one contained 26 highly clustered proteins: C5AR1, CASR, CCL21, CHRM4, CNR1, CORT, CXCL13, EDN1, EDN3, GNRH2, GNRHR, GRM8, HCRTR1, HCRTR2, MTNR1A, NMBR, NPY1R, NTS, NTSR1, OPRD1, OPRK1, P2RY10, P2RY6, PPY, TACR3, and XCR1. In the MCODE networks, the interactions among the 26 hub proteins were colored by MCODE cluster IDs (Figure 3A) and identities of the proteins (Figure 3B). The annotations of the 26 genes are presented in Table 3.

Based on the 26 genes, GO and KEGG enrichment analyses were performed using the gProfiler. Of the top significantly enriched BP, CC, and MF terms, cellular calcium ion homeostasis (GO:0006874) and calcium ion homeostasis (GO:0055074) were noted in the process of myocardial I/R injury (Figure 3C). Thirteen genes were involved, including P2RY10, HCRTR1, HCRTR2, NTSR1, CASR, XCR1, P2RY6, EDN3, EDN1, CNR1, CCL21, C5AR1, and CXCL13. In addition, KEGG analysis showed one enriched pathway, neuroactive ligand-receptor interaction (Figure 3C), in line with the previous results (Figure 2C). All enrichment results are available in Supplementary Figure S3; Supplementary Table S6.

To further reveal the interactions of the 26 hub proteins, PPI network and functional analysis were conducted using the GeneMANIA (Figure 3D). There were substantial and complex relationships within these proteins, including co-expressions, physical interactions, shared protein domains, pathways, predicted interactions, co-localizations, and genetic interactions (Supplementary Table S7). In addition, the functional analysis revealed several gene functions (Supplementary Table S8). Of note, there were three functions of calcium regulation (cellular calcium ion homeostasis, calcium ion homeostasis, and positive regulation of cytosolic calcium ion concentration). In these calcium-related functions, five hub proteins (CCL21, XCR1, CXCL13, EDN1, and CASR) were revealed, which were significantly enriched in two calcium processes (Figure 3C).

Based on the ischemia samples, GSEA revealed several significantly enriched pathways, including neuroactive ligand receptor interaction, GNRH signaling pathway, calcium signaling pathway, WNT signaling pathway, VEGF signaling pathway, and JAK-STAT signaling pathway (Figure 4A). Based on the reperfusion samples, GSEA revealed several significantly enriched pathways, including neuroactive ligand receptor interaction, VEGF signaling pathway, JAK-STAT signaling pathway, and WNT signaling pathway (Figure 4B). All results of GSEA enriched pathways are available in Supplementary Tables S9, S10.

Compared to the sham controls, myocardial I/R injury led to an obvious myocardial infract size (approximately 40%) in the Evans blue/TTC staining (Figure 5A). In addition, HE staining showed remarkable myocardial histological damage following myocardial I/R injury (Figure 5B). The structure of myocardial tissue in the sham group was clear and aligned, with intact myocardial fibers; however, the I/R group showed a damaged tissue structure, swollen cardiomyocytes, disarranged myofibrils, and leukocytes infiltration. The TUNEL assay was performed to assess the level of DNA fragmentation in myocardial apoptotic cells. The results showed a significantly increased TUNEL apoptosis rate induced by myocardial I/R injury (Figure 5C).

Based on the PPI network and hub protein analyses, five key genes (CCL21, XCR1, CXCL13, EDN1, and CASR) were identified. To investigate the role of these genes at the early stage of myocardial I/R injury, the mRNA expression of these genes was determined using the qPCR. The results showed that the mRNA expression of CCL21, XCR1, CXCL13, and CASR were significantly up-regulated in the left ventricular samples of myocardial I/R injury (Figure 6A).

Given that chemokines are small secreted proteins with chemoattractant properties that play a key role in inflammation and cell apoptosis (Jarrah et al., 2018; Zhang et al., 2020), we assessed the protein expression of CCL21, XCR1, and CXCL13. The results showed that the protein expression of CCL21 and CXCL13 was significantly increased following myocardial I/R injury (Figure 6B). Next, we tested the protein expression of GPR174 and CCR7 (known receptors of CCL21) (Chen et al., 2020; Van Raemdonck et al., 2020; Zhao et al., 2020) and CXCR5 (a known receptor of CXCL13) (Hussain et al., 2019; Shen et al., 2020). We found that the protein expression of GPR174 and CCR7, but not CXCR5, was significantly elevated. These results suggested that the CCL21-GPR174/CCR7 chemokine axis was involved in the early stage of myocardial I/R injury (Figure 6C).

We detected the expression and cellular localization of GPR174 and CCR7 using the immunofluorescence staining. The double-labeling immunofluorescence results showed that there was no GPR174 positive cells (green) co-localized with the α-Actinin positive cells (red) in the sham or myocardial I/R mice, suggesting that GPR174 did not present on cardiomyocytes (Figure 7A). Therefore, we further investigated whether GPR174 was expressed in the cardiac fibroblasts. Following myocardial I/R injury, GPR174 positive cells (green) and Vimentin positive cells (red) co-localized in the left ventricular myocardial tissues, and the expression of GPR174 on cardiac fibroblasts was significantly increased (Figure 7B). Similarly, we also found that myocardial I/R-induced activation of CCR7 was also evident on cardiac fibroblasts (Figure 7C). Taken together, these results suggested that the up-regulated GPR174 and CCR7 on the cardiac fibroblasts underlie the pathological course of myocardial I/R injury.