Investigation of hub genes and immune status in heart transplant rejection using endomyocardial biopsies

Abstract T cell‒mediated rejection (TCMR) and antibody‐mediated rejection (ABMR) are severe post‐transplantation complications for heart transplantation (HTx), whose molecular and immunological pathogenesis remains unclear. In the present study, the mRNA microarray data set GSE124897 containing 645 stable, 52 TCMR and 144 ABMR endomyocardial biopsies was obtained to screen for differentially expressed genes (DEGs) between rejected and stable HTx samples and to investigate immune cell infiltration. Functional enrichment analyses indicated roles of the DEGs primarily in immune‐related mechanisms. Protein‐protein interaction networks were then constructed, and ICAM1, CD44, HLA‐A and HLA‐B were identified as hub genes using the maximal clique centrality method. Immune cell infiltration analysis revealed differences in adaptive and innate immune cell populations between TCMR, ABMR and stable HTx samples. Additionally, hub gene expression levels significantly correlated with the degree and composition of immune cell infiltration in HTx rejection samples. Furthermore, drug‐gene interactions were constructed, and 12 FDA‐approved drugs were predicted to target hub genes. Finally, an external GSE2596 data set was used to validate the expression of the hub genes, and ROC curves indicated all four hub genes had promising diagnostic value for HTx rejection. This study provides a comprehensive perspective of molecular and immunological regulatory mechanisms underlying HTx rejection.

diagnostic system, the Molecular Microscope™ Diagnostic System (MMDx), assesses EMBs in terms of their molecular phenotype, including stable (normal), TCMR, ABMR and injury; the microarray data are available in GSE124897. 10,11 In the present study, we use bioinformatics analysis to screen the differences in gene expression between TCMR/ ABMR and stable HTx samples and identify hub genes. Additionally, abnormal immune cell infiltration in TCMR/ABMR samples and the significant relationship between hub gene expression and immune cell populations, further implicates these hub genes in HTx immunity.
Finally, the diagnostic value of hub gene expression in detecting HTx rejection was validated using the GSE2596 data set.
High-throughput microarray technologies have been increasingly used for the comprehensive identification of potential therapeutic target genes or biomarkers in several kinds of heart diseases. [12][13][14] However, few studies have investigated the molecular and immunological regulatory mechanisms of HTx rejection by performing bioinformatics analysis. Therefore, the aim of the present study was to identify candidate hub genes and investigate the immune status in HTx rejection through bioinformatics analysis, which deepened our understanding of the diagnosis and treatment of HTx rejection.

