DIAMOnD Genes Expression in Human Macrophages Infected with Leishmania major

1. Introduction

Leishmaniasis is a spectrum of vector-borne parasitic diseases caused by unicellular protozoa of the genus Leishmania. They remain a significant public health burden, with an estimated 12 million infected individuals worldwide and approximately 350 million people at risk [1]. Despite extensive research efforts and accumulated knowledge, leish-maniasis continues to be classified as an emerging or re-emerging disease with complex control challenges. In response to the ongoing challenges and global needs, the scientific research on leishmaniasis has focused on acquiring new knowledge and developing effective control measures [2]. The intricate immune response to Leishmania infection remains incompletely understood, but the parasite’s remarkable ability to outsmart host immune cells is a key factor in its pathogenicity [3]. Leishmania’s strategies extend beyond mere survival within macrophages. Their interaction triggers the release of antigenic particles, stimulating the host’s adaptive immune response [4]. Macrophages, however, stand as the parasite’s preferred sanctuary while play a central role in shaping the immune response, highlighting their critical role in leishmaniasis pathogenesis. These phagocytic cells serve as the parasite’s primary target and battleground. Notably, L. major is frequently utilized in experimental models, providing well-characterized in vivo and in vitro insights into disease development [5]. Deciphering this complex interplay, including the role of different macrophage phenotypes and their interactions with Leishmania holds immense promise for developing novel therapeutic strategies against this persistent parasitic dis- ease. Leishmania infection and the parasite’s ability to resist the host immune response involve changes in gene expression in the cells that harbor the parasite and in the parasite itself. Within a genome, well-defined gene expression control mechanisms regulate those genes that must be expressed by the cell at a given time in response to internal or external factors [6,7,8]. Gene expression profiling has been used in several studies on pathogenic microorganisms, including parasitic protozoa such as Leishmania [6,9]. This technique has allowed researchers to identify differentially expressed genes at distinct stages of infection, leading to a better understanding of the sequence of gene activation and its association with critical events at each stage. Changes in the expression of genes encoding STAT2, IL-18, and CXCL2, as well as genes important in the steroid and cholesterol biosynthesis pathway, have already been identified in infected cells [9]. As in mice, the gene expression of human macrophages can be differentiated into a susceptibility or resistance profiles when they become infected with Leishmania major [10]. In this regard, the study of the interactome, i.e., the whole set of protein-protein interactions (PPI) in a specific cell or organism, would be a valuable resource for understanding the molecular basis of Leish- maniasis, identifying new targets for developing drugs and designing synthetic biology circuits. However, assembling the interactome is a daunting task, as it requires identifying and characterizing all protein-protein interactions in a given system [11]. Recent advances in high-throughput technologies has uncovered interesting aspects of leishmania infection and details of the underlying protein-protein interactions have brought us closer to the goal of assembling the interactome, a comprehensive map of all protein-protein interactions in a cell or organism [12,13]. However, there are still many challenges to overcome, such as the need to validate protein-protein interactions, the abundance of false negatives and to understand the functional significance of these interactions. Despite these challenges, the interactome is a goal of inestimable value with potential to revolutionize our under- standing of biology and medicine. These methodologies have converged in the search for hypotheses that can be tested in silico, in vitro, and finally in vivo in clinical studies that can support the research data [14]. In view of the increasing availability of data on leishmaniasis, the present work was focused on performing data mining of public databases that provide a wide spectrum of different results and methodologies. The goal was to search for results on leishmaniasis in general, leading to filtering and treatment of the data to assemble a protein-protein interaction network and consequent analysis of the centrality and importance of the genes in the network. To confirm these in silico-generated data, some genes were then selected for biological validation. The parasite-host interaction of leishmaniasis has been studied by many models, including computational models that describe the dynamics of the parasite and macrophages in the early phase of the immune response [15]. Recent findings in the field of data mining have also shown the power of this approach [16]. However, no work has yet been done on database integration mining for leishmaniasis, as in this study, in which we looked for information about leishmaniasis in several databases.

2. Results

2.1. Data Sets, Gene Identification and Enrichment

