Dissection of the macrophage response towards infection by the Leishmania-viral endosymbiont duo and dynamics of the type I interferon response

Introduction

Macrophages, one of the first-line defenders of our immune system, are implicated in many processes ranging from physiological to pathological. One of such processes includes the recognition and ultimately elimination of pathogens. In order to achieve that, macrophages employ a plethora of mechanisms such as oxidative burst, inflammation, or antigen presentation (Shapouri-Moghaddam et al., 2018; Rossi and Fasel, 2018a; Ferrari, 2019; Sheu and Hoffmann, 2022; Reverte et al., 2022). However, these responses are not always straightforward nor unidirectional as pathogens use different strategies to hijack the macrophage machinery and ultimately exploit it to their advantage. This is particularly true for intracellular pathogens such as Leishmania that have evolved over thousands of years to counteract the microbicidal power of macrophages (Gupta et al., 2013; Rossi and Fasel, 2018a; Reverte et al., 2022).

Leishmania are medically relevant protozoan parasites that cause leishmaniasis, a neglected tropical disease that can manifest depending on the infecting species and the immune status of the host. These manifestations can range from skin ulcers to the disfiguration of the nasopharyngeal area to even potentially lethal forms if left untreated (Pearson and Anastacio de Queiroz, 1996; Scorza et al., 2017; Burza et al., 2018). Although faced with a hostile environment when phagocytosed by macrophages, Leishmania have evolved mechanisms to subvert the macrophage response and ultimately prevail. Leishmania can, for example, induce the SOCS (suppressor of cytokine signaling) proteins that are negatively regulating the JAK/STAT pathway and the cytokine release leading to the impairment of the crosstalk of the macrophage with the T cell, resulting in the decrease of interleukin 12 (IL-12) and interferon (IFN) gamma (IFN-γ) productions (Chandrakar et al., 2020). Moreover, it has been shown that Leishmania are able to block AIF1 (allograft inflammatory factor-1) expression in macrophages and degrade CXCL1 to inhibit the pro-inflammatory responses and the recruitment of neutrophils, respectively (Yorek et al., 2019; da Silva et al., 2021). Parasites can also inhibit the inflammasome through the upregulation of the A20 protein (Hartley et al., 2018; Zamboni and Sacks, 2019) or interfere with the oxidative stress response (Reverte et al., 2021; Reverte et al., 2022).

Leishmania RNA virus 1 (LRV1) is a double-stranded RNA (dsRNA) virus of the Totiviridae family living endosymbiotically in some strains of Leishmania (Tarr et al., 1988; Stuart et al., 1992; Zamora et al., 2000; Hartley et al., 2012). The presence of LRV1 could contribute to the exacerbation and chronicity of leishmaniasis, the development of destructive metastatic lesions, and even treatment failure (Ives et al., 2011; Cantanhede et al., 2015; Adaui et al., 2016; Bourreau et al., 2016; Hartley et al., 2016). The presence of LRV1 inside Leishmania results in increased macrophage survival in a mechanism dependent on Akt (Eren et al., 2016); an increased production of pro-inflammatory cytokines and chemokines such as interleukin 6 (IL-6), tumor-necrosis factor α (TNF-α), interleukin 17 (IL-17), type I IFNs (IFN-α and IFN-β) (Ives et al., 2011; Hartley et al., 2016; Rossi et al., 2017); and the inhibition of inflammasomes via NLRP3 (Hartley et al., 2018). This immunophenotype depends on the recognition of LRV1 by the endosomal Toll-like receptor 3 (TLR3) (Ives et al., 2011; Olivier and Zamboni, 2020) and the activation of its downstream type I IFN signaling pathways (Rossi et al., 2017). Interestingly, despite the inflammatory environment triggered by LRV1, mice are more susceptible to Leishmania guyanensis (Lgy) bearing LRV1 (from now on LgyLRV1+) than to Lgy devoid of the virus (i.e., LgyLRV1-). Interestingly, the same exacerbatory effect is observed when an exogenous virus is coadministered with LgyLRV1-, resulting in an increase in parasite burden and footpad swelling (Rossi et al., 2017) underlining the importance of the systemic production of type I IFN and its signaling pathways for the exacerbation of the affliction outcome. Mice lacking Ifnar (Ifnar^-/-), the gene coding for IFNAR, the receptor of type I IFNs, exhibited a similar lesion size and parasitemia when infected with both LRV1+ and LRV1- parasites (Rossi et al., 2017). Additionally, infection with a virus even after the resolution of the Leishmania infection has the potential to reactivate the disease caused by Leishmania (Rossi et al., 2017), raising concerns on leishmaniasis control and management. Thus, viral coinfections or later exposure to any potential trigger of the type I IFN response can be considered a risk factor for leishmaniasis relapses. Despite the relevance of viral infection in the context of leishmaniasis, not many studies exist on the role of virus coinfection with Leishmania, with the exception perhaps of a coinfection with HIV, known to modulate and impair the adaptive immune response and can thus be considered a predictor of a worsened leishmaniasis outcome (Lindoso et al., 2016).

The advent of transcriptomics has contributed to the global comprehension of how macrophages respond toward infection, including infection by Leishmania. A large number of genes modulated by Leishmania were related to the immune response (pro- and anti-inflammatory), glycolysis, lipid metabolism, biogenesis, and phagocytosis (Fernandes et al., 2016; Aoki et al., 2017; Aoki et al., 2019; Shadab et al., 2019; Restrepo et al., 2021; Reverte et al., 2021; Salloum et al., 2021; Chaparro et al., 2022). While these studies have shed light on a plethora of mechanisms potentially perturbed by Leishmania, they have not explored the impact that coexposure to additional agents may have on the immune response mounted toward Leishmania. As mentioned above, this point is particularly relevant when macrophages are concomitantly exposed to the Leishmania and its endosymbiont, the LRV1 duo. Understanding how the host responds toward infection by Leishmania and the contribution of the presence of the virus and its downstream type I IFN response is paramount to design strategies for improving disease prevention, progression, and the outcome. To have a comprehensive and dynamic understanding of the responses mounted by the macrophage toward the Leishmania-LRV1 duo as well as to dissect their individual impact, their synergies or their opposing effects, an RNA sequencing (RNA-Seq) analysis was performed. Here, the different cellular pathways affected by the presence of Leishmania or a virus were disentangled and the genes driving those pathways were identified. In addition, the pathways exacerbated by the presence of LRV1 inside Leishmania that contributed to worsening the disease fallout were spotted. Interestingly, although non-pathogenic for humans, LRV1 had the potential to surpass the effect of Leishmania and ultimately drive the immune response of the macrophage. The dynamic and timely nature of the orchestrated macrophage response was demonstrated by the importance of the genes that had an early, late, or sustained response during infection by the Leishmania-LRV1 duo.

Materials and methods

Animals

Wild-type (WT; C57BL/6JOlaHsd) mice were purchased from Envigo (Gannat, France), Ifnar^-/- (B6.129S2 Ifnar1^tm1Agt/Mmjax) mice were obtained from M. Aguet, Swiss Institute of Experimental Cancer Research, Epalinges, Switzerland. Mice were maintained and bred in GM500 IVC Green Line cages (Tecniplast, Buguggiate, Italy) at the specific pathogen-free (SPF) animal facility at the University of Lausanne (Epalinges site), with a relative humidity 55 ± 10%, 21 ± 2°C, and in 11/13-h dark/light cycles with ad libitum water (local acidified and autoclaved water) and food (Kliba Nafag or Safe). Cardboard or plastic tunnels and igloos, as well as paper tissues, were supplied as enrichment. Bone marrow was extracted from 7-week-old female WT mice (RNA-Seq dataset#1), 7- to 10-week-old female WT mice (RNA-Seq dataset#2), and 9-week- old female Ifnar^-/- mice (RNA-Seq dataset#2) (Table 1). Experiments were performed according to the ethical guidelines set out by the Swiss Federal Food Safety and Veterinary Office (FSVO), and procedures were approved by the Veterinary Commission of the Canton de Vaud (SCAV, Switzerland) under the authorization numbers 2113.3.

TABLE 1

Table 1 Summary table of RNA sequencing samples used for the analysis.

Parasites

The study was performed using lines of Leishmania guyanensis (MHOM/BR/78/M4147), named LgyLRV1+ and LgyLRV1-, bearing LRV1 or not, respectively (Kuhlmann et al., 2017). Both parasites were recovered from footpads of infected C57BL/6 mice and not kept for longer than six passages. Leishmania strains were cultured in a complete Schneider’s medium, containing a Schneider’s Drosophila Medium (Gibco, Thermo Fischer Reinach, Switzerland) supplemented with 20% heat-inactivated fetal bovine serum (FBS; Gibco), 1% penicillin–streptomycin solution (BioConcept, Allschwill, Switzerland) (i.e., 100 IU/ml and 100 µg/ml, respectively), 2% 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES) (BioConcept) (i.e., 20 mM), hemin–folate solution (prepared from Porcine Hemin, Sigma-Aldrich and Folic Acid, Fluka) (5 and 10 µg/ml, respectively), and 6-Biopterin (Sigma-Aldrich, Buchs, Switzerland) (0.6 µg/ml). They were maintained at 26°C and diluted every week. Six-day-cultured (stationary phase) parasites were used for macrophage infection.

Bone marrow–derived macrophages

Macrophages were derived from the bone marrow of C57BL/6 mice (WT or Ifnar^-/-). The mice of each genotype were sacrificed, and hind leg bones were collected (both tibia and femur) and kept in the Dulbecco’s modified Eagle’s medium (DMEM; Gibco) on ice. The same following procedure was applied to each mouse. After cleaning the bones with 70% ethanol and the DMEM, the bone marrow was flushed out with the complete DMEM, containing DMEM, supplemented with 10% heat-inactivated FBS (Gibco), 1% penicillin–streptomycin solution (BioConcept) (i.e., 100 UI/ml and 100 µg/ml, respectively), and 1% HEPES (BioConcept) (i.e., 10 mM). The bone marrow suspension was then passed through a 40-µm cell strainer (Corning) and centrifuged for 10 min, at 1,500 rpm (453 × g), at 4°C. The cell pellet was resuspended with the complete DMEM medium containing 50-ng/ml recombinant mouse macrophage colony stimulating factor (M-CSF, ImmunoTools). Approximately 10 ml were distributed per Sterilin Petri Dish (Thermo Scientific), and six plates were prepared. After 3 days in culture at 37°C and 5% CO₂, 5 ml of the complete DMEM medium containing 50-ng/ml M-CSF were added per plate. The cells were differentiated into macrophages for three additional days, that is, 6 days total, at 37°C and 5% CO₂. After these 6 days, the bone marrow–derived macrophages were detached with 1X Dulbecco's phosphate-buffered saline (DPBS) (Gibco) containing 5 mM ethylenediamine tetraacetic acid (EDTA), washed and resuspended with the complete DMEM medium. The cells were then counted and adjusted to 2.6 million cells per milliliter. They were respectively plated in 12-well plates (TPP) with 1.95 million cells per well. They were finally put back in culture at 37°C and 5% CO₂ for 1 day.

Macrophage infection and treatments

The day after their plating, stationary-phase Leishmania promastigotes were centrifuged, washed with 1X DPBS, resuspended with a complete DMEM medium, counted, and adjusted to the concentration required. Murine bone marrow–derived macrophages (BMDMs) were infected with the Leishmania at multiplicity of infection (MOI) of three and five parasites per macrophage, in the first and the second RNA-Seq datasets, respectively. Cells were also alternatively treated with some synthetic TLR agonist or cotreated with recombinant murine type I IFNs. Poly I:C (HMW), a TLR3 agonist (mimicking dsRNA), was obtained from In vivoGen and used at 2 µg/ml. Murine recombinant IFN-α and IFN-β were used each at 1,000 IU/ml (CellScience). LgyLRV1-infected BMDMs were treated or not with type I IFNs at 6 h postinfection (p.i.). BMDMs in the complete DMEM alone were used as a control. Cells were incubated at 35°C and 5% CO₂ until 8 and 24 h p.i. or posttreatment (Table 1).

RNA isolation and selection

