key: cord-0730213-ch5axwvv
authors: Martinez-Hernandez, J. Eduardo; Hammoud, Zaynab; de Sousa, Alessandra Mara; Kramer, Frank; do Monte-Neto, Rubens L.; Maracaja-Coutinho, Vinicius; Martin, Alberto J. M.
title: Network-Based Approaches Reveal Potential Therapeutic Targets for Host-Directed Antileishmanial Therapy Driving Drug Repurposing
date: 2021-10-20
journal: Microbiology spectrum
DOI: 10.1128/spectrum.01018-21
sha: 5e0f54861345c37bf965a22a0ce0eaea97f9ecf8
doc_id: 730213
cord_uid: ch5axwvv

Leishmania parasites are the causal agent of leishmaniasis, an endemic disease in more than 90 countries worldwide. Over the years, traditional approaches focused on the parasite when developing treatments against leishmaniasis. Despite numerous attempts, there is not yet a universal treatment, and those available have allowed for the appearance of resistance. Here, we propose and follow a host-directed approach that aims to overcome the current lack of treatment. Our approach identifies potential therapeutic targets in the host cell and proposes known drug interactions aiming to improve the immune response and to block the host machinery necessary for the survival of the parasite. We started analyzing transcription factor regulatory networks of macrophages infected with Leishmania major. Next, based on the regulatory dynamics of the infection and available gene expression profiles, we selected potential therapeutic target proteins. The function of these proteins was then analyzed following a multilayered network scheme in which we combined information on metabolic pathways with known drugs that have a direct connection with the activity carried out by these proteins. Using our approach, we were able to identify five host protein-coding gene products that are potential therapeutic targets for treating leishmaniasis. Moreover, from the 11 drugs known to interact with the function performed by these proteins, 3 have already been tested against this parasite, verifying in this way our novel methodology. More importantly, the remaining eight drugs previously employed to treat other diseases, remain as promising yet-untested antileishmanial therapies. IMPORTANCE This work opens a new path to fight parasites by targeting host molecular functions by repurposing available and approved drugs. We created a novel approach to identify key proteins involved in any biological process by combining gene regulatory networks and expression profiles. Once proteins have been selected, our approach employs a multilayered network methodology that relates proteins to functions to drugs that alter these functions. By applying our novel approach to macrophages during the Leishmania infection process, we both validated our work and found eight drugs already approved for use in humans that to the best of our knowledge were never employed to treat leishmaniasis, rendering our work as a new tool in the box available to the scientific community fighting parasites.

L eishmaniases are a group of vector-borne neglected tropical diseases caused by Leishmania parasites (1) . The clinical manifestations range from self-healing skin ulceration (cutaneous leishmaniasis [CL] ) to splenomegaly and hepatomegaly (visceral leishmaniasis [VL] , which is the deadliest form of leishmaniasis) (2) . According to the last report for leishmaniasis emitted by the World Health Organization (WHO; www .who.int/en/news-room/fact-sheets/detail/leishmaniasis), it is estimated that annually there are about 1 million infections for CL and 90,000 for VL. Despite the availability of some drugs for leishmaniasis treatment (3) , drawbacks in their current clinical use have been documented, ranging from high costs, toxicity, and the selection of resistant parasites (4) .

Host-directed therapies (HDTs) are a group of strategies that interfere in the host mechanisms that are necessary for pathogen survival and/or stimulate the immune response to respond to pathogens and eliminate them, bypassing existing limitations with conventional treatments, such as the chance of developing resistance (5, 6) . Recently, HDTs have been proposed for the treatment of diverse bacterial, viral, and parasitic diseases (7), such as tuberculosis (8) , malaria (9) , HIV infections (10, 11) , and most recently COVID-19, caused by SARS-CoV-2 (12, 13) .

To this day, a variety of strategies have been developed to identify new hostdirected therapies for leishmaniasis treatment (5, 14) . Many of these strategies are focused on the improvement of the immune response of the host (14) . In previous work, Murray and colleagues demonstrated that a combination of interleukin 12 and the typical treatment with pentavalent antimony (Sb v ) helped in the recovery of animals infected with Leishmania donovani, proving that a joint therapy between drugs that improved the host's immune response and conventional therapies can be useful for the parasite elimination (15) . Another study revealed that imatinib, an anticancer drug, was useful for reducing the severity of lesions caused by Leishmania amazonensis (16) . Other studies were focused on promoting the production of interferon gamma (IFN-g) (17) , nitric oxide (NO) and interleukin 12 (IL-12) (18) , and reactive oxygen species (ROS) (19) , resulting in the improvement of immune response and promoting healing.

