Computational Prediction of Biomarkers, Pathways, and New Target Drugs in the Pathogenesis of Immune-Based Diseases Regarding Kidney Transplantation Rejection

Background The diagnosis of graft rejection in kidney transplantation (KT) patients is made by evaluating the histological characteristics of biopsy samples. The evolution of omics sciences and bioinformatics techniques has contributed to the advancement in searching and predicting biomarkers, pathways, and new target drugs that allow a more precise and less invasive diagnosis. The aim was to search for differentially expressed genes (DEGs) in patients with/without antibody-mediated rejection (AMR) and find essential cells involved in AMR, new target drugs, protein-protein interactions (PPI), and know their functional and biological analysis. Material and Methods Four GEO databases of kidney biopsies of kidney transplantation with/without AMR were analyzed. The infiltrating leukocyte populations in the graft, new target drugs, protein-protein interactions (PPI), functional and biological analysis were studied by different bioinformatics tools. Results Our results show DEGs and the infiltrating leukocyte populations in the graft. There is an increase in the expression of genes related to different stages of the activation of the immune system, antigenic presentation such as antibody-mediated cytotoxicity, or leukocyte migration during AMR. The importance of the IRF/STAT1 pathways of response to IFN in controlling the expression of genes related to humoral rejection. The genes of this biological pathway were postulated as potential therapeutic targets and biomarkers of AMR. These biological processes correlated showed the infiltration of NK cells and monocytes towards the allograft. Besides the increase in dendritic cell maturation, it plays a central role in mediating the damage suffered by the graft during AMR. Computational approaches to the search for new therapeutic uses of approved target drugs also showed that imatinib might theoretically be helpful in KT for the prevention and/or treatment of AMR. Conclusion Our results suggest the importance of the IRF/STAT1 pathways in humoral kidney rejection. NK cells and monocytes in graft damage have an essential role during rejection, and imatinib improves KT outcomes. Our results will have to be validated for the potential use of overexpressed genes as rejection biomarkers that can be used as diagnostic and prognostic markers and as therapeutic targets to avoid graft rejection in patients undergoing kidney transplantation.


INTRODUCTION
The severity and occurrence of rejection in kidney transplant (KT) patients depend on numerous variables that can affect the magnitude and nature of immune responses. Understanding how genetic and molecular factors affect the effector functions of immune cells and donor-specific antibodies (DSA) can better renal stratification receptors based on their immunological risk and thus help the clinician make better decisions to anticipate adverse events (1)(2)(3)(4)(5).
In the last two decades, high-throughput technologies such as next-generation sequencing (NGS) and microarrays have been developed. In parallel, international public data repositories have been developed, such as the GEO (Gene Expression Omnibus) database of the National Center for Biotechnology Information (NCBI) (6), the Array Express database of the Institute European Bioinformatics (EBI) (7), in order to store and distribute these data. The large amount of data generated by these technologies makes it necessary to use robust bioinformatics tools that allow us to know in silico how molecules interact and regulate the different biological processes that mediate the biological processes of health and disease.
Therefore, functional analysis tools have been developed to explore and identify critical biological processes. These tools are based on biological knowledge databases such as Gene Ontology (GO) or the Kyoto Encyclopedia of Genes and Genomes (KEGG). GO provides a controlled vocabulary of terms (ontologies) to describe gene products in terms of biological processes, cellular components, and associated molecular functions (8). On the other hand, the KEGG is a reference database for the biological interpretation of metabolism and cellular processes (9). Furthermore, the development of genomic and protein databases has led to many bioinformatics tools to predict the properties of proteins and the genome: splice sites, posttranslational modifications, stability, pathogenicity. The integration of these new tools and data offers unprecedented opportunities that promise to accelerate and improve our understanding of biology and the search for biomarkers, and specifically for the field of organ transplantation.
The aim was to search for differentially expressed genes (DEGs) in AMR patients and to analyze the types of populations of leukocytes infiltrated in the graft developing diagnostic models using computational predictions in order to find diagnostic and prognostic markers and as a therapeutic target that allows avoiding graft rejection in patients undergoing kidney transplantation.

Gene Expression Data of Gene Expression Omnibus Repository
For this study, the raw gene expression data studies were downloaded from the Gene Expression Omnibus (GEO) database (6). The characteristics of the GSE cohorts are summarized in Table 1. Immunosuppressive treatment in the different GEO studies was similar. Around 40-60% of the transplant recipients in the different studies had triple therapy of corticosteroids + tacrolimus + MMF, while the rest had other combinations in different proportions: corticosteroids + MMF + cyclosporine, tacro + corticosteroids. Regarding induction therapy (thymoglobulin and baxilysimab) there is no information available in the studies. All the bases analyzed in this study have been reviewed and are comparable.

Venn Diagrams
Venn diagrams (10) were used to determine the intersection and analysis of the four analyzed cohorts' differentially expressed genes (DEGs).