We utilized the following databases: GWAS, UniProt, Malacards, DisGeNET, Hu-manMine, NCBI, and GeO from cells infected with Leishmania, alongside a search for the most cited proteins produced upon leishmania infection relevant articles. This comprehen- sive approach allowed us to amass a large number of proteins and genes associated with leishmaniasis. We conducted a manual search of the known literature. In all databases, searches were made using the keyword "Leishmaniasis". We then obtained proteins and genes which were organized only with human datasets (HDs) according to the origin of the data, the reference, and the proteins found in the results (Supplementary Table S1). The search and assembly of the HDs generated 439 proteins and they had identification codes (ID). After checking the protein coded by each gene, gene symbols were identified and confirmed using UniProt and HGNC databases. The duplicate genes, proteins that were not associated with a gene symbol, and proteins from different databases/studies were filtered out. After filtering, we were left with 287 genes associated with the 439 proteins. The identification procedures and filtering resulted in a high-quality dataset of genes and proteins related to leishmaniasis. After identifying 287 genes (Supplementary Table S2), they were used in the Cytoscape system, for the assembly of PPI networks and gene ontology enrichment analysis. The gene set enrichment analysis returned pathways related to im- mune response. To go deeper into the biological significance and functional relationships of the 287 genes identified in our initial analysis, we performed an enrichment analysis. This revealed an intriguing association of the identified genes with specific biological processes1 (Figure 1) and with leishmaniasis (Supplementary Figure 1). Notably, significant enrichment was observed in pathways related to cytokine signaling. Furthermore, the gene set exhibited enrichment in terms related to broader biotic interactions, such as "response to other organisms" and "response to external biotic stimuli" (Figure 1A). These findings suggest a potential role for these genes in processes involving inflammatory signaling and inter-organismal interactions. To further illuminate the overarching themes within the enriched pathways, we employed a network-based analysis of Gene Ontology (GO) terms, and this facilitates the identification of overarching themes that are not readily apparent from individual pathway enrichments (Figure 1B,D,F). The examination of cellular components reveal the crucial role of the external plasma membrane’s, the cell surface itself, and their combined interaction in constituting key signaling pathways within source genes (Figure 1C,D). Scrutiny of molecular functions within our source genes reveals a striking enrichment for terms associated with “cytokine receptor interactions”, “activity”, and “signaling pathways”. This pattern suggests the emergence of cytokine signaling as a critical overarching functional module within these genes (Figure 1E,F). These findings underscore the pivotal role of cytokine signaling in shaping the overall functional landscape of genes expressed by leishmania-infected cells.

2.2. Assembly of Interactome for DIAMOnD Genes Identification

To identify the relationship between the 287 genes related to the 439 proteins under scrutiny, we used Cytoscape [17] to build the PPI network using the plugins APID2NET [18] and BioGrid [19]. We then obtained two complex interactomics corresponding to the module of our interest originating from two different plugins, BioGrid and APID [20]. The BioGrid interactome was assembled by importing PPI data from the BioGrid database. This interactome contains 259 nodes and 1727 edges, demonstrating the high linkage and high number of interactions between the source genes (Figure 2A). The APID interactome was assembled by importing data from various sources, including PPI databases, expression datasets, and GO annotations. This interactome contains 253 nodes and 418 edges, showing the interactions between the proteins encoded by the source genes (Figure 2B). The two interactomes we assembled provide complementary views of the PPI network between the source genes. The BioGrid interactome is more comprehensive regarding the number of PPIs, while the APID interactome includes additional information, such as expression data and GO annotations. Combining these two interactomes provided a more complete view of the PPI network and its relationship to other molecular states.

2.3. Identification of DIAMOnD Genes