The use of network-based approaches to determine nonobvious biological interactions and their relationship to disease has recently increased with the continuous evolution of high-throughput technologies for transcriptomics, proteomics, and metabolomics assays, as well as the generated data available in public databases (20) . These biological networks are a form of knowledge representation, used to structure different levels of relationship between bioentities (21) . Networks can simplify the complexity and heterogeneity of biological systems and contain a myriad of knowledge eager to be explored through different specific computational approaches. Examples of this complexity are the Gene Regulatory Networks (GRNs), which represent regulatory events between regulatory elements, such as those between transcription factors (TFs) or noncoding RNAs (ncRNAs) and protein-coding genes (22) . GRNs in combination with expression data can be used to infer condition-specific networks, which can be compared with control contexts or steady states, allowing the identification of possible disease markers or potential targets for a pharmacological treatment (23) . However, biological networks represent not only regulatory interactions but also interactions connected by different means, such as nondirected protein-protein interactions or directed metabolic reactions. These aspects are employed in many fields of biomedicine by combining different networks, such as protein-protein interactions (PPIs), GRNs, and epidemiology information (24) . The integration of the information available in the multitude of networks that one can obtain from a biological system is a complex task. To facilitate this process, multilayered graphs emerge as a solution recently introduced in the field of biomedicine (24) . This concept consists of layering the components of the network, nodes, or edges, grouping them by their type based on the heterogeneity of the network nodes or the relations among them. So, network elements of the same type are organized in the same layer and layers are connected by edges linking different types of elements (25) .

Here, we present a combined approach to determine potential therapeutic targets for host-directed antileishmanial therapies. First, we employed GRNs based on transcriptomic data to model Leishmania infection dynamics in human macrophages, which are followed by a multilayered networks approach to map metabolic/signaling pathways and drug-target identification. Through this approach, we identified five potential novel targets according to their direct connection with 11 known drugs that could be repurposed to be used in host-directed antileishmanial therapy.

Global expression patterns in Leishmania-infected macrophages. A global gene expression analysis was performed using a publicly available set of transcriptome sequencing (RNA-seq) data comprising four time points (4, 24, 48, 72 hours postinfection [hpi]) of Leishmania major-infected and non infected human macrophages previously reported by Fernandes et al. (26) . This analysis revealed a high number of differentially expressed genes (DEGs) in paired-sample analysis (non infected against infected macrophage) at the first time points after infection. At 4 hpi, we observed a higher number of DEGs, totaling 4,704. At this time, 2,518 were upregulated and 2,186 downregulated in the infected macrophages compared to control non infected macrophages. We observed that DEGs decreased as the time postinfection increased, with the lowest number of DEGs at 72 hpi, 950, of which 411 were downregulated and 539 were upregulated (Fig. 1A) .

Although the RNA-seq data correspond to a library obtained from poly A tailing RNAs, we identified the presence of a high number of differentially expressed ncRNAs. Most of these ncRNAs correspond to long noncoding RNAs (lncRNAs). We found at 4 hpi the highest number of overexpressed ncRNAs, with 428 upregulated (Fig. 1A) , of which 407 were annotated as lncRNAs, 8 microRNAs (miRNAs), 6 miscellaneous RNAs (miscRNAs), 5 small nucleolar RNAs (snoRNA), and 1 small nuclear RNA (snRNA) ( Table  S1 ). On the other hand, the lowest number of up- (35) and downregulated (112) ncRNAs were identified at 72 hpi ( Fig. 1A ; Table S1 ). Interestingly, we observed 15 upregulated lncRNAs that are shared in all four time points (Fig. 1C , upper right; Data set S1). On the other hand, 20 downregulated lncRNAs were shared in the four time points (Fig. 1C , upper left; Data set S1). However, more analysis must be performed to determine the biological role of these ncRNAs in the host's response to Leishmania infection.

As with lncRNAs, we evaluated the number of protein-coding genes differentially expressed, their relationship with the immune system, and the GO term enrichment at the four time points postinfection. We identified between 470 (72 hpi) to 1,890 (4 hpi) upregulated protein-coding genes (Fig. 1B) , of which 184 to 949 were related to immune system GO groups according to the ShinyGO analysis (27) (Data set S2), and 66 genes related to immune system were constantly upregulated in all four evaluated time points (Fig. 1C , bottom right; Data set S1). Here, we identify genes such as JUN or MYC, which have been reported as overexpressed in visceral leishmaniasis patients (28) . As we show in Fig. 1B , the number of downregulated protein-coding genes ranges from 266 (72 hpi) to 1,768 (4 hpi). Our analysis of functional groups related to the immune system reveals the presence of several proteins that range from 90 (72 hpi) to 643 (4 hpi). Interestingly, we identified only 17 downregulated genes in all four time points (Fig. 1C , left bottom; Data set S1).

We applied a GO enrichment analysis to identify biological processes associated with host-pathogen interaction, response to stress, and immune response. Enriched GO categories for biological processes obtained in upregulated DEGs from infected against non infected macrophages comparison are listed in Data set S3. Genes upregulated at 4 hpi were enriched in GO categories involved in cytokine responses, such as cellular response to cytokine stimulus (GO:0071345), cellular response to interleukin-1 (GO:0071347), an interleukin which promotes Th1 differentiation and inhibits disease progression in Leishmania major infections (29) , and cellular response to interleukin-7 (GO:0071347), which enhances the elimination of amastigotes (30) . On the other hand, terms related to the homeostasis of metals involved in macrophage function, such as zinc ion homeostasis (GO:0055069) (31), were also enriched. Mitogen-activated protein kinase (MAPK) cascade (GO:0000165), which is a determinant for IL-10 production and host susceptibility in Leishmania infections (32) , was enriched in our GO analysis. Moreover, other GO terms, such as regulation of I-kappaB kinase/NF-k B signaling (GO:0043122) and cellular response to tumor necrosis factor (GO:0071356), were enriched in upregulated genes in infected macrophages at 4 hpi (Fig. 1D ).