Functional Analyzes of the KEGG and GO Terms
DEGs in this study were entered into the WebGestalt web tool (6) for the analyzes of Gene Ontology terms (GO) and biological pathways from the English Kyoto Encyclopedia of Genes and Genomes (KEGG) (11). The minimum number of genes to be considered a pathway was adjusted to a value of two.
The Benjamini and Hockberg FDR method adjusted the pvalue obtained from each biological pathway with FDR<0.05 considered significant. Only those from the molecular function and biological processes were taken into account for the analyses of GO terms. GO terms with an adjusted p value<0.05 were considered significant (12). Furthermore, the enrichment ratio (ER) is indicated for each GO term and KEGG pathway in this analysis. The ER is defined as the number of genes observed among the number of genes expected for each GO term or KEGG pathway (12).

Distribution of Leukocytes in Renal Graft Biopsy
The xCell algorithm (13) was used to estimate the relative distribution of leukocytes in graft biopsy, which allows the inference of cell subpopulations from transcriptomic data. This algorithm allows the relative proportion of 64 cell types to be inferred, including lymphoid, myeloid, hematopoietic progenitors, stromal cells, and other cell types, such as epithelial cells or melanocytes. In our study, we focus exclusively on lymphoid and myeloid cells.
The gene expression data of GSE36059, GSE44131, GSE50084, GSE93658 was analyzed using the xCell tool (Figure 1), generating a score for each cell population proportional to the relative frequency of that population. The score obtained from each subpopulation and sample will be used later for comparison by traditional statistical methods between transplants with and without AMR.

Protein-Protein Interaction Networks
The interaction between the proteins produced by the DEGs obtained in the previous sections was used to generate a proteinprotein interaction map (PPI) (14) using the Network Analyst web tool (15). This tool uses the Innate DB database of protein interactions.
The nodes represent the proteins in PPI networks, while the junction lines represent the known interactions between the bound proteins. The constructed network was limited to include only the core proteins (zero-order interactions) to avoid the hairball effect, and better visualize the network. The network's topological properties were evaluated using the betweenness centrality (BC) and the Degree Centrality (GCy). The GI is the number of connections a node has with other nodes, while BC measures the shortest paths that traverse a node. It indicates the degree to which a node is located between other nodes in the network. Nodes with a high BC and DCy index are essential proteins in cell signaling and control information flow.
The modules of the PPI network are tightly connected areas of the network that are considered relatively independent as they often work together to perform a biological function. Therefore, the modules of the PPI networks usually correspond to metabolic pathways or protein complexes. The module explorer of the Network Analyst tool, based on the random-walk algorithm, was used to identify the modules.
Connections within a module and the rest of the PPI network were analyzed. Values of p<0.05 were considered significant. The significance of each module was calculated using the Wilcoxon test. The significant modules were subsequently functionally analyzed using the KEGG and GO term databases. Those GO terms and KEGG pathways with a value of p<0.05 were considered significant.

Analyzes of Protein Kinases and Transcription Factors
In the search for biomarkers and new therapeutic targets in transplantation, it is crucial to study those genes with altered expression and the entire signaling pathway that regulates them. For this purpose, we employ the Expression2Kinases (X2K) tool that infers the protein kinases and transcription factors that form the regulatory networks upstream of the DEGs obtained previously. miRNA, B cells were inferred by uploading the list of these genes to the Expression2Kinase (X2K) web tool (15). X2K identifies the transcription factors that regulate DEGs using the ChEA database. Protein kinases associated with transcription factors were obtained through the Kinase Enrichment Analysis module of X2K. Transcription factors and protein kinases with p values <0.05 were considered significant.

Computational Prediction of New Target Drugs in Overexpressed Genes During AMR Kidney Transplantation
The computational prediction of new target drugs overexpressed genes during AMR kidney transplantation used the Gene2Drug web tool (16,17).
The Gene2drug algorithm searches for drugs or low molecular weight compounds that induce a significant variation in the expression levels of a list of genes of interest. Gene2drug uses gene expression data from the Connectivity Map database (CMap), containing various cell lines treated with 1309 drugs or low molecular weight molecules. CMap identifies drugassociated transcriptomics alterations. CMap show which genes change their expression induced by these compounds (18). The genes introduced in the Gene2drug algorithm had a value of BC≥100 in the topological analysis of the PPI network, presenting a high connection (Table S1). Gene2Drug calculates the English Enrichment Score (ES value) and the p-value corrected by the Benjamini-Hochberg method for the 1309 compounds and drugs of CMap. The ES value ranges from -1 to +1 and represents the degree to which a compound or drug regulates the genes. The ES value uses a generalization of the Kolgomorov-Smirnov method. Compounds with ES values close to +1 will be those that the algorithm predicts will significantly induce the genes of interest, while compounds with values at -1 will be those that will inhibit the genes of interest. Our objective was to search drugs that counteract to focus on those compounds with ES ≤0 values.

Protein-Protein (STRING) and Drug-Protein (STITCH) Interactions Databases
STRING is a database of protein-protein interactions (16), which includes direct (physical) and indirect (functional) associations obtained from experimental data and computational predictions. STRING version 11.0 was used to obtain the interactions between the PDGFRA and PDGFRB receptors and the proteins with the highest BC and GI values from the PPI network previously analyses (GZMB, STAT1, LYN, IRF1, and IRF8).
STITCH is also a database of computational drug-protein interactions (19). Like STRING, it includes direct (physical) and indirect (functional) associations obtained from experimental data and computational predictions (20). STITCH version 5.0 was used to obtain the drug interactions imatinib.