At 8 and 24 h posttreatment, BMDMs cultured in 12-well plates were cleared out of supernatants and lysed with 350 µl of the RLT buffer (Qiagen) supplemented with 40-mM DTT (Dithiothreitol) for RNA extraction. Plates were then frozen at -80°C until the RNA purification. The RNA samples were purified with the RNeasy Plus Mini Kit (Qiagen) following the manufacturer’s instructions. Purified RNA was eluted with 30 µl of RNase-free water (Qiagen), and their RNA concentrations were measured with a spectrophotometer (NanoDrop, ThermoFisher Scientific). The quality of the RNA samples was then evaluated with a Fragment Analyzer (Agilent) at the Lausanne Genomic Technologies Facility (GTF) for both RNA-Seq datasets. For the dataset#2, the RNA quantification was performed with a fluorimetric method (Ribogreen; ThermoFisher Scientific) by the GTF according to the manufacturer recommendations, then used to prepare the RNA dilutions for the libraries. While for the dataset#1, the total integrated concentration (TIC) values given by the Fragment Analyzer were directly used to prepare the RNA dilutions for the libraries. The RNA samples diluted for the libraries were the ones that have an RNA quality number (RQN) higher than 6 (for dataset#2) and 7 (for dataset#1) and from 3 mice out of 4 and 5, respectively. The samples were diluted in 96-well plates (Labgene Scientific). Then, from these plates, the libraries and the sequencing were performed at the GTF.

RNA sequencing data processing

