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.

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.

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.

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).

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.

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)].

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 (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.

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).

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.

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).

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).

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.

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.

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.

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).

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.

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.).

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.

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).

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.

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.

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).

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).

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.