Similar to that in 4 hpi, in 24 hpi timepoint we identified that cellular response to cytokine stimulus (GO:0071345) and cytokine-mediated signaling pathway (GO:0019221) terms were enriched within upregulated DEGs at this time point (Data set S3; Fig. 1D ). Interestingly, we noticed the presence of enriched terms related to the regulation of macrophage differentiation at 24, 48, and 72 hpi ( Fig. 1D ; Data set S3). Antigen processing and presentation of exogenous peptide antigen via major histocompatibility complex (MHC) class I (GO:0042590) and Wnt signaling pathway, planar cell polarity pathway (GO:0060071), a pathway that canonically was related to defense against Leishmania infections (33) , were enriched at 48 hpi but not at other time points (Fig. 1D ). Also at 48 hpi, we identified terms such as cytokine-mediated signaling pathway (GO:0019221) and regulation of transcription from RNA polymerase II promoter in response to stress (GO:0043618) (Data set S3). Finally, our enrichment analysis showed that tumor necrosis factor-mediated signaling pathway (GO:0033209) and regulation of macrophage-derived foam cell differentiation (GO:0010743) were enriched at 72 hpi ( Fig. 1D ; Data set S3). Additionally, terms related to noncoding RNA (ncRNA) metabolism, such as ncRNA processing (GO:0034470), were consistently enriched in DEGs from the 24 and 48 hpi time points (Data set S3).

Context-specific gene regulatory networks of Leishmania-infected macrophages reveal potential new therapeutic targets. The TF-gene reference human GRN obtained after filtering high-quality connections is described in Table S2 . After that, we used the normalized expression values obtained from our previous analysis to filter this gold standard GRN. We obtained eight time-specific GRNs (four to infected and four to non infected macrophages, one for each time point). Each context-specific network presented a different number of nodes and connections, as described in Table S2 . The larger network for infected macrophages was obtained at 4 hpi, and it includes 19,750 nodes and 343,072 connections with 990 TFs. In addition, the smaller network of infected macrophages corresponds to 72 hpi and is made up of 19,718 nodes, of which 974 were identified as TFs, and 339,390 connections.