Purity-filtered reads were adapters and quality trimmed with Cutadapt [v. 1.8 (Martin, 2011)]. Reads matching to ribosomal RNA sequences were removed with fastq_screen (v. 0.9.3 for dataset#1 and v. 0.11.1 for dataset#2). The remaining reads were further filtered for low complexity with reaper [v. 15-065, (Davis et al., 2013)]. Reads were aligned against the Mus musculus.GRCm38.86 genome using STAR [v. 2.5.2b, (Dobin et al., 2013)] for dataset#1 and the Mus musculus.GRCm38.92 genome using STAR [v. 2.5.3a, (Dobin et al., 2013)] for dataset#2. The number of read counts per gene locus was summarized with htseq-count [v. 0.6.1 for dataset#1 and v. 0.9.1 for dataset#2, (Anders et al., 2015)] using Mus musculus.GRCm38.86 gene annotation for dataset#1 and Mus musculus.GRCm38.92 gene annotation for dataset#2. The quality of the RNA-seq data alignment was assessed using RSeQC [v. 2.3.7, (Wang et al., 2012)]. Reads were also aligned to the Mus musculus.GRCm38.86 transcriptome using STAR [v. 2.5.2b, (Dobin et al., 2013)] for dataset#1 and to the Mus musculus.GRCm38.92 transcriptome using STAR [v. 2.5.3a, (Dobin et al., 2013)] for dataset#2. The estimation of the isoform abundance was computed using RSEM [v. 1.2.31, (Li and Dewey, 2011)].

Normalization and data transformation

Statistical analysis was performed for genes independently in the software environment R. Genes with a low number of counts were filtered out according to the rule of one count(s) per million (cpm) in at least one sample. Library sizes were scaled using TMM normalization. Subsequently, the normalized counts were transformed to cpm values, and a log2 transformation was applied, by means of the function cpm with the parameter setting prior.counts = 1 [EdgeR (Robinson et al., 2010)]. For the analysis with Ifnar^-/- and WT genotypes, two sets of data generated separately were merged and this to integrated WT samples treated with IFN-α and IFN-β. To check if this was possible, we plotted the distributions of the log cpm of the shared conditions between the two datasets (LgyLRV1+ and LgyLRV1-). The distributions had similar patterns at 8- and 24-h time points, and the different groups of conditions were well separated (Supplementary Figure S1). Given the fact that with WGCNA, we are looking at expression patterns across conditions and not individual gene differential expression, we found this sufficient to take advantage of having more experimental conditions by integrating a dataset with WT treated with IFN-α and IFN-β.

Weighted gene coexpression network analysis and downstream bioinformatics analysis

Weighted gene coexpression network analysis (WGCNA) (Langfelder and Horvath, 2008) was performed on normalized RNA-Seq data in R (package WGCNA 1.69). For each analysis, an adjacency matrix was calculated to construct a signed hybrid coexpression network using Spearman correlation. A sequence of soft-thresholding powers was tested to reach a free-topology network with a relatively low mean connectivity, and the following thresholds were chosen for each analysis: 14 for WT at 8 h p.i., 9 for WT at 24 h p.i., 7 for WT + Ifnar^-/- at 8 h p.i., and 16 for WT + Ifnar^-/- at 24 h p.i. A topological overlap matrix (TOM) was then calculated from the adjacency matrix, converted to distance and clustered by hierarchical clustering using average linkage clustering. Modules were identified by a dynamic tree cut with a minimum module size = 20. Module eigengenes (MEs) that are the first principal component of the module were calculated, and similar modules were merged together using an ME distance threshold of 0.08 for WT at 8 h p.i., 0 for WT at 24 h p.i., 0.18 for WT + Ifnar^-/- at 8 h p.i., and 0.08 for WT + Ifnar^-/- at 24 h p.i. (Supplementary Figures S2–S5). kWithin that is the intramodular connectivity and kTotal, the whole network connectivity of each gene, were calculated with the WGCNA package. The relationship of MEs with the infection status (independent variable) was assessed with a regression analysis on MEs (dependent variable) for each module separately. The WT Medium condition was used as a reference (intercept). Adjusted R-squared and p-values were used to assess the model performance. ME average predictions were plotted as a heatmap. Adjusted R-squared (R-squared adjusted for the number of predictors in the model) that is a statistical measure that represents the proportion of the variance for a dependent variable (expression) that is explained by an independent variable (infection groups) was calculated by fitting a linear model for each gene expression value. Gene-adjusted R-squared was then plotted against kTotal values. Gene Ontology (GO) enrichment analysis was performed for gene coexpression modules against GO categories using the topGO R package (topGO 2.26.0) (Alexa and Rahnenfuhrer, 2016) and GO database (07.2019)  (Ashburner et al., 2000; Gene Ontology Consortium, 2021). Closeness centrality (CC) (Freeman, 1978) in the coexpression network was calculated for each gene at 8- and 24-h time points using the Networkx Python package (Networkx 2.5). A density map was plotted of CC at 8 h against CC at 24 h. Network visualization was performed with Cytoscape.

Results

The global network analysis of wild-type-infected macrophages highlighted modules correlated to either Leishmania, LRV1 or to Leishmania but further exacerbated by LRV1

To design appropriate strategies to combat leishmaniasis and understand its clinical presentations it was important to unravel the influence of Leishmania or LRV1 on the response mounted by the macrophage upon infection by LgyLRV1+. To this end, transcriptome profiling using RNA-Seq was performed to identify global changes in the WT murine macrophage infected with Lgy, carrying or not LRV1, defined as LRV1+ or LRV1-, at 8 and 24 h p.i. As a control, cells were left non-infected (Medium) or treated with polyinosinic-polycytidylic acid (poly I:C), a TLR3 agonist, mimicking solely the effect of dsRNA viruses such as LRV1. Of the 36,487 genes annotated for Mus musculus (McGarvey et al., 2015), 12,651 and 12,594 genes were detected in this study at 8 and 24 h p.i., respectively, and used for downstream analysis as shown in Figure 1. First, a weighted gene coexpression network analysis (WGCNA) was performed. WGCNA is a systems biology approach that relies on the hypothesis that genes with similar expression patterns may be functionally related, coregulated, or members of the same pathway. Following this approach, a gene–gene similarity network was constructed allowing to identify highly correlated genes that cluster into modules. The expression profile of the module genes was summarized by the ME that is the first principal component of its expression. Upon identification of the modules a regression analysis was performed on modules eigengenes and the average prediction was calculated to identify modules that were associated with the different experimental groups (Medium, LgyLRV1+, LgyLRV1- and poly I:C). The intramodular (kWithin) and whole network (kTotal) connectivity were computed for each gene allowing to identify hubs (drivers). Following the WGCNA, 49 modules at 8 h and 49 at 24 h were defined. The different modules are represented as a heatmap based on the average predictions of the regression analysis in the different experimental groups in Figures 2A, C (p-values shown in Supplementary Figures S6A, B).

FIGURE 1

Figure 1 Workflow of the bioinformatics analysis. WGCNA was performed first on the wild-type (WT) samples alone, then on WT + Ifnar^-/- samples merged. Regression analysis was performed on obtained modules to assess their relationship to phenotypes (infection groups). Gene Ontology (GO) enrichment was performed for each module. WT + Ifnar^-/- network resulting from weighted gene correlation network analysis (WGCNA) was used to calculate the closeness centrality (CC) of genes at 8 and 24 h postinfection (p.i.).

FIGURE 2

Figure 2 The global network analysis of WT-infected macrophages highlights modules and key pathways associated to Leishmania and to Leishmania RNA virus 1 (LRV1). (A) Heatmap of the average predictions of the fitted linear model on each module eigengene (ME) at the 8-h time point. (B) Network generated from selected modules at the 8-h time point associated with Leishmania infection (LgyLRV1+, LgyLRV1-), virus (LgyLRV1+, poly I:C), and “exacerbatory” modules. Only the top five highest connected genes were selected. Node colors indicate the module color they belong to. Edges between genes indicate the correlation between genes. (C) Heatmap of average predictions of the fitted linear model on each ME at the 24-h time point. (D) Network generated from selected modules at the 24-h time point associated with Leishmania infection (LgyLRV1+, LgyLRV1-), virus (LgyLRV1+, poly I:C), and “exacerbatory” modules. Only the top five highest connected genes were selected. Node colors indicate the module color they belong to. Edges between genes indicate the correlation between genes.

The blue to red scale speaks for negative to positive associations with a determined condition based on normalized enrichment scores. On top of each heatmap a dendrogram, representing the hierarchical clustering of the data, is shown. Interestingly, although Leishmania is the common denominator between LgyLRV1+ and the LgyLRV1- groups, LgyLRV1- clusters correlated closer to the group left untreated (mentioned as “Medium” in Figure 2) while LgyLRV1+ clusters correlated closer to the poly I:C group both at 8 and 24 h posttreatment. This observation demonstrates that the presence of LRV1 within Leishmania triggers a drastic macrophage response, shifting its profile closer to the response triggered solely by a virus than to that of Leishmania alone.

Based on the analysis of the heatmaps, the modules were defined as either Leishmania- or LRV1 dependent. Comparatively to the medium condition, the LRV1-modulated group encompassed modules exclusively up- or downregulated by LgyLRV1+ and poly I:C but not modulated by Leishmania devoid of LRV1 (i.e., LgyLRV1-). On the other hand, the Leishmania-modulated group was characterized by modules similarly regulated by LgyLRV1- and LgyLRV1+ comparatively to medium, but not by poly I:C. Among the 49 modules at each time point, and 5 modules were found to be associated with Leishmania at 8 and 24 h p.i., respectively, and 6 and 9 modules were found to be associated with the presence of LRV1 at 8 and 24 h p.i., respectively. Thus, the modules selected for further analysis due to their association with Leishmania were lightcyan (245 genes), grey60 (255 genes), magenta4 (70 genes), floralwhite (112 genes) at 8 h p.i. and mediumpurple3 (116 genes), green (493 genes), plum1 (134 genes), grey60 (185 genes), paleturquoise (148 genes) at 24 h p.i. The modules selected due to their association to LRV1 were blue2 (725 genes), firebrick3 (159 genes), coral3 (644 genes), navajowhite1 (130 genes), lavenderblush3 (145 genes), plum (83 genes) at 8 h p.i. and black (345 genes), skyblue3 (135 genes), lightsteelblue1 (116 genes), blue (1207 genes), brown (1137 genes), darkgrey (168 genes), red (352 genes), tan (237 genes), darkslateblue (86 genes) at 24 h p.i. Given the known exacerbatory role that LRV1 has on the leishmaniasis outcome (Ives et al., 2011; Hartley et al., 2012; Cantanhede et al., 2015; Rossi et al., 2017) another category of modules was defined as “exacerbatory”. These so-called “exacerbatory” modules are, to a certain extent modified by LgyLRV1-, but the presence of LRV1 in LgyLRV1+ further exacerbates this response in a similar direction as that observed with the poly I:C treatment. At 8 h p.i., there were 7 “exacerbatory” modules highlighted: green (1,159 genes), orangered4 (435 genes), coral1 (744 genes), darkseagreen4 (344 genes), bisque4 (1,129 genes), darkgrey (884 genes), and mediumpurple4 (158 genes). At 24 h p.i., 6 “exacerbatory” modules: darkolivegreen (145 genes), pink (294 genes), yellow (639 genes), magenta (287 genes), brown4 (99 genes), and ivory (111 genes) were identified. (For the whole lists of genes with their module membership and their connectivity at each time point, see Supplementary Table S1).

Based on the above-mentioned selection of modules, the presence of Leishmania lead to a concerted modulation of 682 genes at 8 h p.i. and 1,076 genes at 24 h p.i. representing 5.4% and 8.5%, respectively, of the total genes identified. On the other hand, the presence of LRV1 impacted 1,886 genes at 8 h p.i. and 3,783 genes at 24 h p.i. representing 14.9% and 30% of the total amount of genes identified at 8 and 24 h p.i., respectively. As for the “exacerbatory” modules, they comprise a large proportion of the genes with a synchronized modulation, 4853 and 1575 genes representing 38.4% and 12.5% of the genes identified at 8 and 24 h p.i., respectively. This approach gave us a general overview of the macrophage response toward Leishmania, LRV1, or the combination of both agents. The results support the importance of the modulation of the macrophage response by LRV1 both at 8 and 24 h p.i.

The global network analysis of wild-type-infected macrophages highlighted key pathways associated to Leishmania and to LRV1

Using the WGCNA groups of genes with a similar expression were identified and classified into modules; however, to shed light into the biological processes that may be present in each of the modules selected above, a GO enrichment analysis using the topGO R package (topGO 2.26.0) (Alexa and Rahnenfuhrer, 2016) and the GO database (Ashburner et al., 2000; Gene Ontology Consortium, 2021) was performed. The aim was to identify annotated functions overrepresented in the modules. Among the biological process (BP) category of the GO terms identified, a curation of the GOs were conducted and GO terms related to either clearly different types of cells or tissues, or even to completely distinct organs or organisms were removed. The 5 most statistically significant terms (with p-values strictly lower than 0.01) within a module are shown in Supplementary Table S2. The fact that WGCNA analysis resulted in modules enriched for biologically important processes related to infection, including an innate immune response (found in the 8h_coral3 module) and a positive regulation of interleukin-6 production (found in the 8h_bisque4 module), suggest that these modules are a robust feature of the molecular architecture of Leishmania and LRV1 infection. Modules associated to Leishmania were enriched in terms such as the positive regulation of telomere capping, RNA secondary structure unwinding, purine deoxyribonucleotide, and glutathione metabolic process in the grey60, floralwhite, lightcyan, and magenta4 modules, respectively, at 8 h p.i. (Supplementary Table S2). Moreover, at 24 h p.i., the Leishmania-associated modules were enriched in terms such as positive regulation of receptor internalization and asymmetric protein localization in mediumpurple3 and plum1, respectively, and the regulation of transcription, negative regulation of focal adhesion assembly, and cellular response to acid chemical in paleturquoise, grey60 and green, respectively (Supplementary Table S2). Modules associated with LRV1 at 8 h p.i. were enriched in pathways related to tRNA methylation, transcription, translation, GTPase activity, and DNA replication, in the blue2, firebrick3, navajowhite1, and lavenderblush3 modules at 8 h p.i., whereas coral3 and plum were enriched in the immune response (for example, innate immune response and a positive regulation of acute inflammatory response) and phagocytosis pathways. Comparatively, LRV1-associated modules at 24 h p.i., such as skyblue3, blue, lightsteelblue1, darkgrey, and tan, were enriched in pathways related to methylation, translation, and transcription processes, whereas the black, brown, red, tan, and darkslateblue modules were enriched in phagocytosis, cytokine production, and the oxidation–reduction process (Supplementary Table S2). On the other hand, the “exacerbatory” modules were enriched in terms such as the RNA biosynthetic process, oxidation–reduction process and mitosis-related pathways in green, orangered4 and darkseagreen4; and NF-kappaB signaling and cytokine production in coral1 and bisque4 at 8 h p.i. Furthermore, the “exacerbatory” modules at 24 h p.i. were enriched in DNA replication, transcription, and macrophage regulation in the darkolivegreen, yellow, magenta, and brown4 modules (Supplementary Table S2).

Identification of hub genes and network analysis provided insight on how the different modules may interact

In each module, the most central genes, that is, the most interconnected genes measured by their intramodular connectivity, that is, its kWithin, can be further identified as hub genes. Due to their high connectivity, hub genes are thus considered as functionally important genes that are most likely to drive the group of genes and thus the biological processes present within each module. The top 1% of genes with the highest kWithin are shown in Table 2 for each of the modules selected. In order to visualize the interaction between the Leishmania virus and “exacerbatory” modules, the five genes with the highest connectivity for each module (i.e., with the highest kWithin) are displayed as a network for 8 and 24 h p.i. in Figures 2B, D. The five top driver genes associated to Leishmania were Ift27, Cdk5, Surf1, Trf, Ndufb11 (lightcyan), Ift46, Selenbp1, Cyp4f16, Dhrs1, Haus4 (grey60), Dynlrb1, Gm9843, Gm10076, Rps13-ps2, Gm15772 (magenta4), Memo1, Esd, Gm42517, Psmd12, Cct3 (floralwhite) at 8 h p.i.; and Ddx39b, Ptpmt1, Rnh1, Glmp, Wsb2 (mediumpurple3), Arrdc3, Clk1, Map3k8, Hyal1, Slc16a6 (green), Zfx, Sp1, Stag1, Dusp11, Pan3 (plum1), St3gal3, Klf16, Asrgl1, Rdh10, Irf2bp1 (grey60), Hook3, Ubxn7, Ankrd44, Brd1, Nrip1 (paleturquoise) at 24 h p.i. The five top driver genes inside the modules selected due to their association to LRV1 were Cd3eap, Cebpa, Oit3, Ckap4, Chek2 (blue2), Nudcd2, C1qbp, Anp32b, Fam212a, Adprhl2 (firebrick3), 5430427O19Rik, Sdc4, Tmem243, Nupr1, Lpxn (coral3), Tmem131, Lnpep, Man2a1, Sptbn1, Nckap1l (navajowhite1), Jarid2, Mtus1, Oxct1, Abcb1a, Sephs1 (lavenderblush3), Clasp2, Pros1, Hsdl2, Uap1, Rap1gds1 (plum) at 8 h p.i.; and Chek2, Myo5a, Trem2, Cnrip1, Mamdc2 (black), Cc2d1b, Dnase2a, Decr2, Pafah2, Lamtor4 (skyblue3), Rack1, Eif3i, Rps18, Eri3, Rpl6 (lightsteelblue1), Pygl, Gsn, Prkag2, Rpl3, Aldh9a1 (blue), 1600014C10Rik, Rnf19b, Plcl2, Abtb2, Il15ra (brown), Patl1, Prkaa1, Thrap3, Ppp1r13b, Med13 (darkgrey), Tom1l2, Mllt6, Sec24b, Dnajc1, Tspan15 (red), Al504432, Trps1, Dock4, Flnb, Foxn2 (tan), Plin3, Pla2g12a, Ptrh1, Trabd, Ssbp4 (darkslateblue) at 24 h p.i. For the “exacerbatory” modules, the 5 top driver genes were: Arhgap4, Zdhhc9, Mvb12b, Prkar2a, Kdm2b (green), Tsnax, Mtfr1l, Tmem126b, Cmas, Mri1 (orangered4), Phlpp1, Crybg1, Fmnl2, Fam53c, Prkx (coral1), BC031181, Rfc2, Rps24, Cox7a2, Stx8 (darkseagreen4), Irf1, Lhx2, Stx11, Cd40, Cxcl10 (bisque4), Rora, Camsap2, Ehmt1, Yeats2, Sacs (darkgrey) and 9930111J21Rik1, Caprin1, Trappc8, Zdhhc5, Pstpip2 (mediumpurple4) at 8 h p.i.; and Acbd5, Brca2, Slc16a7, Sgpp1, Rabgap1l (darkolivegreen), Nomo1, Marveld1, Pfas, Cad, Tjp2 (pink), Igf1r, Abi2, B3glct, 2610507B11Rik, Cbx6 (yellow), Camsap2, Kank2, Mib1, Rora, Mef2a (magenta), Tmem40, Gtf2f1, Mgea5, Hmgb1-ps5, Igsf11 (brown4), Glce, Rasal2, Tmem135, Wdr6, Tprn (ivory) and Slc31a2, Gpb3, Isg20, Cd274, Nmi (turquoise) at 24 h p.i. (for the whole lists of genes with their module membership and their connectivity at each time point, see Supplementary Table S1).

TABLE 2

Table 2 The top 1% genes with the highest kWithin of the three groups of modules selected in WT analysis at 8 and 24 h postinfection (p.i).

As shown, the top five hub genes of each module were clustered diametrically, depending on whether they were either positively or negatively modulated by a given treatment as shown in Figure 2B. Interestingly, at 8 h, the connection between these two extremities was achieved by a module with an exacerbatory phenotype, the 8h_darkseagreen4, enriched in mitosis-related pathways as mentioned before (Supplementary Table S2). At 24 h, the connection between modules that were positively or negatively modulated is accomplished by two different branches. Interestingly, the connecting modules at 24 h show different modulations for either Leishmania (LRV1+ and LRV1-) or poly I:C (i.e., the molecule used to map virus-dependent responses). On one side, the connection was realized by the 24h_green module consisting of upregulated genes in Leishmania-infected macrophages but downregulated in poly I:C-treated cells. As mentioned previously, the 24h_green module was enriched in the cellular response to the acid chemical (Supplementary Table S2). On the other side, the link between modules with up- and downregulated genes was achieved by the 24h_grey60 and 24h_plum1 showing both the downregulation of genes associated with Leishmania but being either upregulated or non-modulated by poly I:C. As mentioned previously, the 24h_grey60 module was enriched in genes implicated in the negative regulation of focal adhesion assembly and the 24h_plum1 in asymmetric protein localization (Supplementary Table S2). The topology of the networks represented by the five top genes from each module suggested the existence of a dependency between the different cascades triggered by Leishmania and LRV1 infection.

The global network analysis of wild-type-infected macrophages highlighted modules with highly connected genes that explained most of the variance of the data (highly adjusted R-squared)

Following the approach described above four, six, and seven modules were found associated to Leishmania, LRV1, and the combination of both, respectively, at 8 h p.i., whereas five, nine, and six modules were associated to Leishmania, LRV1, and the combination of both, respectively, at 24 h p.i. However, in order to understand which modules had potentially the highest impact on the macrophage transcriptional response, genes were plotted according to their total connectivity (kTotal) and adjusted R-squared (Figure 3). This approach allowed us to identify core modules. Thus, the genes present at the tip of the plot belonging to the bisque4, coral1, and green modules at 8 h and brown, turquoise, red, and blue modules at 24 h were likely to be strongly affected by Leishmania and LRV1 infection. Interestingly, all modules at 8 h could be considered to have an exacerbatory modulation, while at 24 h, three out of the four modules were virus driven (brown, red, and blue) while turquoise was exacerbatory. This observation indicated that while at early time points, the combination of Leishmania and LRV1 had the highest impact on the macrophage response, this effect was mainly achieved later by the presence of LRV1, suggesting that the presence of Leishmania at later time points impacts the macrophage response to a much lesser extent than LRV1. Moreover, out of the seven modules present at the tip, five (8h_bisque4; 8h_coral1; 24h_brown; 24h_turquoise and 24h_red; with a total of 1,872 genes at 8 h and 2,810 genes at 24 h) were enriched in pathways associated with the immune response (Table 3), highlighting the importance of a coordinated immune response and its overall impact when macrophages were challenged by an external agent. Among these immune-enriched pathways, there were, for example, the GO terms: innate immune response, cellular responses to or the regulations of IFNs or cytokines, defense responses to virus or protozoan (Table 3), whereas the other non-immune modules present at the tip were enriched in, for example, the GO terms: RNA biosynthetic and metabolic processes, the regulation of transcription from the RNA polymerase II promoter and oxidation–reduction process (Table 3). These results showed how immune response–related genes shape the network topology, emphasizing its impact on the transcriptomic profile of the macrophage and highlighting the strong and sustained effect that LRV1 had on the macrophage response.

FIGURE 3

Figure 3 The global network analysis of WT-infected macrophages highlights modules with highly connected genes that explain most of the variance of the data (highly adjusted R-squared). (A) Scatter plot of kTotal (whole network connectivity) against adjusted R-squared for all genes in an 8-h network. Genes are colored according to the module (described in Figure 2) they belong to. (B) Scatterplot of kTotal (whole network connectivity) against adjusted R squared for all genes in the 24 h network. Genes are colored according to the module (described in Figure 2) they belong to.

TABLE 3

Table 3 Examples of main Gene Ontology (GO) terms of tip modules in WT analysis (biological process “BP” category and p-value < 0.01) at 8 and 24 h p.i.

Type I interferons played a preponderant and central role in the infection mounted by macrophages toward LgyLRV1+

Given the importance of the type I IFN response to the immunophenotype observed upon infection with LgyLRV1+ (Rossi et al., 2017) coupled with the observation that modules containing genes with the highest kTotal and adjusted R-squared (8h_bisque4 and 24h_turquoise modules) are enriched in the cellular response to IFN-β (Figures 3A, B, Table 3) led us to further explore the systemic role of type I IFNs for gene modulation. To achieve this, a WGCNA approach was conducted in the Ifnar knockout macrophages (Ifnar^-/-). Moreover, to directly test the role of type I IFNs, two additional treatments, IFN-α and IFN-β, coadministered after 6 h post- LgyLRV1-infection, were tested to specifically mimic the type I IFN response induced by the Leishmania virus LRV1 on the modulation of genes. The analysis generated 38 modules at 8 h p.i. and 31 modules at 24 h p.i., comprising a total of 11,442 and 11,426 genes, respectively (Figures 4A, B, p-values shown in Supplementary Figures S6C, D). As expected, the dendrogram displayed on top of each heatmap showed that modules representing the infection of Ifnar^-/- cells with LgyLRV1+ clustered closer to the modules of WT cells infected with LgyLRV1-, than to the modules of WT cells infected with LgyLRV1+, confirming the major impact that the recognition of type I IFNs had on the macrophages infected by LgyLRV1+. At 8 h p.i. (Figure 4A), the addition of type I IFNs to LgyLRV1- did not lead to a drastic effect as these groups clustered closer together with non-treated cells than with LgyLRV1+. This result was not surprising as treatment with type I IFNs was only performed 6 h p.i. with LgyLRV1- and thus cells were only exposed to type I IFNs for 2 h before sample collection. Nevertheless, four modules where the addition of type I IFNs showed an effect by shifting LgyLRV1- closer together to LgyLRV1+ and to poly I:C were identified: 8h_greenyellow (2,226 genes), 8h_mediumpurple3 (1,870 genes), 8h_firebrick4 (479 genes), and 8h_purple (424 genes). These modules contained genes whose expression was promptly modulated by type I IFNs, either exogenously added or produced in response to LgyLRV1+. The heatmap for 24 h (Figure 4B) also showed the role of IFNAR in gene modulation as modules representing Ifnar^-/- cells infected with LgyLRV1+ clustered closer to modules representing WT infected with LgyLRV1- rather than LgyLRV1+. Furthermore, at 24 h, cells infected with LgyLRV1- cluster together with non-treated cells. The addition of type I IFNs to LgyLRV1- at 24 h leads to a new cluster further away from non-treated cells contrary to what was observed at 8 h, demonstrating the long-term effect of the type I IFN response. Modules in which the effect of type I IFNs was clearly observed were 24h_tan (579 genes), 24h_thistle2 (744 genes), 24h_thistle1 (680 genes), 24h_lightyellow (266 genes), 24h_darkturquoise (166 genes), 24h_lightsteelblue1 (1,020 genes), and 24h_magenta (1207 genes). Based on the total number of genes detected at 8 and 24 h (11,442 and 11,426, respectively), the percentages of genes modulated by type I IFNs at 8 and 24 h were 43.7% and 40.8%, respectively, showing the drastic impact that type I IFNs have on the concerted modulation of genes.

FIGURE 4

Figure 4 Type I IFNs play a preponderant and central role in the infection mounted by macrophages toward LgyLRV1+. (A) Heatmap of the average predictions of the fitted linear model on each ME at the 8-h time point in WT + Ifnar^-/- analysis. (B) Heatmap of average predictions of the fitted linear model on each ME at the 24-h time point in WT + Ifnar^-/- analysis. (C) Scatter plot of kTotal (whole network connectivity) against adjusted R-squared for all genes in an 8-h network in WT + Ifnar^-/- analysis. Genes are colored according to the module they belong to. (D) Scatter plot of kTotal (whole network connectivity) against adjusted R-squared for all genes in a 24-h network in WT + Ifnar^-/- analysis. Genes are colored according to the module they belong to.

Following the approach described above, to determine and narrow down the importance of modules, the kTotal (the connectivity of the gene in the whole network) and adjusted R-squared (the proportion of the variance for expression that is explained by infection groups) of each gene were calculated at the two time points and plotted on the x-axis and y-axis, respectively (Figures 4C, D). From this analysis, two preponderant modules could explain the variance of the data at 8 h p.i. (modules with the highest kTotal and the highest adjusted R-squared, localized at the tip of the graph), namely, 8h_greenyellow and 8h_mediumpurple3 (Figure 4C), and four modules were found at the tip of the graph at 24 h p.i., namely, 24h_thistle1, 24h_thistle2, 24h_lightsteelblue1, and 24h_magenta (Figure 4D). Out of these six modules, three were enriched in pathways connected to the immune response, namely, 8h_greenyellow, 24h_thistle1, and 24h_thistle2. These modules were enriched in the following pathways: defense response to virus and protozoan and response to IFN-gamma, -beta, and -alpha; interleukin-6; and tumor necrosis factor productions, whereas the three other ones were enriched in RNA transcription and histone acetylation pathways, as well as in oxidation–reduction and metabolic processes (Supplementary Table S3). Curiously, one of the modules identified, the 8h_greenyellow module, was constituted of 2,224 genes, approximately one-fifth of all genes detected at 8 h in the RNA-Seq conducted, showing not only the importance but also the coordination of the early immune response in macrophages upon a challenge. Given the influence of hub genes in orchestrating the response of the modules, the top 1% genes with the highest kWithin of these two and four modules were examined at 8 and 24 h p.i. (Table 4). Interestingly, as examples Cxcl11, Tgtp1, and Il27 (in the greenyellow module) at the 8-h time point; Ifih1 and Igtp (in the thistle1 module), Il15ra and Jak2 (in the thistle2 module) at the 24-h time point were found (Table 4) (for the whole lists of genes with their module membership and their connectivity at each time point, see Supplementary Table S1).

TABLE 4

Table 4 The top 1% genes with the highest kWithin of the preponderant modules that explain the variance of the data in WT + Ifnar^-/- analysis at 8 and 24 h p.i.

Overlap of highly connected modules at early and late time points uncovered the temporal dynamics of the interferon response

The greenyellow module identified at 8 h that represented a large proportion of the genes identified is also enriched in different Gos, including many related to immunity. However, at 24 h p.i., more than one module was enriched in Gos connected to immunity (Supplementary Table S3). To understand the temporal dynamics of the macrophage response, the overlap between the major module identified at 8 h and the modules identified at 24 h p.i. in the WT + Ifnar^-/- analysis was evaluated. Interestingly, the upset graph (Figure 5A) showed that the 8h_greenyellow module overlapped very strongly with 24h_thistle1 and 24h_thistle2, 24h_tan, 24h_darkmagenta, and 24h_lightyellow. As previously mentioned, the striking GO terms enriched in immune-related Gos were found in 8h_greenyellow, 24h_thistle1, and 24h_thistle2. The predominant terms identified in the 8h_greenyellow module were also found at 24-h split into the 24h_thistle1 module (e.g., defense response to virus and protozoan, innate immune response, and cellular response to IFN-gamma and -beta) and 24h_thistle2 (e.g., a positive regulation of the ERK1 and ERK2 cascades, interleukin-12, and chemokine productions) (Supplementary Table S3). Intriguingly, other Gos that demonstrated a concerted modulation in 8h_greenyellow and were thus associated with initial immune responses were now split into 24h_tan (e.g., a negative regulation of the NF-kappaB transcription factor activity), 24h_darkmagenta (e.g., a positive regulation of the inflammatory response) and 24h_lightyellow (e.g., the regulation of cytokine secretion and the MyD88-dependent Toll-like receptor signaling pathway) (Supplementary Table S3 and data not shown). This latter part was interesting in terms of Gos that can be initially associated with the immune response toward invaders but were then separated. Remarkably, some other GO terms of the 8h_greenyellow module were not found significantly (p-value < 0.01) in any of the five modules cited above and identified in Figure 5A (e.g., a negative regulation of interleukin-10 production). These results pointed toward the dynamic nature of the response, being more concerted at 8 h and then split into the groups of genes with slightly different modulations and thus ending up in different modules.

FIGURE 5

Figure 5 Overlap of highly connected modules at early and late time points uncovers the temporal dynamics of the interferon response. (A) An UpSet plot of the 8h_greenyellow module overlapping with the modules at the 24-h time point in WT + Ifnar-/- analysis. Top five overlapping modules are shown. Intersection size is shown in the y-axis. The bottom-right part shows the total module size. (B) Density plot of the CC of genes in the WGCNA network at 8 h (x-axis) against 24 h (y-axis) in WT + Ifnar-/- analysis. Count unit corresponds to the number of genes in each rectangle. (C) Zoom of the tip of the CC plot. Genes with very high centrality at both 8 and 24 h p.i. (D) Scatter plot of CC of genes in the WGCNA network at 8 h (x-axis) against 24 h (y-axis) in WT + Ifnar-/- analysis. Positions of ISGs are highlighted in red and the names of 10 examples are shown.

To further dissect the temporal dynamic of the macrophage response toward infection, the CC of each gene in the network was calculated at 8 and 24 h p.i., then plotted on x- and y-axes, respectively. This measure allowed us to identify genes that were best placed to influence the entire network as regulators and important signal transducers. The different categories of genes regarding their position on the CC plot are shown in Figure 5B. Strikingly, there were genes with very low centrality (below 0.05) at 8 h (on the y-axis) and at 24 h (on the x-axis), 76 genes versus 1567 genes, respectively (approximately 20 times more), suggesting that these genes were essential at 24 h versus 8 h p.i., respectively. The great difference in the number of genes in these two groups showed that many central events could happen at 8 h, right at the beginning of the infection, losing its importance later on or setting the stage for the later events. At 8 h, the predominant GO terms were related to the RNA polymerase II pathway such as the positive regulation of transcription and transcription initiation from the RNA polymerase II promoter (Figure 6). To sum up, 23 terms were enriched in RNA polymerase II processes in comparison to only 12 terms related to both RNA polymerase I and III pathways (for example, a positive regulation of the transcription from the RNA polymerase I promoter and RNA polymerase III promoter, respectively) (Figure 6). In contrast, no relevant GO term was found enriched in the set of 76 genes with high CC at 24 h p.i. and very low at 8 h p.i. On the CC plot, the blue-to-yellow scale (density) represents the number of genes per position. A “hotspot” with the highest density, localized on the x-axis (i.e., with very low CC at a 24-h time point), with a CC value approximately 0.6 at an 8-h time point was pinpointed. This “hotspot” of genes belonged to the 1,567 genes mentioned before (with very low centrality at 24 h p.i.).

FIGURE 6

Figure 6 Genes from RNA polymerase II pathway are predominantly central at early time point. Density plots of the CC of genes in the WGCNA network at 8 h (x-axis) against 24 h (y-axis) in WT + Ifnar^-/- analysis. Genes belonging to the examples of RNA polymerase I (A), II (B), and III (C) processes are highlighted in red. The lists of GO terms found for RNA polymerase I (Table A), II (Table B), and III (Table C) keywords are listed.

The tip of the CC plot was defined by genes with very high centrality both at 8 and at 24 h p.i. The genes present at the tip of the CC plot belonged to the modules 8h_greenyellow, 8h_lightgreen, and most importantly 8h_mediumpurple3 (Figure 5C). As mentioned previously, this last module was globally enriched in GO terms linked to transcription, meaning that these pathways initially influenced the network the most. While there were three modules from the 8-h modules, the genes at the tip were divided into more than 10 modules from the 24-h modules (higher than four times more). For example, the tip of the CC plot included Akt1, a member of the protein kinase B family, known to be involved in macrophage survival and parasite persistence (Eren et al., 2016).

Given the high and central importance of type I IFNs in the modulation of the transcript response and the outcome of an infection by LgyLRV1+, a non-exhaustive list of genes stimulated by type I IFNs was established from the GO database and some other references and their position on the CC plot evaluated (Figure 5D). Genes from three GO terms were compiled: cellular response to IFN-alpha (GO:0035457), cellular response to IFN-beta (GO:0035458), and cellular response to type I IFN (GO:0071357). In addition, genes cited in different publications were also considered (Zhao et al., 2009; Schoggins and Rice, 2011; Rand et al., 2012; Cheon et al., 2014; Manjunath et al., 2017; Nasr et al., 2017; Pervolaraki et al., 2018; Ashley et al., 2019; Yang et al., 2020). Thus, the list comprised a total of 175 genes analyzed at 8 and 24 h p.i. in WT + Ifnar^-/- analysis. Among them, 139 and 138 genes were detected at 8 and 24 h p.i., respectively. The consistency of detection between early and late time points showed a sustained induction of IFN-stimulated genes (ISGs). However, there was one gene detected at 8 h p.i. but not at 24 h p.i.: Ifi202b (Supplementary Table S4). The membership status of the whole list of detected genes was then investigated, and, as expected, the main modules related to immune pathways were predominant. More than half of the ISGs detected at 8 and 24 h p.i. belong to the 8h_greenyellow (96 genes) and its overlapping 24h_thistle1 and 24h_thistle2 modules (70 + 20 = 90 genes), respectively. Nevertheless, the early gene, Ifi202b, belonged to the 8h_blue2 module (Supplementary Table S4). To confront these ISGs to the CC analysis, the 175 genes were sought, and the detected ones are highlighted (in red in Figure 5D). The five top genes were Myd88, Tgtp1, Mov10, Trim21, and Samhd1 showing the highest CC values at 8 h p.i. (>0.72) and Irak1, Sun2, Ttll12, Mmp12, and Ube2g2 with the highest CC values at 24 h p.i. (>0.37) (Figure 5D, Supplementary Table S4). Ifi202b has a CC value of 0.60 at 8 h p.i. and is localized on the x-axis; thus, it is part of the “hotspot” already defined before (Figure 5B). Even if detected at 24 h p.i., another ISG with a very low CC value (below 0.05) at 24 h p.i. belonged to this “hotspot” and was found on the x-axis: Ifitm1 (8h_CC value of 0.61) (Figures 5B, D, Supplementary Table S4). Ifitm1 was part of the 8h_plum1 module, while, as expected at 24 h p.i., it was included in 24h_grey, a label reserved to genes that was not part of any module (Supplementary Table S4). Therefore, these two genes, Ifi202b and Ifitm1, seemed to be very central at the early time point, suggesting that they were both important at 8 h p.i. and not at 24 h p.i.

Discussion

Many studies have addressed, using RNA-Seq, the transcriptional changes in the macrophage response upon Leishmania infection (Salloum et al., 2021). Traditionally, these studies have investigated differentially expressed genes with significantly altered expression across the groups of samples. However, moving from a gene-by-gene analysis to function is very challenging. In that regard, network analysis methods, such as WGCNA, have the advantage to provide a comprehensive view at the system level. As an unsupervised algorithm, WCGNA can establish and detect the relationship between gene expression and phenotypic traits. Moreover, WGCNA has the advantage over other methods in analyzing multivariate and complex datasets (Langfelder and Horvath, 2008). Indeed, WGCNA allows to condense a large network of genes into a limited number of variables (i.e., modules). Furthermore, WGCNA permits the identification of hub genes that play a central role in driving the system of interest. While an WGCNA approach has already been conducted in the blood of patients (Gardinassi et al., 2016) and the popliteal lymph nodes of dogs (Sanz et al., 2021) infected with L. infantum as well as in the cutaneous lesions of L. braziliensis–infected patients (Christensen et al., 2016) and L. major–infected human dendritic cells (Zhao et al., 2019), no WGCNA investigation has been piloted on Leishmania-infected macrophages. Aside from being, to our knowledge, the first study where WGCNA was conducted in Leishmania-infected macrophages in vitro, this study also represents the first comprehensive transcriptional study addressing the response mounted by macrophages upon infection by LgyLRV1+, thus addressing the impact the presence of LRV1 has on the macrophage response (Cantacessi et al., 2015; Salloum et al., 2021). The present WGCNA identified biologically relevant groups of transcripts, classified into modules, that were modified upon macrophage infection by LgyLRV1- and LgyLRV1+ after 8 or 24 h of infection. These modules were further analyzed and key pathways and hub genes associated with Leishmania, with LRV1 or with the exacerbatory effect of LRV1 in LgyLRV1+, were identified. Additionally, the impact of type I IFNs on the transcriptional profile of macrophages upon LgyLRV1+ infection was evaluated.

The major impact of LRV1 and the relevance of type I interferons

Globally, this study showed that the presence of LRV1 had a major impact in the transcriptomic profile of macrophages both at 8 and 24 h p.i. The LRV1 and the “exacerbatory” modules, that is, modules reflecting the effect of the presence of LRV1 on top of the presence of Leishmania, were identified and shown to represent a very considerable part of the overall transcriptomic change (representing together 53.3% and 42.5% at 8 and 24 h p.i.). The immunophenotype associated to the presence of LRV1, or to other viruses, was shown to be highly dependent on the systemic production of type I IFNs (Rossi et al., 2017; Rossi and Fasel, 2018b; Rath et al., 2019; Heirwegh et al., 2021). Although type I IFNs, such as IFN-α and IFN-β, are crucial for viral clearance, their influence during a Leishmania infection appears to be highly contextual. Type I IFNs have a worsening effect on the outcome of L. guyanensis, while in the case of L. major, the addition of type I IFNs promotes the resolution of infection (Hartley et al., 2014). Here, we show the massive role that the production of type I IFNs downstream of the recognition of LRV1 has at a transcriptional level influencing 43.7% and 40.8% of all transcriptional changes at 8 and 24 h p.i. Type I IFNs are produced by different cell types including macrophages and can act both in an autocrine and paracrine manner. The recognition of type I IFN by the IFNAR receptor (a complex composed of two subunits: IFNAR1 and IFNAR2) initiates a signaling cascade, leading to the expression of a vast panel of ISGs (IFN-stimulated genes) with a positive feedback loop and also direct antiviral effects (Sadler and Williams, 2008). Many studies have investigated the induction of some ISGs in response to a viral infection, for example, adenovirus, morbillivirus, retrovirus, vesiculovirus, herpesvirus, and paramyxovirus infections (Zhao et al., 2009; Manjunath et al., 2017; Nasr et al., 2017; Pervolaraki et al., 2018; Ashley et al., 2019; Yang et al., 2020). Interestingly, some of the hub genes found in the WT and WT + Ifnar^-/- analyses (Figure 2, Table 4) were part of the ISGs such as 9930111J21Rik1, Cxcl10, Cxcl11, Ifih1, Igtp, Irf1, Jak2, and Tgtp1 in WT + Ifnar^-/- analysis (Supplementary Table S4). The chemokine CXCL10 has been demonstrated to be linked to the endosymbiotic LRV1 (Ives et al., 2011; Kariyawasam et al., 2017) and the Leishmania evasion from the host (Antonia et al., 2019).

TNF-α and IL-6 regulation as evidence of the crosstalk between modules

Hyperinflammation is a hallmark of an LgyLRV1+ infection (Ives et al., 2011). Two major mediators of inflammation contributing to this immunophenotype are the cytokines IL-6 and TNF-α. Surprisingly, these two molecules were found present in two distinct modules: IL-6 was found in the 8h_coral 3 (an LRV1-associated module), while TNF-α was found associated to the 8h_bisque4 modules (an “exacerbatory” module of WT analysis). In addition to TNF-α, CXCL10 (one of the top five hub genes), and CCL5, two chemokines also found to be the hallmarks of the inflammation caused by LgyLRV1+ infection were also found to be present in the 8h_bisque4. Interestingly, even though IL-6 is not found in the 8h_bisque4, this module was shown to be enriched in the GO term “positive regulation of IL-6.” Further analysis of both 8h_bisque4 and 8h_coral3 modules revealed that the members of the Src-family kinase (SFK) were also found in these two modules. Src and Hck were also found present in the 8h_bisque 4 together with TNF-α, whereas the other members of the SFK-family such as Lyn and Lck were found present in the 8h_coral3 module, the module that contains IL-6. This specific distribution of SFK is particularly relevant since, from our previous studies, it was shown that KB-SRC4, a specific inhibitor of c-Src, blocks IL-6 but not TNF-α (Brandvold et al., 2012; Reverte et al., 2021). On the contrary, PP2, which efficiently blocks Fyn, Hck, and Lck, inhibits TNF-α but not IL-6 secretion (Hanke et al., 1996; Reverte et al., 2021). Thus, it is interesting to underline not only the relevance of SFKs in TNF-α and IL-6 secretion but also the crosstalk between two modules containing genes upregulated in LgyLRV1+-infected cells. Interestingly, Csk, the kinase that regulates the intramolecular inhibition of SFKs, was also identified in the 8h_coral3 module. This module also includes important genes, such as IL-1 or IFN-γ previously associated with the presence of LRV1 (Ives et al., 2011; Hartley et al., 2016; Rossi et al., 2017; Kariyawasam et al., 2017; Castiglioni et al., 2017) as well as other genes Il18, Ccl8, Socs 2, Casp3, Il10, Jun, Jak3, Bcl10, and Ly6e implicated in the response to Leishmania infection. Overall, these data highlight a physiological and phenotypically relevant crosstalk between two different modules. Moreover, since the role of IL-6 has been considered, by Osero et al., one of the greatest oxymora in the leishmaniasis outcome due to different studies showing discordant roles (Osero et al., 2020), this crosstalk also pointed toward the relevance that a module associated in this study to an LRV1 module may contribute to the leishmaniasis outcome.

IL-15 as a possible important hub

As mentioned above, the WGCNA analysis identifies hub genes that play a central role in a given module, that is, genes most likely to influence the whole module and thus drive its enriched pathways. Among the hubs identified, one gene, found in both analyses, is Il15ra, encoding for the receptor of the interleukin, IL-15 (Figure 2, Table 4). This cytokine, in combination with IL-12, stimulates the cellular immune response in dogs and humans with visceral leishmaniasis (Milano et al., 2002; Costa et al., 2020). Similarly, another cytokine-related gene was observed as a driver gene but only in WT + Ifnar^-/-, for example, Il27. IL-27 is associated with resistance or susceptibility to Leishmania infection (Jafarzadeh et al., 2020) and to the antiviral response to several RNA viruses, such as HIV, Hepatitis C Virus (HCV), or Chikungunya Virus (CHIKV) (Valdés-López et al., 2022). Interestingly, Il15, Il15ra, Il12a, Il12b, and Il27 genes were all present in the same “exacerbatory” module at 8 h p.i. (8h_bisque4), then separated into the 24h_brown (for Il15 and Il15ra) and 24h_turquoise (for Il12b and Il27) modules, associated to LRV1 and exacerbatory phenotypes in WT analysis.

Activation of specific transcription factors in response to LRV1

Given the described importance of LRV1 in the leishmaniasis outcome, progression, relapse, and treatment failure, hub genes associated to the presence of LRV1 are of utmost interest when designing strategies to fight leishmaniasis (Hartley et al., 2014; Adaui et al., 2016; Bourreau et al., 2016; Rossi et al., 2017). Several relevant transducers and transcription factors were present in the 8h_coral3 and 8h_bisque4 modules in WT analysis (e.g., Myd88, Keap1, Tank, Ikbke, Stat1, Stat2, Irf1, Irf5, Irf7, Ets2, and Xbp1). For example, Irf1 was part of the “exacerbatory” 8h_bisque4 module. IRF1 participates in the expression of cytokines such as CCL5 and CXCL10. The lack of IRF1 shows a dramatically exacerbated leishmaniasis disease in L. major–infected mice associated with a decrease of IFN-γ and IL-12 productions (Lohoff et al., 1997). IRF1 is also involved in antiviral responses against dsRNA viruses, for example, HCV (Pflugheber et al., 2002; Feng et al., 2017). Furthermore, Igf1r was a driver gene in WT analysis and belonged to the 24h_yellow module and the “exacerbatory” group. Similarly, the insulin-like growth factor 1 receptor (IGF1R) has previously been shown to be involved in arginase (Arg1) expression in visceral leishmaniasis (Osorio et al., 2014) and also targeted by rotaviruses, dsRNA viruses, to manipulate the PI3K/Akt pathway and block autophagy (Zhou et al., 2018). More intriguingly, several transcription factors were found among the list of these hub genes in WT analysis, such as Sp1. Sp1 that belonged to the 24h_plum1 is involved in the evasion mechanisms of Leishmania. Different species can induce the binding of Sp1 on Il10, Ucp2, and Hdac1 promoters in macrophages, resulting in the production of the inhibitory IL-10 cytokine, the suppression of oxidative burst, and a decrease of iNOS expression, respectively (Yang et al., 2007; Basu Ball et al., 2014; Calegari-Silva et al., 2018).

Interferon-stimulated genes and temporal expression

Some studies have analyzed different time points following an infection challenge and revealed that IFN-stimulated genes (ISGs) can be temporally expressed: earlier or later, in order to specifically target and counteract different steps of the virus life cycle (Zhao et al., 2009; Manjunath et al., 2017; Nasr et al., 2017; Pervolaraki et al., 2018). The analysis of the closeness centrality of some ISGs brought two genes to the forefront, Ifi202b and Ifitm1, which were shown to be central at an 8-h time point only (Figure 5, Supplementary Table S4). Therefore, these two candidates could be decisive in the early events happening during the cellular response and could represent new targets to further investigate. On one hand, Ifi202b, also known as IFN-activated gene 202B and encoding for the P202 protein, belongs to the pyrin and HIN domain–containing (PYHIN) proteins, as AIM2 and IFI16 (Wang et al., 2018). In several studies, Ifi202b has been associated with sex differences in autoimmune diseases (Panchanathan et al., 2011; Cao et al., 2018) and this question has recently been addressed in the Leishmania inflammation context (Lockard et al., 2019; Snaka et al., 2022). Likewise, Ifi202b has been shown to be upregulated upon L. major footpad infection in C57BL/6 mice (Ehrchen et al., 2010) and its corresponding P202 protein to exert an inhibitory effect on the AIM2 inflammasome (Wang et al., 2018). On the other hand, Ifitm1, encoding the IFN-induced transmembrane protein 1 (IFITM1), is known to have a role in the restriction of many RNA viruses, such as HIV (Chutiwitoonchai et al., 2013) or HCV (Narayana et al., 2015).

Toward strategies to face Leishmania-virus coinfections and beyond

The apparent disproportional response toward a non-pathogenic element such as LRV1 has the potential to change the outcome of the infection by a pathogenic element such as Leishmania. Although LRV1 is only present in some strains of Leishmania, the results presented here can have a broader impact: the modulation of leishmaniasis by endogenous or exogenous viruses is particularly relevant as the vector of Leishmania, the phlebotomine sand flies, also carries phleboviruses such as Toscana virus shown previously to cause the same phenotype as the presence of LRV1 when coadministered with Leishmania. Additionally, the effect that virus can have go beyond LRV1 or the coadministration of phlebovirus during a sand fly bite. Mice challenged with a virus, such as LCMV, upon the resolution of Leishmania infection, led to a relapse of the Leishmania infection showing that exacerbation and the metastatic phenotype are not always linked to the presence of LRV1 (Valencia et al., 2022) but could also be related to the presence of other viruses (Rossi et al., 2017; Rath et al., 2019; Heirwegh et al., 2021). By exposing the modules and pathways affected by LRV1 as well as the hub genes responsible for driving such events, this study can contribute to the future design of strategies to deal with a Leishmania-viral coinfection.

The RNA-Seq data used in this manuscript are available to all communities through the website https://amelbek.shinyapps.io/fasel_lab_data/ where individual or groups of genes can be visualized in different genotypes of macrophages (as the data for this study were part of larger datasets containing additional mouse genotypes IFNAR-/-, IFNg-/-, iNOS-/-, NLRX1-/-, NOX2-/-, and PRX5-/-). The data can be downloaded for further statistical analysis. Thus, this study and these datasets also constitute a resource that can be further explored to study the impact of Leishmania on the macrophage response and on the importance of viral coinfections or type I IFNs and beyond in the outcome of leishmaniasis.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/, GSE201120, https://www.ncbi.nlm.nih.gov/geo/, GSE203088.

Ethics statement

The animal study was reviewed and approved by Veterinary Commission of the Canton de Vaud (SCAV, Switzerland).

Author contributions

NI, FT, TS, CD, DK, FP, and MR performed experiments. AB analyzed the data. FT, AB, and NI wrote the manuscript. FT, AB, NI, IX, and NF interpreted and discussed the data. TS, SC, IX, and NF reviewed and edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the grant from the Swiss National fund for research to NF (Grant No. 310030_173180) and the grant from Foundation Pierre Mercier pour la science to FT.

Acknowledgments

The authors thank the veterinary and animal facility staff at CIIL for ensuring animal welfare and ethical standards. Library preparation, sequencing and raw data normalization were performed at the Lausanne Genomic Technologies Facility (GTF), University of Lausanne (UNIL), Switzerland (https://wp.unil.ch/gtf/). The authors specially thank the following members of the facility for their support: Leonore Wigger, Sandra Calderon, Julien Marquis, Hannes Richter, Johann Weber, Roberto Sermier, and also the former members, Sylvain Pradervand, Floriane Consales Barras and Keith Harshman. The authors thank Matteo Rossi for his help with the initiation of the study and Alexander Miesch for his precious assistance and patience in making the data available. The authors are very thankful to Slavica Masina for her critical reading.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2022.941888/full#supplementary-material

References

Adaui, V., Lye, L. F., Akopyants, N. S., Zimic, M., Llanos-Cuentas, A., Garcia, L., et al. (2016). Association of the endobiont double-stranded RNA virus LRV1 with treatment failure for human leishmaniasis caused by leishmania braziliensis in Peru and Bolivia. J. Infect. Dis. 213 (1), 112–121. doi: 10.1093/infdis/jiv354

PubMed Abstract | CrossRef Full Text | Google Scholar

Alexa, A., Rahnenfuhrer, J. (2016). topGO: Enrichment analysis for gene ontology. r package version. Bioconductor 2(0). doi: 10.18129/B9.bioc.topGO

CrossRef Full Text | Google Scholar

Anders, S., Pyl, P. T., Huber, W. (2015). HTSeq–a Python framework to work with high-throughput sequencing data. bioinformatics 31 (2), 166–169. doi: 10.1093/bioinformatics/btu638

PubMed Abstract | CrossRef Full Text | Google Scholar

Antonia, A. L., Gibbs, K. D., Trahair, E. D., Pittman, K. J., Martin, A. T., Schott, B. H, et al. (2019). Pathogen evasion of chemokine response through suppression of CXCL10. Front. Cell. Infect. Microbiol. 9, 280–280. doi: 10.3389/fcimb.2019.00280

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoki, J. I., Muxel, S. M., Zampieri, R.A ., Laranjeira-Silva, M. F., Muller, K. E., Nerland, A. H, et al. (2017). RNA-Seq transcriptional profiling of leishmania amazonensis reveals an arginase-dependent gene expression regulation. PloS Negl. Trop. Dis. 11 (10), e0006026. doi: 10.1371/journal.pntd.0006026

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoki, J. I., Muxel, S. M., Zampieri, R. A., Muller, K. E., Nerland, A. H., Floeter-Winter, L. M. (2019). Differential immune response modulation in early leishmania amazonensis infection of BALB/c and C57BL/6 macrophages based on transcriptome profiles. Sci. Rep. 9 (1), 19841. doi: 10.1038/s41598-019-56305-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M, et al. (2000). Gene ontology: tool for the unification of biology. the gene ontology consortium. Nat. Genet. 25 (1), 25–29. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashley, C. L., Abendroth, A., McSharry, B. P., Slobedman, B. (2019). Interferon-independent innate responses to cytomegalovirus. Front. Immunol. 10, 2751. doi: 10.3389/fimmu.2019.02751

PubMed Abstract | CrossRef Full Text | Google Scholar

Basu Ball, W., Mukherjee, M., Srivastav, S., Das, P. K. (2014). Leishmania donovani activates uncoupling protein 2 transcription to suppress mitochondrial oxidative burst through differential modulation of SREBP2, Sp1 and USF1 transcription factors. Int. J. Biochem. Cell Biol. 48, 66–76. doi: 10.1016/j.biocel.2014.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bourreau, E., Ginouves, M., Prevot, G., Hartley, M. A., Gangneux, J. P., Robert-Gangneux, F., et al. (2016). Presence of leishmania RNA virus 1 in leishmania guyanensis increases the risk of first-line treatment failure and symptomatic relapse. J. Infect. Dis. 213 (1), 105–111. doi: 10.1093/infdis/jiv355

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandvold, K. R., Steffey, M. E., Fox, C. C., Soellner, M. B. (2012). Development of a highly selective c-src kinase inhibitor. ACS Chem. Biol. 7 (8), 1393–1398. doi: 10.1021/cb300172e

PubMed Abstract | CrossRef Full Text | Google Scholar

Burza, S., Croft, S. L., Boelaert, M. (2018). Leishmaniasis. Lancet 392 (10151), 951–970. doi: 10.1016/S0140-6736(18)31204-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Calegari-Silva, T. C., Vivarini, Á.C., Pereira, R. d. M. S., Dias-Teixeira, K. L., Rath, C. T., Pacheco, A. S. S., et al. (2018). Leishmania amazonensis downregulates macrophage iNOS expression via histone deacetylase 1 (HDAC1): a novel parasite evasion mechanism. Eur. J. Immunol. 48 (7), 1188–1198. doi: 10.1002/eji.201747257

PubMed Abstract | CrossRef Full Text | Google Scholar

Cantacessi, C., Dantas-Torres, F., Nolan, M.J., Otranto, D. (2015). The past, present, and future of leishmania genomics and transcriptomics. Trends Parasitol. 31 (3), 100–108. doi: 10.1016/j.pt.2014.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Cantanhede, L. M., da Silva Junior, C. F., Ito, M. M., Felipin, K. P., Nicolete, R., Salcedo, J. M, et al. (2015). Further evidence of an association between the presence of leishmania RNA virus 1 and the mucosal manifestations in tegumentary leishmaniasis patients. PloS Negl. Trop. Dis. 9 (9), e0004079. doi: 10.1371/journal.pntd.0004079

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Wang, L., Wang, C.-Y., Ye, J., Wang, Y., Li, T., et al. (2018). Sex differences in correlation with gene expression levels between Ifi200 family genes and four sets of immune disease-relevant genes. J. Immunol. Res. 2018, 1290814–1290814. doi: 10.1155/2018/1290814

PubMed Abstract | CrossRef Full Text | Google Scholar

Castiglioni, P., Hartley, M. -A., Rossi, M., Prevel, F., Desponds, C., Utzschneider, D. T, et al. (2017). Exacerbated leishmaniasis caused by a viral endosymbiont can be prevented by immunization with its viral capsid. PloS Negl. Trop. Dis. 11 (1), e0005240–e0005240. doi: 10.1371/journal.pntd.0005240

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandrakar, P., Parmar, N., Descoteaux, A., Kar, S. (2020). Differential induction of SOCS isoforms by leishmania donovani impairs macrophage-T cell cross-talk and host defense. J. Immunol. 204 (3), 596–610. doi: 10.4049/jimmunol.1900412

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaparro, V., Graber, T.E., Alain, T., Jaramillo, M. (2022). Transcriptional profiling of macrophages reveals distinct parasite stage-driven signatures during early infection by leishmania donovani. Sci. Rep. 12 (1), 6369. doi: 10.1038/s41598-022-10317-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheon, H., Borden, E. C., Stark, G. R. (2014). Interferons and their stimulated genes in the tumor microenvironment. Semin. Oncol. 41 (2), 156–173. doi: 10.1053/j.seminoncol.2014.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, S. M., Dillon, L. A., Carvalho, L. P., Passos, S., Novais, F. O., Hughitt, V. K, et al. (2016). Meta-transcriptome profiling of the human-leishmania braziliensis cutaneous lesion. PloS Negl. Trop. Dis. 10 (9), e0004992. doi: 10.1371/journal.pntd.0004992

PubMed Abstract | CrossRef Full Text | Google Scholar

Chutiwitoonchai, N., Hiyoshi, M., Hiyoshi-Yoshidomi, Y., Hashimoto, M., Tokunaga, K., Suzu, S. (2013). Characteristics of IFITM, the newly identified IFN-inducible anti-HIV-1 family proteins. Microbes Infect. 15 (4), 280–290. doi: 10.1016/j.micinf.2012.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, S. F., Gomes, V. O., Dos Santos Maciel, M. O., Melo, L. M., Venturin, G. L., Bragato, J. P, et al. (2020). Combined in vitro IL-12 and IL-15 stimulation promotes cellular immune response in dogs with visceral leishmaniasis. PloS Negl. Trop. Dis. 14 (1), e0008021–e0008021. doi: 10.1371/journal.pntd.0008021

PubMed Abstract | CrossRef Full Text | Google Scholar

da Silva, R. L., Elizondo, D. M., Brandy, N. Z. D., Haddock, N. L., Boddie, T. A., de Oliveira, L. L, et al. (2021). Leishmania donovani infection suppresses allograft inflammatory factor-1 in monocytes and macrophages to inhibit inflammatory responses. Sci. Rep. 11 (1), 946. doi: 10.1038/s41598-020-79068-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, M. P. A., van Dongen, S., Abreu-Goodger, C., Bartonicek, N., Enright, A. J. (2013). Kraken: a set of tools for quality control and analysis of high-throughput sequence data. Methods (San Diego Calif.) 63 (1), 41–49. doi: 10.1016/j.ymeth.2013.06.027

CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S, et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29 (1), 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Ehrchen, J. M., Roebrock, K., Foell, D., Nippe, N., von Stebut, E., Weiss, J. M., et al. (2010). Keratinocytes determine Th1 immunity during early experimental leishmaniasis. PloS Pathog. 6 (4), e1000871–e1000871. doi: 10.1371/journal.ppat.1000871

PubMed Abstract | CrossRef Full Text | Google Scholar

Eren, R. O., Reverte, M., Rossi, M., Hartley, M. A., Castiglioni, P., Prevel, F, et al. (2016). Mammalian innate immune response to a leishmania-resident RNA virus increases macrophage survival to promote parasite persistence. Cell Host Microbe 20 (3), 318–328. doi: 10.1016/j.chom.2016.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, H., Lenarcic, E. M., Yamane, D., Wauthier, E., Mo, J., Guo, H., et al. (2017). NLRX1 promotes immediate IRF1-directed antiviral responses by limiting dsRNA-activated translational inhibition mediated by PKR. Nat. Immunol. 18 (12), 1299–1309. doi: 10.1038/ni.3853

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandes, M. C., Dillon, L. A., Belew, A. T., Bravo, H. C., Mosser, D. M., El-Sayed, N. M.. (2016). Dual transcriptome profiling of leishmania-infected human macrophages reveals distinct reprogramming signatures. mBio 7 (3), 00027–ndash;16. doi: 10.1128/mBio.00027-16

CrossRef Full Text | Google Scholar

Ferrari, C. K. B. (2019). “Oxidative stress and antioxidants in host defense in leishmaniasis,” in Oxidative stress in microbial diseases (Singapore: Springer), 245–256.

Google Scholar

Freeman, L. C. (1978). Centrality in social networks conceptual clarification. Soc. Networks 1 (3), 215–239. doi: 10.1016/0378-8733(78)90021-7

CrossRef Full Text | Google Scholar

Gardinassi, L. G., Garcia, G. R., Costa, C. H., Costa Silva, V., de Miranda Santos, I. K. (2016). Blood transcriptional profiling reveals immunological signatures of distinct states of infection of humans with leishmania infantum. PloS Negl. Trop. Dis. 10 (11), e0005123. doi: 10.1371/journal.pntd.0005123

PubMed Abstract | CrossRef Full Text | Google Scholar

Gene Ontology C. (2021). The gene ontology resource: enriching a GOld mine. Nucleic Acids Res. 49 (D1), D325–D334. doi: 10.1093/nar/gkaa1113

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, G., Oghumu, S., Satoskar, A. R. (2013). Mechanisms of immune evasion in leishmaniasis. Adv. Appl. Microbiol. 82, 155–184. doi: 10.1016/B978-0-12-407679-2.00005-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanke, J. H., Gardner, J. P., Dow, R. L., Changelian, P. S., Brissette, W. H., Weringer, E. J, et al. (1996). Discovery of a novel, potent, and src family-selective tyrosine kinase inhibitor. study of lck- and FynT-dependent T cell activation. J. Biol. Chem. 271 (2), 695–701. doi: 10.1074/jbc.271.2.695

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartley, M. A., Ronet, C., Zangger, H., Beverley, S. M., Fasel, N. (2012). Leishmania RNA virus: when the host pays the toll. Front. Cell. Infect. Microbiol. 2. doi: 10.3389/fcimb.2012.00099

CrossRef Full Text | Google Scholar

Hartley, M.-A., Drexler, S., Ronet, C., Beverley, S.M., Fasel, N. (2014). The immunological, environmental, and phylogenetic perpetrators of metastatic leishmaniasis. Trends Parasitol. 30 (8), 412–422. doi: 10.1016/j.pt.2014.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartley, M.-A., Bourreau, E., Rossi, M., Castiglioni, P., Eren, R. O., Prevel, F, et al. (2016). Leishmaniavirus-dependent metastatic leishmaniasis is prevented by blocking IL-17A. PloS Pathog. 12 (9), e1005852–e1005852. doi: 10.1371/journal.ppat.1005852

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartley, M. A., Eren, R. O., Rossi, M., Prevel, F., Castiglioni, P., Isorce, N., et al. (2018). Leishmania guyanensis parasites block the activation of the inflammasome by inhibiting maturation of IL-1beta. Microb. Cell 5 (3), 137–149. doi: 10.15698/mic2018.03.619

PubMed Abstract | CrossRef Full Text | Google Scholar

Heirwegh, E., MacLean, E., He, J., Kamhawi, S., Sagan, S. M., Olivier, M. (2021). Sandfly fever Sicilian virus-leishmania major co-infection modulates innate inflammatory response favoring myeloid cell infections and skin hyperinflammation. PloS Negl. Trop. Dis. 15 (7), e0009638. doi: 10.1371/journal.pntd.0009638

PubMed Abstract | CrossRef Full Text | Google Scholar

Ives, A., Ronet, C., Prevel, F., Ruzzante, G., Fuertes-Marraco, S., Schutz, F, et al. (2011). Leishmania RNA virus controls the severity of mucocutaneous leishmaniasis. Science 331 (6018), 775–8. doi: 10.1126/science.1199326

PubMed Abstract | CrossRef Full Text | Google Scholar

Jafarzadeh, A., Nemati, M., Chauhan, P., Patidar, A., Sarkar, A., Sharifi, I, et al. (2020). Interleukin-27 functional duality balances leishmania infectivity and pathogenesis. Front. Immunol. 11, 1573–1573. doi: 10.3389/fimmu.2020.01573

PubMed Abstract | CrossRef Full Text | Google Scholar

Kariyawasam, R., Grewal, J., Lau, R., Purssell, A., Valencia, B. M., Llanos-Cuentas, A, et al. (2017). Influence of leishmania RNA virus 1 on proinflammatory biomarker expression in a human macrophage model of American tegumentary leishmaniasis. J. Infect. Dis. 216 (7), 877–886. doi: 10.1093/infdis/jix416

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhlmann, F. M., Robinson, J. I., Bluemling, G. R., Ronet, C., Fasel, N., Beverley, S. M.. (2017). Antiviral screening identifies adenosine analogs targeting the endogenous dsRNA leishmania RNA virus 1 (LRV1) pathogenicity factor. Proc. Natl. Acad. Sci. U.S.A. 114 (5), E811–E819. doi: 10.1073/pnas.1619114114

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., Horvath, S. (2008). WGCNA: an r package for weighted correlation network analysis. BMC Bioinf. 9, 559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

Li, B., Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinf. 12 (1), 323. doi: 10.1186/1471-2105-12-323

CrossRef Full Text | Google Scholar

Lindoso, J. A. L., Cunha, M. A., Queiroz, I. T., Moreira, C. H. V. (2016). Leishmaniasis-HIV coinfection: current challenges. HIV/AIDS (Auckland N.Z.) 8, 147–156. doi: 10.2147/HIV.S93789

CrossRef Full Text | Google Scholar

Lockard, R. D., Wilson, M. E., Rodriguez, N. E. (2019). Sex-related differences in immune response and symptomatic manifestations to infection with leishmania species. J. Immunol. Res. 2019, 4103819. doi: 10.1155/2019/4103819

PubMed Abstract | CrossRef Full Text | Google Scholar

Lohoff, M., Ferrick, D., Mittrücker, H. -W., Duncan, G. S., Bischof, S., Röllinghoff, M, et al. (1997). Interferon regulatory factor-1 is required for a T helper 1 immune response in vivo. Immunity 6 (6), 681–689. doi: 10.1016/S1074-7613(00)80444-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Manjunath, S., Mishra, B. P., Mishra, B., Sahoo, A. P., Tiwari, A. K., Rajak, K. K, et al. (2017). Comparative and temporal transcriptome analysis of peste des petits ruminants virus infected goat peripheral blood mononuclear cells. Virus Res. 229, 28–40. doi: 10.1016/j.virusres.2016.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J. 17 (1), 10–12. doi: 10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

McGarvey, K. M., Goldfarb, T., Cox, E., Farrell, C. M., Gupta, T., Joardar, V. S, et al. (2015). Mouse genome annotation by the RefSeq project. Mamm. Genome Off. J. Int. Mamm. Genome Soc. 26 (9-10), 379–390. doi: 10.1007/s00335-015-9585-8

CrossRef Full Text | Google Scholar

Milano, S., Di Bella, G., D'Agostino, P., Barbera, C., Caruso, R., La Rosa, M, et al. (2002). IL-15 in human visceral leishmaniasis caused by leishmania infantum. Clin. Exp. Immunol. 127 (2), 360–365. doi: 10.1046/j.1365-2249.2002.01749.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Narayana, S. K., Helbig, K. J., McCartney, E. M., Eyre, N. S., Bull, R. A., Eltahla, A, et al. (2015). The interferon-induced transmembrane proteins, IFITM1, IFITM2, and IFITM3 inhibit hepatitis c virus entry. J. Biol. Chem. 290 (43), 25946–25959. doi: 10.1074/jbc.M115.657346

PubMed Abstract | CrossRef Full Text | Google Scholar

Nasr, N., Alshehri, A. A., Wright, T. K., Shahid, M., Heiner, B. M., Harman, A. N, et al. (2017). Mechanism of interferon-stimulated gene induction in HIV-1-Infected macrophages. J. Virol. 91 (20), e00744–17. doi: 10.1128/JVI.00744-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Olivier, M., Zamboni, D. S. (2020). Leishmania viannia guyanensis, LRV1 virus and extracellular vesicles: a dangerous trio influencing the faith of immune response during muco-cutaneous leishmaniasis. Curr. Opin. Immunol. 66, 108–113. doi: 10.1016/j.coi.2020.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Osero, B. O., Aruleba, R. T., Brombacher, F., Hurdayal, R. (2020). Unravelling the unsolved paradoxes of cytokine families in host resistance and susceptibility to leishmania infection. Cytokine X 2 (4), 100043. doi: 10.1016/j.cytox.2020.100043

PubMed Abstract | CrossRef Full Text | Google Scholar

Osorio, E. Y., Travi, B. L., da Cruz, A. M., Saldarriaga, O. A., Medina, A. A., Melby, P. C. (2014). Growth factor and Th2 cytokine signaling pathways converge at STAT6 to promote arginase expression in progressive experimental visceral leishmaniasis. PloS Pathog. 10 (6), e1004165. doi: 10.1371/journal.ppat.1004165

PubMed Abstract | CrossRef Full Text | Google Scholar

Panchanathan, R., Duan, X., Arumugam, M., Shen, H., Liu, H., Choubey, D.. (2011). Cell type and gender-dependent differential regulation of the p202 and Aim2 proteins: implications for the regulation of innate immune responses in SLE. Mol. Immunol. 49 (1-2), 273–280. doi: 10.1016/j.molimm.2011.08.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Pearson, R. D., Anastacio de Queiroz, S. (1996). Clinical spectrum of leishmaniasis. Clin. Infect. Dis. 22 (1), 1–11. doi: 10.1093/clinids/22.1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Pervolaraki, K., Rastgou Talemi, S., Albrecht, D., Bormann, F., Bamford, C., Mendoza, J. L, et al. (2018). Differential induction of interferon stimulated genes between type I and type III interferons is independent of interferon receptor abundance. PloS Pathog. 14 (11), e1007420. doi: 10.1371/journal.ppat.1007420

PubMed Abstract | CrossRef Full Text | Google Scholar

Pflugheber, J., Fredericksen, B., Sumpter, R., Jr., Wang, C., Ware, F., Sodora, D. L, et al. (2002). Regulation of PKR and IRF-1 during hepatitis c virus RNA replication. Proc. Natl. Acad. Sci. United States America 99 (7), 4650–4655. doi: 10.1073/pnas.062055699

CrossRef Full Text | Google Scholar

Rand, U., Rinas, M., Schwerk, J., Nohren, G., Linnes, M., Kroger, A., et al. (2012). Multi-layered stochasticity and paracrine signal propagation shape the type-I interferon response. Mol. Syst. Biol. 8, 584. doi: 10.1038/msb.2012.17

PubMed Abstract | CrossRef Full Text | Google Scholar

Rath, C. T., Schnellrath, L. C., Damaso, C. R., de Arruda, L. B., Vasconcelos, P., Gomes, C, et al. (2019). Amazonian Phlebovirus (Bunyaviridae) potentiates the infection of leishmania (Leishmania) amazonensis: Role of the PKR/IFN1/IL-10 axis. PloS Negl. Trop. Dis. 13 (6), e0007500. doi: 10.1371/journal.pntd.0007500

PubMed Abstract | CrossRef Full Text | Google Scholar

Restrepo, C. M., Llanes, A., Herrera, L., Ellis, E., Lleonart, R., Fernandez, P. L. (2021). Gene expression patterns associated with leishmania panamensis infection in macrophages from BALB/c and C57BL/6 mice. PloS Negl. Trop. Dis. 15 (2), e0009225. doi: 10.1371/journal.pntd.0009225

PubMed Abstract | CrossRef Full Text | Google Scholar

Reverte, M., Eren, R. O., Jha, B., Desponds, C., Snaka, T., Prevel, F, et al. (2021). The antioxidant response favors leishmania parasites survival, limits inflammation and reprograms the host cell metabolism. PloS Pathog. 17 (3), e1009422. doi: 10.1371/journal.ppat.1009422

PubMed Abstract | CrossRef Full Text | Google Scholar

Reverte, M., Snaka, T., Fasel, N. (2022). The dangerous liaisons in the oxidative stress response to leishmania infection. Pathogens 11 (4), 409. doi: 10.3390/pathogens11040409

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinf. (Oxford England) 26 (1), 139–140. doi: 10.1093/bioinformatics/btp616

CrossRef Full Text | Google Scholar

Rossi, M., Castiglioni, P., Hartley, M. A., Eren, R. O., Prevel, F., Desponds, C, et al. (2017). Type I interferons induced by endogenous or exogenous viral infections promote metastasis and relapse of leishmaniasis. Proc. Natl. Acad. Sci. U.S.A. 114 (19), 4987–4992. doi: 10.1073/pnas.1621447114

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossi, M., Fasel, N. (2018a). How to master the host immune system? leishmania parasites have the solutions! Int. Immunol. 30 (3), 103–111. doi: 10.1093/intimm/dxx075

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossi, M., Fasel, N. (2018b). The criminal association of leishmania parasites and viruses. Curr. Opin. Microbiol. 46, 65–72. doi: 10.1016/j.mib.2018.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Sadler, A. J., Williams, B. R. G. (2008). Interferon-inducible antiviral effectors. Nat. Rev. Immunol. 8 (7), 559–568. doi: 10.1038/nri2314

PubMed Abstract | CrossRef Full Text | Google Scholar

Salloum, T., Tokajian, S., Hirt, R. P. (2021). Advances in understanding leishmania pathobiology: What does RNA-seq tell us? Front. Cell Dev. Biol. 9, 702240. doi: 10.3389/fcell.2021.702240

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanz, C. R., Miro, G., Sevane, N., Reyes-Palomares, A., Dunner, S. (2021). Modulation of host immune response during leishmania infantum natural infection: A whole-transcriptome analysis of the popliteal lymph nodes in dogs. Front. Immunol. 12, 794627. doi: 10.3389/fimmu.2021.794627

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoggins, J. W., Rice, C. M. (2011). Interferon-stimulated genes and their antiviral effector functions. Curr. Opin. Virol. 1 (6), 519–525. doi: 10.1016/j.coviro.2011.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Scorza, B. M., Carvalho, E. M., Wilson, M. E. (2017). Cutaneous manifestations of human and murine leishmaniasis. Int. J. Mol. Sci. 18 (6), 1296. doi: 10.3390/ijms18061296

CrossRef Full Text | Google Scholar

Shadab, M., Das, S., Banerjee, A., Sinha, R., Asad, M., Kamran, M., et al. (2019). RNA-Seq revealed expression of many novel genes associated with leishmania donovani persistence and clearance in the host macrophage. Front. Cell Infect. Microbiol. 9, 17. doi: 10.3389/fcimb.2019.00017

PubMed Abstract | CrossRef Full Text | Google Scholar

Shapouri-Moghaddam, A., Mohammadian, S., Vazini, H., Taghadosi, M., Esmaeili, S. A., Mardani, F., et al. (2018). Macrophage plasticity, polarization, and function in health and disease. J. Cell. Physiol. 233 (9), 6425–6440. doi: 10.1002/jcp.26429

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheu, K. M., Hoffmann, A. (2022). Functional hallmarks of healthy macrophage responses: Their regulatory basis and disease relevance. Annu. Rev. Immunol. 40 (1), 295–321. doi: 10.1146/annurev-immunol-101320-031555

PubMed Abstract | CrossRef Full Text | Google Scholar

Snaka, T., Bekkar, A., Desponds, C., Prevel, F., Claudinot, S., Isorce, N., et al. (2022). Sex-biased control of inflammation and metabolism by a mitochondrial nod-like receptor. Front. Immunol. 13, 882867. doi: 10.3389/fimmu.2022.882867

PubMed Abstract | CrossRef Full Text | Google Scholar

Stuart, K. D., Weeks, R., Guilbride, L., Myler, P. J. (1992). Molecular organization of leishmania RNA virus 1. Proc. Natl. Acad. Sci. U.S.A. 89 (18), 8596–8600. doi: 10.1073/pnas.89.18.8596

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarr, P. I., Aline, R.F., Jr., Smiley, B. L., Scholler, J., Keithly, J., Stuart, K., et al. (1988). LR1: a candidate RNA virus of leishmania. Proc. Natl. Acad. Sci. U.S.A. 85 (24), 9572–9575. doi: 10.1073/pnas.85.24.9572

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdés-López, J. F., Fernandez, G. J., Urcuqui-Inchima, S. (2022). Synergistic effects of toll-like receptor 1/2 and toll-like receptor 3 signaling triggering interleukin 27 gene expression in chikungunya virus-infected macrophages. Front. Cell Dev. Biol. 10, 812110–812110. doi: 10.3389/fcell.2022.812110

PubMed Abstract | CrossRef Full Text | Google Scholar

Valencia, B. M., Lau, R., Kariyawasam, R., Jara, M., Ramos, A. P., Chantry, M., et al. (2022). Leishmania RNA virus-1 is similarly detected among metastatic and non-metastatic phenotypes in a prospective cohort of American tegumentary leishmaniasis. PloS Negl. Trop. Dis. 16 (1), e0010162. doi: 10.1371/journal.pntd.0010162

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P.-H., Ye, Z. -W., Deng, J. -J., Siu, K. -L., Gao, W. -W., Chaudhary, V. (2018). Inhibition of AIM2 inflammasome activation by a novel transcript isoform of IFI16. EMBO Rep. 19 (10), e45737. doi: 10.15252/embr.201845737

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Wang, S., Li, W. (2012). RSeQC: quality control of RNA-seq experiments. Bioinformatics 28 (16), 2184–2185. doi: 10.1093/bioinformatics/bts356

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Arslan, M., Liu, X., Song, H., Du, M., Li, Y., et al. (2020). IFN-gamma establishes interferon-stimulated gene-mediated antiviral state against Newcastle disease virus in chicken fibroblasts. Acta Biochim. Biophys. Sin. (Shanghai) 52 (3), 268–280. doi: 10.1093/abbs/gmz158

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z., Mosser, D. M., Zhang, X. (2007). Activation of the MAPK, ERK, following &lt;em<Leishmania amazonensis&lt;/em< infection of macrophages. J. Immunol. 178 (2), 1077. doi: 10.4049/jimmunol.178.2.1077

PubMed Abstract | CrossRef Full Text | Google Scholar

Yorek, M. S., Poudel, B., Mazgaeen, L., Pope, R. M., Wilson, M. E., Gurung, P. (2019). Leishmania major degrades murine CXCL1 - an immune evasion strategy. PloS Negl. Trop. Dis. 13 (7), e0007533. doi: 10.1371/journal.pntd.0007533

PubMed Abstract | CrossRef Full Text | Google Scholar

Zamboni, D. S., Sacks, D. L. (2019). Inflammasomes and leishmania: in good times or bad, in sickness or in health. Curr. Opin. Microbiol. 52, 70–76. doi: 10.1016/j.mib.2019.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zamora, M., Guilbride, L., Sacks, L., Stuart, K. (2000). Phylogenetic analysis of the 5' subterminal region of isolates of leishmania RNA virus-1. Ann. Trop. Med. Parasitol. 94 (2), 123–133. doi: 10.1080/00034980057464

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, H., Boije, H., Granberg, F., Pettersson, U., Svensson, C. (2009). Activation of the interferon-induced STAT pathway during an adenovirus type 12 infection. Virology 392 (2), 186–195. doi: 10.1016/j.virol.2009.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Dang, Z., Wei, R., Gui, W., Zhang, Y., Chong, S. (2019). The effects of CCR7 and related signaling pathways on leishmania major -infected human dendritic cells. J. Cell. Physiol. 234 (8), 13145–13156. doi: 10.1002/jcp.27985

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Geng, P., Liu, Y., Wu, J., Qiao, H., Xie, Y., et al. (2018). Rotavirus-encoded virus-like small RNA triggers autophagy by targeting IGF1R via the PI3K/Akt/mTOR pathway. Biochim. Biophys. Acta (BBA) - Mol. Basis Dis. 1864 (1), 60–68. doi: 10.1016/j.bbadis.2017.09.028

CrossRef Full Text | Google Scholar