To investigate the connectivity patterns of protein interactions, we leveraged the well- established connection of the plug-ins with databases of protein interactions. We have opted to employ the DIAMOnD [21] predictive algorithm. This algorithm facilitates the assessment of the significance of the connectivity between immediate neighbors of proteins by iteratively expanding the module; it then adds the most significantly connected node to the module increasing the module by one node at a time. We used DIAMOnD to identify special modules of genes that are involved in leishmaniasis, namely, modules that are highly interconnected in the two assembled interatomics (APID and BioGrid interactomes). The output of DIAMOnD was selected to be 300 genes, which is compatible with the literature stated by Ghiassian19. This number is reasonable considering that our original dataset of interest contained 287 genes. However, the two results did share some common genes. Applying DIAMOnD to the genes derived from the BioGrid plug-in resulted in an increased number of genes and edges within the network, resulting in 300 nodes and 11,741 edges (Figure 2C). This led to the formation of a denser interactome. Similarly, applying DIAMOnD to the APID-derived interactome also yielded an augmented number of nodes and edges, identifying 300 nodes and 2,929 edges (Figure 2D). When we overlapped the results of the 300 genes obtained from both interactomes, we identified 217 genes with higher centrality in both the APID and BioGrid interactomes. We dubbed them "Diamond genes". The two plug-ins are different from each other, but not so much, which makes the result of 217 genes (Supplementary Table S2) compatible with what was expected. Further study of the Diamond genes is needed to understand their role in leishmaniasis.

2.4. The DIAMOnD Genes Are Expressed by Human Macrophages Infected with Leishmania major

After identifying the DIAMOnD genes from both the APID and BioGrid plugins, we compared the categorizations generated by the results of the two plugins. APID DIAMOnD and BioGrid DIAMOnD generated ranked results based on the centrality values of these genes in the networks, generating a ranking of the 217 genes. The comparison of the results from the two plugins allowed us to select the genes of interest more precisely. The shared genes found in the DIAMOnD module suggest that DIAMOnD G genes are strong candidates for the category of central genes regulating the response to leishmaniasis. To confirm the importance and the centrality of these genes in the control of leishmaniasis (with the macrophage as the reference cell), further biological validation of the Diamond Genes as core genes was required. We then chose to select the 22 most prominent genes for biological validation (Supplementary Table S3). To perform a biological validation that was consistent with our gene selection validation, we used a model of in vitro human macrophage infection, THP-1 cells [22]. The parasite used was the Leishmania major, a Leishmania species most used in experimental infection [23,24]. After separation of the metacyclic infective forms of Leishmania major, THP-1 cells were infected with 10 parasites for each macrophage. THP-1 cells were effectively infected3 (Figure 3A,B). with the metacyclic parasites of Leishmania major, thus demonstrating that THP-1-derived macrophages are functionally hosts for the parasite [25]. To confirm the biological validation data, we analyzed first the expression of four cytokine genes related to innate immune response, well known to be produced by resistant macrophages to L. major (IL-1alfa, IL-1beta, TNF-alfa, and IL-6). The INF-gamma and LPS were used as a positive control for activated macrophages. As expected, the results showed that THP-1 cells stimulated with IFN-gamma and LPS for 24 hours had higher expression of inflammatory cytokines compared to no-infected cells. However, the THP-1 infected with L. major for 24 hours had no increase in the expression of those cytokines, demonstrating THP-1 cells as model for susceptibility to L. major infection3 (Figure 3C).

We opted to perform RT-qPCR experiments at 12 and 24 hours post-infection. The results obtained at 12 hours post L. major infection, show the expression of several of the "Diamond Genes" during the early stage of the infection. However, none of these genes exhibited significant differential expression compared to non-infected cells (either over- or under-expressed) at the assessed time point (Figure 4A). We then investigated the expression of the DIAMOnD genes in THP-1 cells infected with L. major at 24 hours post infection and found that these genes were still expressed, but at lower levels than those found after 12 hours of infection. No significant differences in the expression of the analyzed genes (neither over- nor under-expressed) were found between control and infected cells at this time point (Figure 4B). Upon analyzing the expression of the same 22 genes at 24 hpi under IFN-gamma plus LPS stimulation (Figure 4C), we confirmed the reactivity of THP-1 cells to the stimulus. We found statistically significant differences in expression for some genes (specifically, CASP8, IL4R, NCK2, Ripk2, TANK, and TNFRSF10A), compared to no stimulated cells and to L. major infected cells. Again, we observed that L. major can induce decreased gene expression in THP-1 cells, showing a profile compatible with susceptibility to leishmania infection. In short, the identified DIAMOnD genes can be a good target to modulate the immune response to leishmania infection.

3. Discussion