Statistical Analysis
The graphs and statistical analyzes were carried out in the Statistical Package for the Social Sciences (SPSS, version 22, Chicago, IL) and GraphPad Prism (version 6, San Diego, CA) software as well as in the R programming language, used for the latter the environment Integrated Development R Studio version 3.4.
The results have been expressed as the mean ± standard deviation (SD) for quantitative data or percentages for categorical data. For the comparison of categorical variables, the c 2 test or Fisher's exact test was used. The verification of the normality of the data was carried out using the Kolmogorov-Smirnov test. The Mann-Whitney U test was used to compare two groups with variables that did not adjust to normality. The Kruskal Wallis test and Dunn's post hoc test with Bonferroni correction for multiple comparisons compared three or more groups. Correlation analyzes were carried out using the Spearman index, as previously published (5,21,22). For the longitudinal comparison of two related groups, the Wilcoxon non-parametric test for related samples was used. The Friedman test with Wilcoxon post hoc was used to compare three or more related groups. The evaluation of the sensitivity and specificity of a biomarker was carried out by constructing ROC curves. The discriminatory capacity was evaluated using the AUC ( Figure 5). To obtain the optimal cut-off value that maximizes sensitivity and specificity, we use the Youden index. To correct the p-value in multiple comparisons, the Benjamini-Hochberg or Bonferroni method was used. Values of p<0.05, or p-corrected<0.05 in the case of multiple comparisons, were considered significant for all statistical tests.

Description of the Databases Used in the Different Bioinformatic Analyzes
The different databases used in the different bioinformatic analyzes designed to search for biomarkers in KT is shown in Figure 1.
Firstly, a search was conducted in the GEO database for studies performed with kidney graft biopsy samples from transplant recipients with AMR or without acute rejection (NR) ( Table 1). Next, the gene expression data in the GEO2R application to obtain the differentially expressed genes (DEGs) in the AMR samples were analyzed. The DEGs list was used in various bioinformatics analyses to obtain essential characteristics at the molecular and cellular level of the AMR, create diagnostic models, and search for new therapeutic options.

Differentially Expressed Genes in Biopsies With AMR
The gene expression data of the selected studies (GSE36059, GSE44131, GSE50084, and GSE93658) were analyzed using the GEO2R web tool establishing two groups for each cohort, one with biopsies from transplants NR and AMR.
The GEO2R tool obtained the Fold-change and adjusted p values using the Benjamini and Hochber FDR method, establishing DEGs between the NR and AMR groups with an FDR value <0.05. The number of DEGs obtained for each study is shown in Table 1. In the study, a total of 1778, 2901, 1387, 2855 DEGs were obtained in GSE36059, GSE44131, GSE50084, and GSE93658, respectively.
We selected those DEGs present in these four studies for further analyses and ruled out genes of little relevance to rejection. In all studies, the deregulation of these genes indicates that they play an essential role in the rejection process and are therefore potential biomarkers. Venn diagram analysis ( Figure 2) shows 340 overexpressed and 26 underexpressed DEGs in the AMR group. The complete list of DEGs obtained is shown in Table 2.

Biological and Functional Pathway Analyzes of Differentially Expressed Genes
The selected DEGs were subjected to functional analyzes of GO terms and biological pathways of the KEGG. The main biological pathways overrepresented in DEGs are shown in Table S2A, where two biological pathways related to AMR stand out, graftversus-host disease (GVHD; hsa05332, FDR<0.0001) and allograft rejection (hsa05330, FDR<0.0001).
The functional analyzes of GO terms for biological processes (BP) and molecular functions (MF) are summarized in Table 3.
For over-expressed DEGs, the most overrepresented BPs are those related to the immune response activation (GO:0002253, FDR<0.0001) and immune effector processes (GO:0002252, FDR<0.0001). For molecular functions, the most critical terms were   Table 3A) GO term with an FDR <0.05 was obtained for underrepresented DEGs due to the low number of genes included in the analysis (Table 3B). According to the unadjusted p values, the main BPs were lipid metabolism (GO:0006629, p=0.0001) and fatty acid oxidation (GO:0019395, p=0.0001). Regarding the molecular functions, they were the oxidoreductase activity of CH-OH groups (GO:0016614, p=0.0002) and the activity of GA-3P DHase (GO: 0043878, p=0.005).