| Data collection and identification of differentially expressed genes (DEGs)
The mRNA microarray data set GSE124897 deposited by Chang et al was downloaded from the GEO database (http://www.ncbi.nlm. nih.gov/geo/), which contained EMBs from HTx recipients, including 645 stable, 52 TCMR and 144 ABMR samples (NCT02670408, an INTERHEART study). The inclusion criteria for selecting the microarray data set were set as follows: (a) the samples detected are heart transplant tissues from homo sapiens, (b) tissues are diagnosed with rejected and stable heart transplant tissues, (c) gene expression profiling of mRNA, and (d) the sample size is greater than 200. The available information on clinical patient/biopsy characteristics is obtained from the results of NCT02670408 (Table 1). The GPL15207 [PrimeView] Affymetrix Human Gene Expression Array was utilized to obtain gene expression profiles. DEGs between TCMR, ABMR and stable HTx samples were screened using the R 3.5.0 software and the LIMMA package. Adjusted P values <.05 and |log 2 fold change (FC)|> 1 were set as cut-off standards and indicate statistical significance. 15

| Functional enrichment analyses of DEGs
Functional enrichment analyses of DEGs were performed using Metascape, a powerful tool that can annotate large numbers of genes and perform enrichment analysis. 16 This tool integrates several authoritative functional databases such as Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome.

| Protein-protein interaction (PPI) network construction and analysis
Search Tool for the Retrieval of Interacting Genes/Proteins (STRING 11.0; http://string.embl.de/) was used to trace and predict the proteinprotein interaction (PPI) network after the DEGs were imported into the database, with a combined score >0. 4. 17 An open-source visualization software Cytoscape 3.6.1 was used to construct and visualize the PPI network. Subsequently, the Molecular Complex Detection (MCODE) plug-in was used to screen the most significant modules in the PPI network, with node score cut-off = 0.2, K-Core = 2, max depth = 100, and degree cut-off = 4 set as the selection criteria. Furthermore, the cyto-Hubba plug-in was utilized with the maximal clique centrality (MCC) method to identify hub genes in the PPI network.

| Immune cell infiltration evaluation and analysis
CIBERSORT is an analytical tool used to estimate the relative proportion of 22 human immune cells using gene expression data. 18 We uploaded the normalized gene expression profiles to the CIBERSORT web portal and set the algorithm to 100 permutations. Samples that met the CIBERSORT P < .05 requirements were considered eligible for continued analysis. In each sample, the combined proportions of all immunocyte types were set to equal 1.

| Prediction of drug-gene interactions
Drugs approved by the Food and Drug Administration (FDA) were selected through the Drug-Gene Interaction Database (DGIdb; http://www.dgidb.org/search_inter actions) based on the hub genes identified as promising targets. Drug-gene interactions were constructed and visualized using Cytoscape.

| Validation of expression levels and diagnostic value of hub genes
The GSE2596 data set, which includes 27 stable and 16 rejection HTx recipients, was used as a verification data set to evaluate the expression levels of the identified hub genes. Furthermore, receiver operating characteristic (ROC) analysis was used to assess the performance of hub genes in this research using GraphPad Prism 7.0 software. The area under the ROC curve (AUC) was then computed.

| Statistical analysis
GraphPad Prism 7.0 software was used for statistical analysis: (a) Two-tailed Student's t test was used to analyse the differences between immune cell populations in eligible rejection and stable HTx samples. (b) ROC curve analysis was used to determine the diagnostic effectiveness of hub genes in the verification GSE2596 data set.

| Identification of differentially expressed genes (DEGs) in heart transplantation (HTx) rejection endomyocardial biopsies (EMBs)
The normalization process of the data performed by the R LIMMA package is shown in Figure 1A. A total of 18,835 genes were detected in HTx EMBs. Compared with the stable group, 740 and 231 DEGs were identified in the TCMR and ABMR groups, respectively ( Table S1). The most differentially expressed gene in each of the TCMR and ABMR groups was C-X-C motif chemokine ligand 9 (log 2 FC = 4.04) and the C-X-C motif chemokine ligand 11 (log 2 FC = 3.80), respectively. The volcano plots for DEGs, TCMR vs stable and ABMR vs stable, are shown in Figure 1A,B, respectively, and the expression levels of the top 50 DEGs in the TCMR and ABMR groups are represented by heat maps in Figure 1C,D, respectively.

| Functional enrichment analyses of DEGs associated with HTx rejection reveal roles in the immune response
Functional enrichment analyses using the Metascape tool revealed that DEGs in both TCMR and ABMR groups were enriched mainly in pathways related to the immune response ( Figure 2A

| Protein-protein interaction (PPI) network analysis identified hub genes common to both TCMR and ABMR HTx
The PPI networks of DEGs in the TCMR and ABMR groups were constructed using the STRING database and visualized using Cytoscape software. PPI networks were constructed with 679 nodes and 10 511 edges in the TCMR group ( Figure 3A) and 184 nodes and 2146 edges in the ABMR group ( Figure 3B). The three top-ranked modules from each PPI network, TCMR and ABMR groups, respectively, were identified using the MCODE plug-in (Figure 3). Functional enrichment analyses revealed that DEGs in these modules were mainly associated with immune-related processes, suggesting that immune dysregulation plays a critical role in HTx rejection (Table 2).
To identify hub genes, the cytoHubba plug-in was applied to identify the top 10 DEGs in the PPI network ranked by the MCC

| Immune cell infiltration is altered in TCMR and ABMR samples
According to the CIBERSORT algorithm, 143 stable, 52 TCMR and 142 ABMR HTx EMBs that matched the requirements of CIBERSORT P < .05 were identified after filtering of data (Table S2)

| Significant correlations between hub gene levels and immune cell proportions
Correlation analysis revealed relationships between hub gene levels and immune cell proportions in TCMR and ABMR groups ( Figure 6).
In the TCMR group, ICAM1 expression significantly correlated with

| Validation of hub genes
The external GSE2596 data set was used to validate the expression levels and diagnostic value of the hub genes.

| D ISCUSS I ON
T cell-mediated rejection and ABMR have emerged as major risk factors for limiting heart allograft survival. [19][20][21] However, the differences in gene expression and immune status of rejection and stable heart graft samples are still unknown. Bioinformatics analysis identified a total of 740 and 231 DEGs in the TCMR and ABMR groups, respectively. Functional enrichment analysis revealed that DEGs in the TCMR and ABMR groups belong mainly to processes associated with the immune response, indicative of the abnormal immune regulation that occurs during cardiac allograft rejection. 22 In the PPI network construction, the top-ranked DEGs were se- which contributes to cardiac fibrosis and dysfunction and heart failure. 38 CD44 is also crucial in the development of cardiac remodelling and myocardial fibrosis, which accelerates the progression to heart failure. 39 In addition, HLA-A and HLA-B also play potential roles in heart failure. 40 These findings suggest that it is worth investigating the exact role of these hub genes in heart failure following HTx and whether the expression of these hub genes in different myocardial regions has different patterns to affect the development of heart failure following HTx.
The CIBERSORT algorithm defined the immune cell infiltration characteristics in the HTx samples and revealed abnormal adaptive and innate immunity in the TCMR and ABMR groups. Significant differences in the proportions of several types of immune cells were found in our study, which were consistent with other studies: (a) The abundance of B cells, including CD20 + CD27 + memory B cells and CD20 -CD138 + plasma cells, in heart graft infiltrates is closely associated with cardiac allograft vasculopathy. 41

(b)
Activated CD4 + /CD8 + T cell responses play a critical role in heart graft rejection. 19,42,43 (c) Activated pro-inflammatory follicular helper T cells help to generate donor-specific antibody (DSA) responses after organ transplantation, which is an attractive target for improving the effects of immunosuppressive drugs. 44,45 (d) Treg therapy is a novel tool for delaying graft rejection following solid organ transplantation. 46 Expanding the frequency of anti-inflammatory FoxP3 + Tregs prevents heart graft rejection and prolongs graft survival in mice. [47][48][49] (e) Activated γδ T cells promote heart graft rejection in mice by inducing IL-17 production, which can be reversed by depletion of γδ T cells. 50 macrophages are essential in heart allograft rejection, depletion of which can promote graft acceptance in mice. 56,57 Analysis in this study revealed that the expression levels of hub   Inevitably, the present study had several limitations. Some key information on clinical and laboratory characteristics of the patients in the present study such as comorbidities, medications, cardiac function and haemodynamics, and laboratory biomarkers were unavailable from the results of NCT02670408. What's more, some patient/biopsy characteristics may potentially affect on our results, which needs further exploration. In addition, the present study lacked further experimental verification as a solid foundation. However, the results of the current study provide a strong basis on which to develop future studies to obtain a more detailed understanding of the molecular and immune mechanisms underlying HTx rejection, and thus improving strategies to diagnose and prevent HTx rejection.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Formal analysis (equal); Funding acquisition (equal); Project administration (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
The datasests generated and/or analyzed during the current study are available from the Gene Expression Omnibus GSE124897 and