Leishmania is an obligate intracellular parasite that can subvert the host immune re- sponse in its favor affecting macrophage response mechanisms such as apoptosis, oxidative damage, subversion of the immune response, and modification of nutrient availability [26,27,28]. To evaluate the genes involved in the parasite-host relationship, new methodologies are needed. In the present study, we searched databases for genes related to leishmaniasis and then filtered the results only for those human genes associated with macrophages. We thus turned to understanding the parasite-host response from the point of view of gene expression by macrophages in an attempt to build an omics model of the parasite-host relationship. Similar relatively recent studies have shown satisfactory results in the liter- ature of data mining [29]. Due to the pivotal role of protein–protein interactions (PPI) in the pathogenesis of diseases, the manipulation of protein-protein complexes holds signifi- cant clinical relevance. The effective design of compounds influencing PPI necessitates a comprehensive understanding of the molecular intricacies of the involved protein-protein system. Although there has been progress in data mining and in the use of interactomes as study tools, the biological validation of these data remains crucial. When we started this step, it was expected that many DIAMOnD genes analyzed would be differentially expressed by infected cells when compared to uninfected cells. Surprisingly, the genes analyzed at both infection times showed the same behavior in infected and uninfected cells. This fact could be linked with the susceptibility to leishmania. The successful infection and ultimate dissemination of parasites throughout the body may depend, at least in part, on the quite early interactions between the parasite and the host cell [30]. Genotypic differences between the same Leishmania species can lead to detectable differences in macrophage responses [31]. However, the expression of certain genes in the macrophage in response to a pathogen is also described for some genes, among which, some included in our DIAMOnD genes list, such as TNFR1A. This gene, which encodes a cytokine receptor, is described to be much more expressed in macrophages when they are infected by a pathogen and are able to react [31]. Leishmania infection has already been well described and characterized in THP-1 cells. This model has been used for more than two decades [23,32]. THP-1 cells differentiate into human macrophages and are therefore a good model for the study of parasites that infect this cell type, such as Leishmania. These cells have also been used in studies on leishmanicidal drugs and inhibitors against Leishmania [32,33]. Another use of these cells was in studies that demonstrated modification in gene expression of proteins related to gene transcription, RNA editing, histones, DNA repair and replication, and also proteins involved in cell survival and signal transduction [34], in addition to proteomic and other omics studies even demonstrating their altered response to infection with different pathogens [22]. All these successful applications demonstrate the acceptance and practicality that this model offers. Therefore, susceptible and resistant macrophages show a different expression profile when infected with the same species of Leishmania10. The inherent transcriptional heterogeneity among different THP-1 macrophage cell lines in response to Leishmania infection, with consequent variations in functional pathway and gene mechanisms, provides valuable context for interpreting our current results [35]. This raises the questions as to when and how Leishmania species affect and modify host cell expression. Our results confirm that there is a large number of genes involved in this process, considering that the parasite-host interaction induces changes in gene expression in several models already studied [36]. L. major then seems to trigger normal gene expression mechanisms in THP-1 macrophages, i.e., L. major seems to induce a gene expression pattern in macrophages equivalent to the uninfected state, thus preventing the development of a strong immune response against the parasite. This suggests an escape mechanism favoring the survival of the parasite and spread even within a cell of the immune system [37]. The lack of differential expression in Diamond genes at both 12h and 24h suggests a potential Leishmania-mediated manipulation of host cell gene expression. This might involve main- taining baseline expression levels of crucial immune response genes to establish a favorable environment for parasite survival. The observed decrease in Diamond gene expression at 24h might indicate additional regulatory mechanisms employed by the parasite.

4. Conclusions

This study’s novel methodology yielded profound results. First, we displayed capacity for public data acquisition and processing, employing enrichment analysis and Cytoscape for comprehensive filtering and analysis. This approach culminated in identifying 217 " DIAMOnD Genes", exhibiting significant relevance to leishmaniasis. Subsequent biological validation revealed that L. major maintains basal expression of key immune response-related genes within the utilized model. This observation warrants further investigation into the parasite’s regulatory mechanisms governing the expression of these genes. Additionally, our findings suggest a highly efficacious intracellular escape mechanism employed by L. major, contributing to its persistence within host cells. This discovery opens exciting avenues for exploring the parasite’s immune evasion strategies and developing novel control measures.