Identification of Interaction Networks Between Proteins
The PPI network obtained was divided into three subnets, of which we analyzed in greater detail the largest one composed of 98 nodes and 154 interactions ( Figure 3). The degree of interconnection shows that 25 (25.5%) of the nodes were equal to one, while 73 (74.5%) show a degree ≥2. Of those 73, there are four nodes (STAT1, IRF1, LYN, and IRF8) that show a high GI (≥10): STAT1 (GI=21), IRF1 (GI=19), LYN (GI=10), and IRF8 (GI=10). Regarding the BC, we observe a range from 15 to 2229.9 for a total of 62 (63.2%) of the nodes.
The nodes with the highest levels of BC were STAT1 (BC=2229.9), IRF1 (BC=997.3), LYN (BC=796.9), and GZMB (BC=781.8). The genes encoding the above proteins were among those over-expressed in the AMR group. Of the proteins encoded by genes under-expressed in AMR, only one was represented in the PPI network as a central gene, ATP5B (GI=2, BC=153.7).
The results showed that proteins with a high degree of GI also show high values of BC, which suggests that, possibly, they are the essential proteins in regulating intracellular signaling. Therefore, they can be essential in the search for therapeutic targets or biomarkers. A. over-expressed genes in the AMR group (n=340)
The modules of the PPI network are densely connected areas of the network, which are considered relatively independent since they often work together to perform a biological function. Therefore, the modules of the PPI networks usually correspond to metabolic pathways or protein complexes. In our network, we obtained a total of seven modules ( Figure 3), although module 6 was not significant (p=0.233) and therefore was not considered for the functional analysis. The proteins of each module are shown in Table 4. The graphic representation of the protein connections of each of the modules is shown in Figure 4.
Next, a functional analysis was performed with terms from the KEGG database to obtain the main BPs in each module. The results of the functional analysis of the PPI network modules are shown in Table 5. The analysis shows the strong relationship of the modules with effector functions of the immune system such as antigen presentation (Module 0, p=0.0002), NK cell-mediated cytotoxicity (Module 1, p=0.0048), phagocytosis mediated by Fcg receptors (Module 3, p=0.0017), and trans-endothelial migration of leukocytes (Module 4, p<0.0001).

Identification of DEG-Associated Transcription Factors and Protein Kinases
Regarding protein kinases, in Table S3, we observed that the most important in the over-expressed genes are MAPK3, CSNK2A1, MAPK14, ERK1, and MAPK1. On the other hand, for the underexpressed genes, the essential protein kinases were DNAPK, CDK1, MAPK14, ERK1, and MAPK1. The MAPK14, MAPK1, and ERK1 protein kinases appear important for over-expressed and underexpressed gene signaling pathways, revealing commonalities upstream of the biological pathways activated during rejection.
For the transcription factors, the most important ones associated with over-expressed genes in the AMR group are SPI, IRF8, RUNX1, RELA, and IRF1, while for under-expressed genes, they were MYC, KLF4, USF1, SOX2, and TCF3 (Table  S4), showing a clear divergence in pathways between over-and under-expressed genes. The transcription factors SPI and RUNX1 are the ones that target the highest number of DEGs, with 17.9% (61/340) and 17.3% (59/340), respectively, of the total over-expressed DEGs. Regarding the under-expressed DEGs in the AMR group, the transcription factor USF1 targets 23% (6/26) and KLF4 19.2% (5/26) of the total.

The Cellular Infiltrate in the Graft With Antibody-Mediated Rejection
To evaluate the abundance of immune cells infiltrated in the graft during an AMR, we used the xCell tool (13). This tool is capable of estimating the relative abundance of 64 cell types from transcriptomic data. In our study, we focused on analyzing the 33 immune-type subpopulations that xCell can infer (Figures 5,  6). The transcriptomic data of the cohort GSE36059, composed of 281 biopsies of transplants NR and 65 of transplants with AMR, were entered into xCell. The results obtained of GSE36059 and the rest of the cohorts analyzed were similar are summarized in Table 6. These genes were entered in the Gene2Drug web tool to infer those molecules and drugs with the most remarkable capacity to decrease the expression levels of genes over-expressed in AMR ( Figure 7). Gene2Drug searches the CMap database, made up of more than 7,000 drug-induced gene expression patterns and low molecular weight molecules. Gene2Drug assigns each drug an ES (Enrichment Score) score ranging from -1 to +1. Those drugs that induce alterations in gene expression similar to those observed in rejection will have ES values close to 1, while producing an alteration in the opposite direction will have values close to -1. For our purpose, we will choose those drugs with ES values closest to -1 since they will be those that best counteract the expression levels of the altered genes in the AMR. The compounds with statistical significance are shown in     Table 7, and only imatinib is currently approved by the Food and Drug Administration (FDA) and the European Medicines Agency (EMA); that is why we decided to deepen the interactions of this drug. Imatinib is an inhibitor of kinases such as the PDGFRA and PDGFRB receptors, currently used to treat some cancers such as chronic myeloid leukemia or gastrointestinal tumors. As we see in Table 7, imatinib obtains better ES scores than other drugs commonly used in KT patients, such as Tacrolimus (ES= -0.227), Mycophenolic Acid (ES= -0.180), and Methylprednisolone (ES=0.173). Figure 8 shows that imatinib acts directly on CSF1R, a receptor for colony-stimulating factor 1, overexpressed in AMR. Furthermore, indirectly, it can alter the expression of STAT1 and LYN through blocking PDGFRA and PDGFRB, which, as we have seen previously, are central proteins that regulate the expression of a large number of genes involved in AMR. Therefore, imatinib may be a promising drug for the prevention or treatment of AMR.