We compared each infected macrophage network against the control non infected macrophage networks obtained from the same time point. We found that non-TF genes show few alterations in the GRNs that represent macrophages in the first hours after infection, as indicated by F1 metric calculated by LoTo [90] for the presence or absence of network motifs in each compared GRN between 0.96 and 0.98. We also analyzed all TFs that had a significant change of regulation according to F1 values; for this, we used an F1 cutoff of 0.95 (Table 1 ). These results indicate that many TFs related to immune response show alterations in their regulations; in addition, not all these TFs were identified as differentially expressed according to our analysis. Therefore, we selected all nodes identified as differentially expressed and their connections present in the infected context networks and absent in non infected macrophage networks ( Fig Figure 2B shows the subnetwork after filtering all the edges present only in the infected macrophage network at 4 hpi. In this subnetwork, there are 160 connections and 63 nodes, of which 11 were TFs. Additionally, as listed in Table S3 , we found 244 nodes for 24 hpi network comparison (Fig. S1B) , 155 for 48 hpi (Fig. S2B) , and 52 for 72 hpi (Fig. S3B ).

To select a list of initial candidates, we evaluated all genes recovered from these comparisons to remove all non-DEGs, as well as genes not related to immune response, host-pathogen interaction, and/or stress. This first selection contains a total of 373 genes, of which only 113 passed all filters. Interestingly, we identified several TFs that act in the activation of the immune response, such as JUN, or the negative regulation of immune response, such as MYC. We also determined with this analysis that some effectors of the immune response, such as IL-16 (a proinflammatory interleukin), were present in our list of possible therapeutic targets. In Data set S4, we summarize all the identities, expression values, and functions of all 113 candidate genes.

Finally, we used the 113 selected genes as seeds to obtain the biological pathways in which their coded product participates. We obtained a list of 313 Reactome IDs of 0.936 Immune system process/leukocyte activation a TF, transcription factor; F1 represents the harmonic mean between precision and recall, ranging from 0 to 1, in which 1 represents a higher similarity of node X in both networks. b Downregulated gene. c Upregulated gene. hpi. Raw reads were analyzed using an in-house-developed pipeline that takes raw reads as input, and as output we obtained bona fide read counts per gene. Then, counts were used to obtain a normalized counts matrix and detect the differentially expressed genes. Next, we filtered a reference human GRN using normalized data to contextualize the GRN and get infected and non infected contexts simultaneously. After that, we applied a (Continued on next page) Reactome pathway entries linked to 909 UniProt IDs, which we downloaded and used to generate 313 mully multilayered graphs. After mapping the protein layers to our list of selected candidates, the mully graphs were filtered by deleting all nonprotein layers while adding the transitive edges to preserve the connections. Then, we deleted all protein nodes that were present in the mully graphs but not included in our list of candidates. Next, from these filtered graphs, we extracted the drug connections from UniProt IDs of our selected list and DrugBank, adding a drug layer to each filtered pathway graph. We combined all drug-protein connections from the different graphs to obtain the final list of drug targets (Fig. 3B) . A final set of 21 gene-pathway-drug interactions were obtained. In total, we identified 124 different biological pathways and 331 different drugs that have a direct connection with these 21 genes.

To reduce the number of drug-target interactions, we filtered this list using as selection criteria those drugs that were annotated as "approved," as well as "approved-investigational" or "approved-vet_approved" combinations. A total of 195 approved drugs were finally selected and presented direct interaction with 13 different genes. After this selection, we evaluated these 13 genes to identify their relationship with Leishmania infection. We found 8 genes that were confirmed in the literature as genes involved in Leishmania infection (34) (35) (36) (37) (38) . In Table S4 , we included the list of 8 genes that are potential candidates as host-directed therapeutic targets. We found that these 8 genes have a direct connection to 145 different drugs. This information was cross-linked with the expression data and other related metadata, such as drug-target interaction, evidence of the previous usage as antileishmanial drugs, current usage, or side effects, to determine the best set of drugs that could be repurposed as host-directed therapy for the treatment of leishmaniasis. After this analysis, we identified the coding product of five genes, androgen receptor (AR), C-Jun (JUN), platelet-derived growth factor receptor alpha (PDGFRA), prostaglandin endoperoxide synthase 2 (PTGS2), and vascular endothelial growth factor A (VEGFA), as potential therapeutic targets. In the pathway analysis, we identified that the product of these five genes has direct participation in pathways such as detoxification of reactive oxygen species, regulation of the apoptosome activity (JUN), or SCF/c-kit signaling pathway (PDGFRA) (Table  S4) . Furthermore, 11 drugs were selected by their direct interaction with the products of these five genes as well as their meta-features. Clascoterone and adapalene, which have previously been used as drugs to treat skin conditions such as acne, and other anti-inflammatories drugs, such as tolfenamic acid and flufenamic acid, were proposed as potential drugs for leishmaniasis chemotherapy. Moreover, the antipsychotic acetophenazine and the antineoplastic ripretinib were also selected. Table 2 displays more information about the potential targets, their interaction with drugs, their usage, and other relevant information. Our results indicated that several new therapeutic targets could be identified from the changes in gene regulation that occur during the infective process of Leishmania major of human macrophages, coupled with the integration of data related to pathways and drug-target direct connections.

Conventional therapies to tackle leishmaniasis are typically composed of only a few different drugs (39, 40) . Moreover, these therapies present several drawbacks, affecting mainly their efficiency, toxicity, and the ability to select resistant parasites (4) . Recently, several approaches for the treatment of parasitic diseases were focused on host-directed therapies (5), aiming to repurpose drugs previously used for other diseases (41) . This strategy search identifies new uses for drugs used or candidates in advanced clinical phases through network screening, followed by phase II and III pairwise comparison of infected against non infected contextualized networks to obtain the nodes and connections present in a disease condition. Next, we used the list of nodes to keep only genes involved in processes related to immune response, response to stress, or host-pathogen interaction and that were evidenced as differentially expressed. (B) Schematic workflow was applied to identify the drug targets using the Multipath package. With the filtered list, we mapped the gene set of interest to their gene products and related biological pathways in which these proteins participate and obtained the drug-gene product direct connection. Finally, drug-target interactions were literature filtered to select the best candidate targets for host-directed antileishmanial treatment. (42) . Here, we employed network-based approaches to identify novel therapeutic strategies by determining drug target proteins that play key roles in the host during leishmaniasis. Our results demonstrate how network-based approaches allow for the effective identification of new therapeutic targets for host-directed treatment in parasitic diseases. Compared to the traditional process, drug repositioning has some advantages, in that its cost and development time are reduced. In this study, we found five genes whose encoded products have the potential to be new therapeutic targets. We identified a direct connection between these gene products and 11 FDA-approved drugs (Table 2, Fig. 4 ). Many of these drugs have previously been used to treat skin disease (e.g., acne) as anti-inflammatories or antineoplastics.

Drugs such as clascoterone are antagonists of the androgen receptor (AR) for the treatment of androgen-dependent skin diseases, including androgenetic alopecia and acne. Previous reports have shown that AR has a relevant role in the immune response derived from parasitic infections (38) . Two different works, reported by Sánchez-García and colleagues and Qiao and colleagues, discovered that hormones such as dihydrotestosterone (DHT) and testosterone that interact with the androgen receptor can alter parasites' development and survival or inhibit apoptosis in Leishmania-infected macrophages (43, 44) . These effects seem to be exerted by specific receptors for androgen in macrophages, leading to greater replication of the parasite and an increased rate of infection, with increased proliferation of the parasite. Acetophenazine is a moderately potent antipsychotic and antagonist of the dopamine D2 receptor and the androgen receptor. Previous studies show that dopaminergic receptor antagonists can inhibit the growth and multiplication of toxoplasmosis parasites (45) .

Esculin, a phenolic compound, typically used as a vasoprotective agent, was tested to evaluate their inhibitory activity against Leishmania infantum arginase (ARGLi), demonstrating a null inhibitory activity in in vitro assay (46) . On the other hand, this compound promotes a reduction of transmission and viability of L. infantum and Leishmania mexicana in a vector model assay (47) . Notwithstanding that our findings showed AR as a possible target for leishmaniasis treatment, studies with tamoxifen, an anticancer drug that acts as a nonsteroidal estrogen receptor modulator, presents antileishmanial activity (48, 49) by inhibiting the parasite's inositol phosphorylceramide (50) with reduced chances of selecting resistant parasites (51) . It is worth noting that tamoxifen treatment resulted in scrotal swelling that leads to infertility in cutaneous experimental leishmaniasis caused by L. major in male mice (52) . Androgen is converted to estrogen by an aromatase enzyme that can be inhibited by testosterone (53) . In humans, combined therapy using oral tamoxifen and meglumine antimoniate in patients suffering from cutaneous leishmaniasis resulted in cure rates similar to those of conventional schemes (54) . Initially used to treat breast cancer, tamoxifen is a classic example of a repurposed drug with different targets in the parasite that was promising during in vitro investigations but presented limitations when following the drug development pipeline.

Regarding the viability of other potential targets in our list, we found that prostaglandin endoperoxide synthase 2 (PTGS2) was previously reported as a biomarker for the response to infections for L. major (34) . Our results showed a consistent misregulation of this gene that is consistent with these previous studies. Also, we identified a great number of direct connections between PTGS2 and different drugs, such as tolfenamic acid, a nonsteroidal anti-inflammatory drug (NSAID) usually used for migraine pain. Also, this drug shows antibacterial activity against Staphylococcus aureus (55), Burkholderia pseudomallei (56), Liberibacter asiaticus (57) , and in lesser potency, Vibrio cholerae (58) . Additionally, in a recent study, it is shown to be effective to kill schistosomes in in vitro and in vivo assays (59) . In this regard, Leishmania braziliensis prostaglandin F2a synthase (LbrPFG2S) was associated with host-parasite interaction playing a crucial role in pathogenicity through proinflammatory lipid synthesis (60, 61), revealing and validating a potential intervention in a parasite's prostaglandin biosynthesis pathway.

Our approach also found different anticancer drugs among those targeting our list of candidates. This type of drug has been used before to treat leishmaniasis. For example, miltefosine was developed as an anticancer drug, but today it is a choice for visceral leishmaniasis treatment (62, 63) . Ripretinib was described as a promising drug for the treatment of gastrointestinal stromal tumors, and more recently its drug family has been proposed to be repositioned for COVID-19 treatment (64). Our approach also indicated that ripretinib might be useful for the treatment of leishmaniasis. This is because ripretinib acts as a kinase inhibitor and presents an inhibition activity over that of platelet-derived growth factor receptor (PDGFRA), a tyrosine kinase receptor previously targeted for antileishmanial therapies that presents a significant reduction in parasitic survival (65) . Indeed, the anticancer drug sunitinib, a broad-spectrum receptor tyrosine kinase inhibitor, blocked progressive splenomegaly and improved immunity as adjuvant therapy in murine experimental leishmaniasis (66) . Flufenamic acid, a nonsteroidal anti-inflammatory, showed good antimalarial activity, significantly retarding the intraerythrocytic growth of Plasmodium falciparum, but its antileishmania potential still needs to be tested (67) .

Another candidate target found by our method is vascular endothelial growth factor A (VEGFA), a gene whose expression is induced during Leishmania major infection (68, 69) . The product of this gene promotes lymphangiogenesis, a process that is involved in the inflammatory response and lesion healing (68) . The product encoded by VEGFA is used as a target in antiangiogenesis therapies to reduce the vascularization of tumors (70) . Minocycline is an antibiotic from the tetracycline family that interacts with the VEGFA protein and inhibits angiogenesis (71) . This antibiotic has been proposed to treat parasitic infection in the late 1980s, mainly against Giardia lamblia, due to a potent activity that reduced the survival of the parasite in vitro (72) . The use of minocycline in the treatment of malaria has shown a decrease in T-cell-mediated brain inflammation and a reduction in gene expression independent of the antiparasitic property (73) . Nonetheless, a more recent study reported that a patient diagnosed with leishmaniasis was treated with minocycline (74) . However, the treatment was stopped shortly after and nonconclusive results were obtained. Pidolic acid is an active form of 5-oxoproline and occurs in fundamental biological processes such as intracellular stage differentiation, host cell infection, and resistance to various stresses in protozoa belonging to the genera Trypanosoma and Leishmania. Fargnoli et al. showed in an unprecedented way that L-proline uptake has been proposed as a chemotherapeutic target for Trypanosoma cruzi, opening a new horizon in the development of new chemotherapeutics against Chagas disease and other parasitic diseases (75) .

The last target in our list is JUN, a gene that encodes the c-Jun protein that is a central part of the AP-1 transcription factor. This transcription factor is crucial for the inactivation of macrophages during Leishmania infection (76) . Adapalene is a retinoid that acts as a comedolytic and anti-inflammatory agent that has been widely used due to its milder effects compared to those of other substances of the same class, mainly used for acne treatment that targets c-JUN via the AP-1 transcription factor (77, 78) . Interestingly, adapalene has been used as an anticancer drug in an in vitro assay promoting the apoptosis of colorectal cells (79) . The antileishmanial activity of adapalene was tested in vitro and in vivo in animal models, showing high antileishmanial activity and promoting healing in hamsters infected with Leishmania panamensis (80). This drug not only affects the parasite but also promotes the host's immune response, showing that it could be a very effective drug for leishmaniasis treatment.

We have created a new strategy for the search for therapeutic targets. We applied our approach to finding host-directed antileishmanial therapeutic targets and showed that it provides a significant number of novel potential targets. This new strategy helps to bypass common issues arising from conventional antiparasitic therapies, such as the fast appearance of resistance and strong side effects that preclude generalized drug usage. Furthermore, we demonstrate in this study that many drugs could be repositioned for leishmaniasis treatment. Importantly, all drugs selected in our work need to be experimentally tested to confirm their potential as host-directed antileishmanial therapies. Notably, after demonstrating the potential of our new strategy for the identification of therapeutic targets, it is important to highlight that it requires only transcriptomic data that is integrated with other available data, making it an easy tool to adapt for other diseases.

Our results indicate the possibility of repurposing several drugs that could be useful as antileishmanial therapies. However, even with strong in silico evidence supporting our list of new therapeutic targets, our approach still requires in vivo testing to confirm that these repurposed drugs would be useful as treatments for Leishmania infections.

RNA-seq data sets from Leishmania-infected macrophages. A set of 43 RNA sequencing (RNAseq) samples derived from human macrophages infected with Leishmania major (BioProject accession number PRJNA290995) were downloaded from NCBI SRA (81) . This data set is composed of 4 time points, 4 h postinfection (hpi), 24 hpi, 48 hpi, and 72 hpi, five biological replicates per time point, and six biological replicates for uninfected macrophages as a control at the same time points as infected macrophages.

RNA-seq data analysis. RNA-seq raw data quality control inspection was performed using FastQC version 0.11.8 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Low-quality reads were removed with Trimmomatic version 0.36 (82) using a Phred cutoff value (Q) of 30. Then, we mapped the remaining high-quality reads with Hisat2 version 2.1.0 (83) to the human genome (GRCh38), downloaded from Ensembl (84) . Expression values were calculated using HTSeq-count version 0.7.2 (85) . The resulting counting reads matrix was normalized and used to identify the differentially expressed genes through DESeq2 version 1.24.0 (86). Genes with an adjusted P value of #0.5 and absolute log fold change of $0.5 were considered differentially expressed.

Gene regulatory network analyses. A human reference GRN was obtained from DoRothEA release 1.3.3 (87) , considering only high-quality connections (A, B, and C evidence codes reported in their original work). This reference GRN was then filtered using RNA-seq normalized data in the same way as that described by Santander et al. (88) , with some modifications related to the expression threshold employed to contextualize the GRN. Here, interactions between a transcription factor (TF) and its target gene were saved only if DESeq2 normalized expression values (median ratios) for that TF were at least 10 (89) . Then, we obtained contextualized networks for each condition (8 GRNs, 4, 24, 48 , and 72 hpi, infected and non infected human macrophages) and applied a pairwise comparison of infected against non infected for all four serial times using LoTo (90) . LoTo identifies genes whose regulatory environment varies based on binary classification metrics calculated with the presence or absence of network motifs in each compared GRN, determining differentially regulated nodes. Once these comparisons were obtained, we evaluated the F1 values of all TF and non-TF genes. F1 represents the harmonic mean between precision and recall, ranging from 0 to 1, in which 1 represents a higher similarity of node X in both networks (90) . Next, we crossed the list of differentially regulated genes in each comparison with the list of differentially expressed genes from DESeq2 analysis. Finally, gene candidates with F1 lower than 0.95 and an absolute value of logFC greater than 0.5 were considered in further analyses. The resulting list was functionally enriched using EnrichR version 3.0 (91-93) and ShinyGO v0.61 (27) and then filtered to keep only genes associated with processes related to immune response, response to stress, or host-pathogen interaction.

Gene-protein, protein-pathway, and drug-target interaction mapping. Multipath version 1.0.3 is an R package used to generate integrated reproducible pathway knowledge (94) . Using Multipath, BioPAX-encoded pathways (95) can be modeled into multilayered graphs, where the biological pathways components are embedded into different layers based on their biological type. The built graphs are reproducible, i.e., all modifications applied to the graphs are stored. Multipath is also used to integrate influencing pathway knowledge from external databases like drugs from DrugBank version 5.1.8 (96) . We used this package to query pathway knowledge databases and fetch relevant information needed in our computational analysis. To map the gene set of interest to their gene products from UniProt, Multipath uses UniProt.ws to fetch the UniProt IDs of the corresponding proteins, which were mapped to the list of candidates. Then, we got a list of biological pathways from Reactome version 73 Released (97) in which these proteins participate. These Reactome IDs were downloaded to generate mully multilayered graphs (98) . Next, we filtered all proteins that were not coded by genes in our previous list. Finally, we extracted the drug targets from UniProt release 2020_05 and DrugBank version 5.1.8 and added a drug layer to each filtered pathway graph, preserving only those genes for which drug direct connection was identified (Fig. 3B) .

Selection of potential host-directed therapeutic targets for leishmaniasis treatment. After drugtarget direct interactions were obtained, we filtered this network to reduce the number of potential targets. To do so, we first selected those drug-target connections in which the drug had at least one of the following labels: approved, approved and investigational, or approved and vet_approved. After that, literature mining was applied to identify all genes that were playing a role in Leishmania infection and were identified as differentially expressed in infected macrophages in the transcriptomic analysis. Next, we used this information to refine the filtering of selected drug-target interactions. Finally, we evaluated all filtered drugs to obtain the type of drug-target interaction, evidence of previous usage as antileishmanial drugs, actual usage, side effects, pharmacological action, and any other relevant information for the best target selection. Figure 3 summarizes the entry pipeline employed in this work.

Supplemental material is available online only. fellow (312965/2020-6). The funders had no role in the design of the study, collection, analysis, and interpretation of data, or writing the manuscript.

Powered@NLHPC: this research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02); and by the computing infrastructure of the Centro de Genómica y Bioinformática, Universidad Mayor.