5. Materials and Methods

5.1. Data source and Filtering

We selected the following databases: GWAS (Genome-Wide Association Studies), UniProt (Universal Protein Resource), Malacards (Leishmaniasis Knowledgebase), Dis- GeNet (Disease-Gene Network), Humanmine (Human Disease-associated Genes and Muta- tions), and NCBI (National Center for Biotechnology Information). To conduct our searches, the keyword "leishmaniasis" was employed, yielding results encompassing microarray and RNAseq data concerning gene expression, DNA methylation, miRNAs, and immune response. The data collected for analysis encompassed genes, proteins, and transcripts associated with “leishmaniasis” word in general. However, we exclusively focused on data about humans. By adhering to these criteria, we effectively selected results that align with our research objectives, subsequently delving into exploring expressed genes and networks induced by leishmaniasis in human macrophages. Tata from each selected result underwent filtration following the aforementioned criteria. Subsequently, the UniProt and HGNC databases were utilized to convert the results into protein and gene symbols.

5.2. Enrichment Analysis

This approach assesses whether the gene set is statistically enriched with genes belong- ing to specific pathways or functional categories, such as those defined by Gene Ontology (GO) terms [38]. To perform the enrichment analysis, we employed the ShinyGO [39] tool alloweing to interrogate the 287 identified genes and assessing their enrichment in Gene Ontology (GO) terms across the three main categories: Biological Process, Molecular Function, and Cellular Component. Within the Biological Process, we focused on terms describing ordered assemblies of molecular functions that contribute to the accomplishment of specific biological functions. We deemed them enriched when the GO term had a false discovery rate (FDR) 0.05. On enrichment analysis, the FDR was calculated based on the nominal p-value derived from the hypergeometric test. Fold enrichment captured the proportion of genes within your list belonging to a given pathway, expressed as a percentage relative to the background set. FDR quantified the probability of observing any given level of enrich- ment by chance, with larger pathways generally exhibiting lower FDRs due to enhanced statistical power. Fold enrichment, conversely, served as an effect size measure, indicating the degree to which genes in a particular pathway are overrepresented. Pathway enrich- ment analysis was restricted to a user-defined size range. Following analysis, pathways were first filtered based on a specified FDR threshold, and successively selected for the top pathways according to FDRs. Finally, highly redundant pathways sharing over 95% gene overlap were represented solely by the most statistically significant one. The network ap- proach employed by the ShinyGO tool, visually depicts the relationships between enriched pathways, with intersecting genes represented by links between pathway nodes. This network visualization allowed the identification of the most biologically relevant processes emerging from the enrichment analysis based on their interconnectedness and shared-gene content. As pathways with greater overlap in gene membership are linked more closely, GO pathways (nodes) are connected if they share 20% or more genes. Darker nodes are more significantly enriched gene sets. Bigger nodes represent larger gene sets. Thicker edges represent more overlapped genes.(Figure 1).

5.3. Network Analysis

To integrate and assemble PPIs, we used the Cytoscape [17] software. We chose two different plugins, APID2NET [18] and BioGRID [19], to measure and analyze the PPIs, inves- tigating the genes. The Apid (Agile Protein Interactomes DataServer) reports information on experimental validation, allowing selection and filtering at different quality levels for each protein-protein relationship. This provides access to the interactome of specific species, we focused on human data, more than 90,000 different proteins, and more than 670,000 unique interactions18. We also used the BioGRID Cytoscape’s plugin to investigate protein- protein interaction networks. It currently has more than 70,000 publications, more than 1 million protein and gene interactions, more than 20,000 chemical associations, and more than 800,000 post-transcriptional modifications from diverse species. By using these two plugins together, we meant to assemble an interactome with a large quantity and quality data resulting from our HDs. The DIAMOnD [21] (Disease Module Detection Algorithm) tool is used to identify "DIAMOnD genes" within interactome modules. DIAMOnD genes are genes that play a crucial role in the interactome and are highly likely to act as central modulators of the response within the interactome. The DIAMOnD algorithm operates by assessing the significance of connectivity among all immediate neighbors of the proteins at each step. The node with the most significant connectivity (i.e., the lowest p-value) is integrated into the module, thus expanding the module by one node at each refining step. In this manner, the DIAMOnD tool facilitates the identification of DIAMOnD genes that hold potential importance within the assembled PPI network and, consequently, within the host immune response process.