DISCUSSION
The diagnosis of graft rejection in KT patients is made by evaluating the histological characteristics of biopsy samples. The evolution of omics sciences and bioinformatic techniques has contributed to the development of new diagnostic approaches by combining histopathological characteristics with the gene expression profile of the kidney graft, thereby achieving a more precise diagnosis. It is still necessary to improve our knowledge of the molecular mechanisms underlying kidney rejection to achieve non-invasive diagnostic methods or new therapeutic objectives that improve the clinical management of transplanted patients. For this reason, the central genes and biological pathways that regulate the AMR process were studied to understand better the molecular mechanisms that control it and propose possible biomarkers and new therapeutic options.
Transcriptomics studies have proven to be a valuable tool for understanding the molecular processes that govern graft rejection. One of the main limitations of these studies is, on the one hand, the small number of individuals per study due to the high cost of these techniques. On the other, the variability between different analyzes platforms, challenging to interpret and compare results between different research groups.
Although in some biomedical fields, such as oncology, transcriptomics studies are quite widespread, in transplantation, there are still few investigations with this methodology and, for this reason, a meta-analysis on gene expression is scarce in the literature. In transplantation, the meta-analyses of gene expression carried out have focused mainly on the search for expression profiles associated with graft rejection and tolerance (3,23), in addition to the search for therapeutic objectives, such as for kidney graft fibrosis (24).
In the context of KT, few studies have focused on studying the specific expression patterns of AMR. Kim et al. (25) performed a meta-analysis where they obtained that CXCL10, CXCL9, and GBP1 were the most over-expressed genes with AMR. Unlike the study carried out by Kim et al. (25), our work has focused exclusively on samples derived from biopsies, excluding those obtained from blood and PBMCs, intending to increase the homogeneity of the study, since the expression profiles in blood and graft may differ, leading to discrepancies between studies, making them difficult to compare.
Our study consisted of four cohorts with 466 renal biopsy samples, 137 from transplants with AMR and 329 from transplants with good graft function. The results showed 340 over-expressed genes in the AMR group, while only 26 were under-expressed. DEGs were functionally classified based on the terms GO for annotation of MFs and BPs, as well as based on the biological pathways of the KEGG.
Genes over-expressed in the AMR group were associated with BPs such as immune response activation, immune effector processes, leukocyte-mediated immunity, and molecular functions such as MHC binding or CXCR3 receptor binding. KEGG metabolic pathway analyzes showed biological pathways directly related to transplantation, such as GVHD and allograft rejection. Increased expression of genes related to MHC binding such as HLA, CD74, or TAP1 genes reflects the antigenpresenting solid activity in the allograft. It is known that in rejection processes, there is an increase in the maturation of antigen-presenting cells (APCs) such as DCs (23), which lead to an increase in the expression of genes related to antigen processing pathways such as CD74 or TAP1 (3,26).
On the other hand, we have obtained an increase in genes related to the CXCR3 receptor. CXCR3 is a chemokine receptor that polarizes CD4 T lymphocytes towards Th1, promoting infiltration of cytotoxic lymphocytes into the graft. Blocking this receptor is effective in inhibiting macrophage infiltration in acute rejection (27). However, some studies indicate that the infiltration of cells that express CXCR3 is only increased in celltype rejections and not in humoral ones (28). In our study, the over-expression of the CXCR3 gene was only obtained in two of the cohorts (GSE44131 and GSE93658), which indicates that it is not a good biomarker for AMR.
In contrast, the CXCL9 and CXCL10 genes, which code for CXCR3 ligands, were over-expressed in all cohorts in our analyzes, thus being better indicators of rejection. A possible explanation for this discrepancy may lie in the different abundance of receptor and ligand-producing cells into the graft. While the expression of the CXCR3 receptor is mainly limited to cells of the immune system such as lymphocytes and NK cells (29), the CXCL9 and CXCL10 ligands can be produced in large quantities by epithelial cells of the renal tubule in response to certain stimuli such as IFN-g (30). Previous studies have shown that the CXCL9 and CXCL10 genes' over-expression was associated with rejection processes, regardless of the tissue's type of rejection and origin (3,31,32). In heart transplant models, CXCL10 inhibition slows rejection and increases graft survival (33) The data from our analyzes confirm those obtained by other researchers, showing that the ligands that regulate the activation of the CXCR3 receptor are a potential biomarker and therapeutic target. In addition, among the over-expressed genes, we also obtained other biological pathways unrelated to transplantation, such as infection by Leishmaniasis and Staphylococcus aureus. These biological pathways are present because the genes involved in the immune response to the allograft are common to those involved in response to infections, such as the HLA or receptor genes of the immunoglobulin Fc fragment (30,34).
Among the under-expressed genes in the AMR group, the BPs related to the oxidation of fatty acids stand out: catabolic processes of fatty acids, monocarboxylic acids, and lipid oxidation. These results are compatible with those obtained in the analyzes of biological routes of the KEGG where two significantly decreased routes were obtained: degradation of fatty acids and degradation of isoleucine, leucine, and valine, although the differences were not significant after correction for multiple comparisons, probably due to the low number of underexpressed genes obtained in the analyzes. During the activation of the immune system, the proliferation and differentiation processes of lymphocytes lead to metabolic reprogramming. After the lymphocyte activation, there is an inhibition of boxidation and an increase in lipid biosynthesis, necessary for cell division (35,36). These observations are compatible with the results obtained in our study. During allograft rejection, a robust immune response is triggered with increased leukocyte proliferation. This can only be sustained if the metabolic pathways involved in biosynthetic processes are favored, thereby inhibiting the catabolic pathways of essential elements such as lipids and amino acids. Inhibition of lipid synthesis has been shown to interfere with lymphocyte proliferation and slow rejection (37). However, some studies indicate that the microenvironment conditions determine the metabolic status of the effector lymphocytes. In murine models of GVHD, it has been observed that fatty acid oxidation is necessary for effector cell function (38), although most research indicates that carbohydrates remain the primary fuel in immune effector processes (39,40). This fact becomes evident after immune activation, where the activation of mTORC1 promotes the increase of transporters such as Glut1 to increase the entry of glucose into the cell, necessary to sustain the increase in glycolytic activity. In addition, there is also   an increase in other transporters, such as leucine, during lymphocyte activation. Our data suggest that the inhibition of amino acid catabolism, together with the increase in amino acid transporters after immune activation observed in other studies, is a cellular response to maintain the intense biosynthetic activity that occurs during rejection (40, 41). Leucine antagonists or their transporters exert effects similar to those caused after inhibition of mTORC1, altering the differentiation of T cells (41,42). Our results show that during the processes that mediate AMR, there are changes in immune metabolism compatible with the processes of proliferation and activation of the immune system, and suggest that the pharmacological alteration of the catabolism pathways of some amino acids such as leucine could be a new therapeutic target for the inhibition of lymphocyte proliferation. Protein-protein interaction networks (PPI) were built to evaluate the relationship between the proteins produced by DEGs and to be able to deduce which are the most important core proteins that control information pathways at the molecular level. The PPI network obtained shows that four proteins control the main biological pathways related to AMR: STAT1, IRF1, LYN, and IRF8. STAT1 is a transcription factor that regulates the expression of genes for cellular response to IFNs and is closely related to the differentiation of T cells towards Th1 (43). The expression of STAT1 has been associated with inflammatory processes in KT related to ischemia-reperfusion damage and BK virus nephropathy (44,45). Furthermore, STAT1 has been linked to rejection processes in multiple organs in studies conducted in humans (46) and murine models (47). Although STAT1 expression was initially associated with cell-type rejections given its role in Th1 polarization, scientific evidence shows that IFN-related signaling pathways are intimately involved with AMR. That activation of cytotoxic-type effector functions also occurs with intensity in AMR (48). Based on these results, it has been suggested that STAT1 inhibition could be beneficial in  slowing down the immune response to graft. Blocking STAT1 with oligonucleotides in murine models of heart transplantation attenuates the recipient's immune response against the graft, promoting anti-inflammatory effects by inhibiting the expression of MCP-1 (49). Similar results have been obtained in murine models of bone marrow transplantation, where STAT1-deficient mice had decreased IFN-g production capacity, blocking T cell polarization towards Th1 and reducing rejection events (50). Furthermore, the STAT1 blockade appears to have more profound effects since it can promote the expansion of Treg (51). However, in other studies, it has been observed that the signaling pathways activated by IFN-g are also necessary to maintain Treg function, enhancing their ability to prevent rejections. More studies are necessary to understand the true impact of STAT1 blockade for the treatment of rejection.
In this sense, IFN-inducible factors (IRFs, Interferon Regulatory Factor) are transcription factors expressed in response to various stimuli, such as IFN, and have immunoregulatory properties (52). IRF1 is expressed in various organs and various types of immune cells such as NK cells and DCs. IRF1 increases the expression of proinflammatory cytokines such as TNF-a or IL-2 and chemokines such as CXCL10 (30,53,54). IRF1 has been shown to act as an early proinflammatory signal during ischemic damage in murine liver transplant models, promoting the secretion of cytokines, such as TNF-a, IFN-g, or IL-15 (55,56). Disruption of IRF-1 has been shown to decrease MHC expression in Knock-out (KO) pattern, although all-response rates remained similar to wild-type (57). IRF-1 mediates resistance to necrosis during rejection events, an effect associated with a decrease in IFN-g (4). In transplant rejection models, the decrease in IRF-1 has been associated with a decrease in the rejection event (58). Other studies carried out with gene expression data obtained from the GEO database have already indicated that the IRF1/STAT1 pathway could be important in rejection processes (46).
On the other hand, the IRF1/STAT1 pathway has also been associated with mechanisms that mediate tolerance in liver transplantation models through the activation of apoptosis (59). On the other hand, IRF8 has a close relationship with B cells differentiation and the secretion of proinflammatory cytokines by the innate immune system (60,61). IRF8 decreases its expression in liver rejection, although it is a study with a low number of samples (62). In KT, the function of IRF8 has been poorly studied. In a microarray study, IRF8 overexpression has been associated with kidney graft loss in those patients without rejection (63).
Additionally, LYN is a protein kinase of the Src kinase family expressed mainly in hematopoietic cells and played an essential role in regulating the responses of the innate and adaptive immune systems. The study of LYN in the context of KT has not yet been addressed, although the relationship between IFN expression and LYN activation could indicate an important role of the latter in rejection. As we have said previously, IFN pathways appear to be especially active during rejection processes. Plasmacytoid DCs and type I INF-producing APCs increase LYN expression during their activation, acting as a positive regulator in the induction of IFN and the production of cytokines (64). On the other hand, some studies indicate that LYN is involved in the signaling pathways involved in the expression of immunoglobulin Fc receptors such as CD16, or the expression of PIK3CG, a regulator of NK chemotaxis and cytotoxicity, suggesting that the increase in LYN expression could be associated with increased antibody-mediated cytotoxicity (65). Some studies indicate that LYN blockade could prevent the cytotoxic action of NK cells (66). In addition to DCs and NK cells, LYN is also expressed in B cells and is involved in the initiation of BCR signaling and the negative feedback responsible for signal inhibition (67). Inhibitors of this kinase are currently being tested to treat some diseases such as chronic B lymphocytic leukemia (68). Therefore, it could be a possible therapeutic target in transplantation, given its relevant role in activating the B cell and the effector functions of NK cells, although the redundancy of function of the Src family kinases, to which LYN belongs, could limit their activities. In addition, LYN-deficient murine studies produce myeloid hyper-activation, with increased BAFF secretion and a heightened inflammatory state (69). The increase in BAFF has been associated with an increase in the production rate of anti-HLA antibodies. Therefore, the LYN blockade could be counterproductive in the context of KT and should be further studied.
The analysis of transcription factors and protein kinases reveals the importance of IRF1 and IRF8 in regulating genes associated with AMR, like the analysis of PPI networks. However, it also shows the importance of other transcription factors such as RUNX1, RELA, and SPI1. These transcription factors are closely related to the differentiation, proliferation, and activation of cells of the immune system. First, RUNX1 is a hematopoiesis controller that controls the expression of essential genes such as NF-kB and, therefore, inflammatory responses. Second, RELA is the p65 subunit of NF-kB and is essential for its activation (70). SPI1 is a transcription factor of the Ets family expressed in B cells and myeloid cells, crucial in B differentiation and neutrophil activity (71). Bioinformatic analyzes have already reported on the relevance of SPI1 in the regulation of AR processes (70,72). Thus, PU.1 is a transcription factor encoded by the SPI1 gene and is vital in differentiating and developing macrophages and B cells (73).
Regarding protein kinases, the MAPK family shows a central role in controlling metabolic pathways of both over-and underexpressed genes. The analysis shows that MAPK1 and MAPK14 are important in the expression of genes of both groups. Among the kinases, it is worth mentioning ERK1, a serine-threonine kinase that participates in the Ras-MEK-ERK signaling cascade, involved in multiple processes such as cell adhesion, cell cycle progression, proliferation, and cell survival (72).
During AMR, vasculopathy occurs at the level of the endothelium of the graft vessels. The binding of anti-HLA class I antibodies to the endothelium increases endothelial, smooth muscle proliferation through an ERK1 process, also demonstrated by the generation of anti-non-HLA antibodies such as anti-AT1R (74). The pharmacological blockade of the MEK-ERK pathway effectively reduces chronic allograft nephropathy and the immune response in murine models of KT (75). In cardiac models, ERK inhibition was shown to inhibit the alloresponse, decreasing the leukocyte infiltrate to the graft and IFN production (76). Therefore, ERK1 kinase could be a potential therapeutic target to slow down antibody damage.
The xCell analyzes based on the inference of immune system cells from gene expression data show the important association of NK cell and monocyte infiltration with AMR processes. Therefore, inference data from cell subpopulations infiltrating the graft indicate that NK cells and monocytes are major cellular mediators of graft damage. Therefore, the STAT/IRF molecular pathway acquires great importance since it is directly related to the formation of cytokines, such as IFN and IL-15, or chemistries such as CXCL10, essential for the activation of these cells and their migration towards the graft.
Similar results have been obtained in heart transplantation, where NK cells and monocytes were also associated with AMR (77). Furthermore, several studies carried out in animal transplantation models corroborate the role of these cells in AMR (48,78). After transplantation, the function of NK cells is exerted through both antibody-dependent mechanisms, such as antibody-dependent cellular cytotoxicity, and non-antibodydependent mechanisms, through activation by ligands such as osteopontin (79). Although the significance is lost after statistical correction for multiple comparisons, we have observed a clear tendency for patients with AMR to have elevated levels of mast cells. Although mast cells have not been directly related to rejection processes, numerous studies link high levels with fibrosis and chronic graft damage (80,81). However, its role in transplantation remains controversial since high levels of mast cell-related transcripts have also been found in grafts from tolerant patients (82).
On the other hand, the levels of immature DCs were decreased in patients with AMR. It is known that, during rejection, there is an increase in the antigenic presentation of alloantigens and the maturation and infiltration of APCs (23). Therefore, the maturation of DCs would explain the low relative frequency of immature DCs obtained in our analyzes. In transplantation animal models, inhibition of its maturation has already been shown to suppress graft rejection (83). Therefore, NK cells, monocytes, and DCs are postulated as candidates for their use as biomarkers of AMR. The inhibition of these cells' maturation and/or the effector functions could be potential therapeutic targets to stop the immune response against the graft effectively.
Despite efforts in searching for new drugs in transplantation, it is still necessary to search for practical strategies for treating or preventing rejection events. In our study, the DEGs derived from the meta-analysis were integrated into bioinformatics applications to suggest new therapeutic options for AMR prevention. Our objective was to find new drugs capable of reversing the expression of the main genes altered during AMR. Of the candidate drugs, only imatinib is currently FDA approved. Imatinib is a tyrosine kinase inhibitor used primarily to treat chronic myeloid leukemia and gastrointestinal tumors (84). With a view to possible use in KT, different studies in murine models have shown the renoprotective effects of imatinib in the long and short term, limiting the progression of glomerulonephritis the prevention of chronic graft nephropathy (85,86). However, it has not been tested for this purpose in humans. In several clinical cases, imatinib has been used to treat chronic myeloid leukemia after KT has shown to be a safe drug without major side effects (87,88).
On the other hand, piperlongumine (PL), a natural alkylamide obtained from the Piper longum plant, caught our attention, which has recently attracted the attention of researchers due to its immunomodulatory properties. PL is effective in inhibiting STAT3 and inhibiting proliferation progression in breast cancer cell lines (89). It has also been observed to decrease the proliferation and survival of cell lines derived from hematological neoplasms by decreasing the expression of crucial transcription factors such as NK-kB, STAT1, or MYC (90,91). In autoimmune diseases, PL has also been observed to suppress DC maturation, decrease proinflammatory cytokines, and inhibit NF-kB (92,93). Therefore, although PL is still in preclinical studies, its immunomodulatory capacities postulate it as a potential drug for use in transplantation. It should be mentioned that, although our study does not cover the full spectrum of available drugs, limiting themselves to 1309 included in the current version of CMap (18), the strategy used is beneficial for identifying new indications for existing drugs. Future studies in animal models and in vitro will be necessary to confirm the hypotheses arising in this study.
Our study has several limitations; in the first place, the strategy of identifying DEGs by selecting those statistically significant in all the selected GEO studies generates a bias as they are studies with a very uneven number of samples. While the GSE36059 study included 346 samples, the GSE44131 study only 23 samples, so the statistical power of this study will be lower and will act as a limiting factor when identifying DGEs and other genes of interest that may have been lost. Therefore, other statistical methods are necessary to integrate and normalize the data from different studies. Another aspect of assessing is the type of platform used in the studies. In our case, all the studies were carried out using Affymetrix microarrays. Although the GSE36059 was carried out under another version, the same platform allows the data to be more homogeneous and comparable.
The low number of under-expressed genes in the AMR group implies that the functional and topological bioinformatics analyzes have lower power, so the conclusions should be assumed with caution. Finally, inference studies of the immune populations in the graft show the relevant role of NK cells and monocytes, but these algorithms do not discriminate between cell subtypes. It is known that NK cells, in addition to their cytotoxic function, also possess immunoregulatory functions through the secretion of cytokines and chemokines. Therefore, using other tools such as flow cytometry is necessary to identify further the phenotype of infiltrating NKs in the graft. Another of the shortcomings of the study has been not to include samples from cellular rejections. Future studies with cellular rejection samples will be necessary to determine which expression profiles are common to rejection processes and differentiate between subtypes. Furthermore, the selected GEO studies do not provide information on the acute or chronic nature of the rejection, which means that the analyzes are based on a heterogeneous cohort regarding the post-transplant period in which graft damage occurs. On the other hand, only samples taken at the time of rejection have been studied. In the future, therefore, it would be interesting to carry out longitudinal studies that make it possible to correlate gene expression profiles in the early stages of transplantation with chronic damage, as other previous own studies (94,95), which allow us to develop prognostic models for the better classification of patients and thereby improve their clinical management. Finally, the results obtained in silico must be confirmed by in vivo and in vitro analysis in the future.
In conclusion, in the present study, we have been able to identify the gene expression profile associated with AMR and the central infiltrating leukocyte populations in the allograft through different bioinformatics tools. The functional analysis shows how the processes related to immune activation and antigenic presentation are fundamental in AMR mechanisms. Furthermore, topological analyses using PPI networks show the importance of the IRF/ STAT1 pathways of response to IFN in controlling the expression of genes related to humoral rejection. On the other hand, inference analyzes of leukocyte populations show a crucial role of NK cells and monocytes in graft damage during rejection. Future studies will be necessary to validate the results obtained and assess the potential use of over-expressed genes as rejection biomarkers. In addition, the results suggest that the central proteins obtained in the PPI networks could be potential therapeutic targets to improve the results of kidney transplantation.

AUTHOR CONTRIBUTIONS
RA, MM, HM-B, and IL participated in designing and contributed to data processing and supervision, optimization, data analysis, and manuscript writing. IL and RA participated in data analysis for gene expression assays and contributed to writing the manuscript. VJ-C, HM-B, JG, CB, MM-Q, and AP participated in the data processing. RA, SL, and MM-P contributed to gene expression data generation, organization, and manuscript writing. AM provided the study samples, participated in discussions in data analysis strategies. All authors contributed to the article and approved the submitted version.

FUNDING
Our work was possible thanks to support from Instituto de Salud Carlos III (ISCIII), Spanish Ministry of Economy and Competitiveness. Grant Number PI15/01370, P19/01194, and PI20/00050 and co-funding the European Union with the European Fund of Regional Development (FEDER) with the principle of "A manner to build Europe".