We declare no conflicts of interest.

Advances in leishmaniasis

Recent advances and new strategies on leishmaniasis treatment

A review of current treatments strategies based on paromomycin for leishmaniasis

Host-directed drug therapies for neglected tropical diseases caused by protozoan parasites

Hostdirected therapies for bacterial and viral infections

Host-directed therapies for infectious diseases: current status, recent progress, and future prospects

Potential of immunomodulatory agents as adjunct host-directed therapies for multidrug-resistant tuberculosis

Red cells from ferrochelatase-deficient erythropoietic protoporphyria patients are resistant to growth of malarial parasites

Coaxing HIV-1 from resting CD4 T cells: histone deacetylase inhibition allows latent viral expression

Histone deacetylase inhibitors for purging HIV-1 from the latent reservoir

Reducing mortality from 2019-nCoV: host-directed therapies should be an option

Identification of SARS-CoV-2-induced pathways reveals drug repurposing strategies

Immune checkpoint targets for host-directed therapy to prevent and treat leishmaniasis

Interleukin-12 regulates the response to chemotherapy in experimental visceral leishmaniasis

The Abl and Arg kinases mediate distinct modes of phagocytosis and are required for maximal leishmania infection

Ibrutinib is an irreversible molecular inhibitor of ITK driving a Th1-selective pressure in T lymphocytes