5.4. Biological Validation

In this study, we performed biological validation through cellular and molecular methodologies to corroborate the data obtained in silico. We cultured THP1 cells, a human monocyte lineage, as described in the literature [40,41,42]. THP1 cells were cultured in RPMI 1640 medium supplemented with 10% inactivated fetal bovine serum (FBS), 10 mM sodium pyruvate, 10 mM L-glutamine, and gentamicin at 37°C and 5% CO2. To induce macrophage differentiation, THP1 cells were plated in a 24-well plate covered by glass slides and were treated with 10 ng/mL phorbol 12-myristate 13-acetate (PMA) for 48 h. Successively, the adherent cells were washed to remove excess PMA and infected with metacyclic Leishmania major in a 10:1 ratio (10 parasites per macrophage). The infection was carried out at 32°C for 3 h. The cells were then washed to remove non-phagocytized Leishmania and the culture was maintained for 12 h or 24 h. Leishmania major (WHO MHOM/IL/80/Friedlin) was grown in Grace’s medium supplemented with 20% FBS, 2 mL glutamine, and gentamicin at 25°C. Metacyclic promastigote forms were separated by Ficoll gradient on day 5 of culture. In a 15 mL tube, 2 mL of 20% Ficoll was added, followed by 2 mL of 10% Ficoll. The Leishmania suspension was then added carefully to the wall of the tube, forming a third phase. The triphasic mixture was centrifuged at 800 g for 10 min at 4°C. The ring formed and all the supernatant above it, corresponding to the 10% Ficoll phases, was collected as it contained the metacyclic forms. The Leishmania was then washed three times and resuspended in complete RPMI. At 12 and 24 h post-infection (hpi), the cells were lysed with Trizol, and RNA was extracted following the manufacturer’s protocol (Invitrogen). RNA quality was assessed using the NanoDrop equipment and evaluating the 260/280 and 260/230 ratios. To evaluate the expression of the genes of interest, we performed quantitative reverse transcription PCR (RT-qPCR). We performed this technique in three steps: (i) conversion of RNA into complementary DNA (cDNA) through reverse transcription, (ii) amplification of cDNA using the polymerase chain reaction (PCR), and (iii) quantification of the amplification of gene expression products in real time [43]. For statistical analysis of RT-qPCR results, we used the Comparative delta-delta CT Method for Relative Quantification (ddCT). This method uses the arithmetic formula 2-CT to achieve the result for relative quantification. The amplification efficiency of the target of interest and the amplification efficiency of the reference (endogenous control) must be approximately equal for the comparative CT method to be validated. To perform RT-qPCR, we employed specific primers for the DIAMOnD Genes (Supplementary Table S3), which were designed based on information from the PrimerBank database44. We also used control (endogenous) expression genes, namely RNA-18s and GAPDH, as referenced in the work by Sikand et al. (2012). The Applied Biosystems 7900HT System equipment, along with Applied Biosystems SDSv2.4 software, was used for data visualization and preliminary analysis of the obtained results. The cellular lysis process involved TRIzol Reagent (Invitrogen) following the manufacturer’s protocol. After lysing THP-1 cells infected with Leishmania major, RNA isolation was performed using chloroform to separate RNA from phenol. Isopropanol was then used to isolate RNA, resuspended in RNAse-free water. The quality of the RNA was assessed using NanoDrop before proceeding with cDNA synthesis. For cDNA synthesis, the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher) was used to generate the necessary cDNA for RT-qPCR. To form the RT-qPCR mix, we used the PowerUp SYBR Green Master Mix (Thermo Fisher) and the primers (listed in Supplementary Table S4), following the manufacturer’s protocol. This mix, along with the samples, was plated on the MicroAmp Optical 96-Well Reaction Plate (Applied Biosystems) sealed using the MicroAmp Optical Adhesive Film (Applied Biosystems). The subsequent RT-qPCR was conducted with the 7900HT System (Applied Biosystems).