Berberine chloride mediates its anti-leishmanial activity via differential regulation of the mitogen activated protein kinase pathway in macrophages

The leishmanicidal activity of oleuropein is selectively regulated through inflammation-and oxidative stress-related genes

A paradigm shift in medicine: a comprehensive review of network-based approaches

A guide to conquer the biological network era using graph theory

Gene-centered regulatory network mapping

Network medicine in the age of biomedical big data

Multilayer networks: aspects, implementations, and application in biomedicine

Multilayer and multiplex networks: an introduction to their use in veterinary epidemiology

Dual transcriptome profiling of leishmania-infected human macrophages reveals distinct reprogramming signatures

ShinyGO: a graphical gene-set enrichment tool for animals and plants

Expression profiling of Sudanese visceral leishmaniasis patients pre-and posttreatment with sodium stibogluconate

Interleukin 1a promotes Th1 differentiation and inhibits disease progression in leishmania majorsusceptible BALB/c mice

Interleukin-7 enhances antimicrobial activity against Leishmania major in murine macrophages

The role of zinc and zinc homeostasis in macrophage function

Activation of the MAPK, ERK, following Leishmania amazonensis infection of macrophages

Wnt5a signaling promotes host defense against Leishmania donovani infection

In search of biomarkers for pathogenesis and control of leishmaniasis by global analyses of -infected macrophages

Growth factor and Th2 cytokine signaling pathways converge at STAT6 to promote arginase expression in progressive experimental visceral leishmaniasis

Profiling gene expression of antimony response genes in Leishmania (Viannia) panamensis and infected macrophages and its relationship with drug susceptibility

Immunotherapeutic potential of eugenol emulsion in experimental visceral leishmaniasis

Androgen-dependent immune modulation in parasitic infection

Emerging therapeutic targets for treatment of leishmaniasis

Current and emerging medications for the treatment of leishmaniasis

Drug repurposing and human parasitic protozoan diseases

Drug repurposing: progress, challenges and recommendations

Dihydrotestosterone enhances growth and infectivity of Leishmania mexicana

Testosterone inhibits apoptosis of Leishmania donovani-infected macrophages

Contribution of dopamine neurotransmission in proconvulsant effect of Toxoplasma gondii infection in male mice

Leishmania infantum arginase: biochemical characterization and inhibition by naturally occurring phenolic substances

Antileishmanial activity of furoquinolines and coumarins from Helietta apiculata

Tamoxifen is effective against Leishmania and induces a rapid alkalinization of parasitophorous vacuoles harbouring Leishmania (Leishmania) amazonensis amastigotes

Tamoxifen induces apoptosis of Leishmania major promastigotes in vitro

Tamoxifen inhibits the biosynthesis of inositolphosphorylceramide in Leishmania

Leishmania is not prone to develop resistance to tamoxifen

Leishmania major: activity of tamoxifen against experimental cutaneous leishmaniasis

Aromatase and regulation of the estrogen-to-androgen ratio in synovial tissue inflammation: common pathway in both sexes

Tamoxifen and meglumine antimoniate combined therapy in cutaneous leishmaniasis patients: a randomised trial

Search of antimicrobial activity of selected non-antibiotic drugs

Immune modulation as an effective adjunct post-exposure therapeutic for B. pseudomallei

Identification of the tolfenamic acid binding pocket in PrbP from

Cellular mechanisms underlying the inhibitory effect of flufenamic acid on chloride secretion in human intestinal epithelial cells

Phenotypic screening of nonsteroidal anti-inflammatory drugs identified mefenamic acid as a drug for the treatment of schistosomiasis

Leishmania braziliensis prostaglandin F synthase impacts host infection. Parasit Vectors

Leishmania infantum lipophosphoglycan induced-prostaglandin E production in association with PPAR-g expression via activation of Toll like receptors-1 and 2

Leishmaniasis treatment: update of possibilities for drug repurposing

Miltefosine: a review of its pharmacology and therapeutic efficacy in the treatment of leishmaniasis

Diarylureas: repositioning from antitumor to antimicrobials or multitarget agents against new pandemics. Antibiotics (Basel) 10

Activity of anti-cancer protein kinase inhibitors against Leishmania spp

Inhibition of receptor tyrosine kinases restores immunocompetence and improves immune-dependent chemotherapy against experimental leishmaniasis in mice

Inhibition of eryptosis and intraerythrocytic growth of Plasmodium falciparum by flufenamic acid

Infection induces macrophage vascular endothelial growth factor A production in an ARNT/ HIF-dependent manner

Leishmania major infection-induced VEGF-A/VEGFR-2 signaling promotes lymphangiogenesis that controls disease

VEGF as a therapeutic target in cancer

Minocycline inhibits angiogenesis in vitro through the translational suppression of HIF-1a

Tetracyclines as antiparasitic agents: lipophilic derivatives are highly active against Giardia lamblia in vitro

Minocycline prevents cerebral malaria, confers neuroprotection and increases survivability of mice during Plasmodium berghei ANKA infection

Diagnosis of Leishmania infection in a nonendemic area of South America

Targeting L-proline uptake as new strategy for anti-chagas drug development

Leishmania-induced inactivation of the macrophage transcription factor AP-1 is mediated by the parasite metalloprotease GP63

Application of ultrasound-mediated adapalene-coated lysozyme-shelled microbubbles in UVA-induced skin photoaging

The integral role of topical and oral retinoids in the early treatment of acne

The synthetic retinoid adapalene inhibits proliferation and induces apoptosis in colorectal cancer cells in vitro

In silico screening of potential drug with antileishmanial activity and validation of their activity by in vitro and in vivo studies

International Nucleotide Sequence Database Collaboration

Trimmomatic: a flexible trimmer for Illumina sequence data

HISAT: a fast spliced aligner with low memory requirements

HTSeq-a Python framework to work with high-throughput sequencing data

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Benchmark and integration of resources for the estimation of human transcription factor activities

Transcriptional profiling of embryos lacking the lipoprotein receptor SR-B1 reveals a regulatory circuit governing a neurodevelopmental or metabolic decision during neural tube closure

Differential expression analysis for sequence count data

LoTo: a graphlet based method for the comparison of local topology between gene regulatory networks

Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool

Enrichr: a comprehensive gene set enrichment analysis web server 2016 update

Gene set knowledge discovery with Enrichr

Multipath: an R package to generate integrated reproducible pathway models

The BioPAX community standard for pathway data sharing

DrugBank 5.0: a major update to the DrugBank database

The reactome pathway knowledgebase

mully: an R package to create, modify and visualize multilayered graphs

Antimicrobial peptide-induced apoptotic death of leishmania results from calcium-dependent, caspase-independent mitochondrial toxicity

Antimonial-induced increase in intracellular Ca21 through non-selective cation channels in the host and the parasite is responsible for apoptosis of intracellular Leishmania donovani amastigotes