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

Abstract

Leishmania RNA virus 1 (LRV1) is a double-stranded RNA virus found in some strains of the human protozoan parasite Leishmania, the causative agent of leishmaniasis, a neglected tropical disease. Interestingly, the presence of LRV1 inside Leishmania constitutes an important virulence factor that worsens the leishmaniasis outcome in a type I interferon (IFN)-dependent manner and contributes to treatment failure. Understanding how macrophages respond toward Leishmania alone or in combination with LRV1 as well as the role that type I IFNs may play during infection is fundamental to oversee new therapeutic strategies. To dissect the macrophage response toward infection, RNA sequencing was performed on murine wild-type and Ifnar-deficient bone marrow-derived macrophages infected with Leishmania guyanensis (Lgy) devoid or not of LRV1. Additionally, macrophages were treated with poly I:C (mimetic virus) or with type I IFNs. By implementing a weighted gene correlation network analysis, the groups of genes (modules) with similar expression patterns, for example, functionally related, coregulated, or the members of the same functional pathway, were identified. These modules followed patterns dependent on Leishmania, LRV1, or Leishmania exacerbated by the presence of LRV1. Not only the visualization of how individual genes were embedded to form modules but also how different modules were related to each other were observed. Thus, in the context of the observed hyperinflammatory phenotype associated to the presence of LRV1, it was noted that the biomarkers tumor-necrosis factor α (TNF-α) and the interleukin 6 (IL-6) belonged to different modules and that their regulating specific Src-family kinases were segregated oppositely. In addition, this network approach revealed the strong and sustained effect of LRV1 on the macrophage response and genes that had an early, late, or sustained impact during infection, uncovering the dynamics of the IFN response. Overall, this study contributed to shed light and dissect the intricate macrophage response toward infection by the Leishmania-LRV1 duo and revealed the crosstalk between modules made of coregulated genes and provided a new resource that can be further explored to study the impact of Leishmania on the macrophage response.

Keywords: Leishmania RNA virus 1 (LRV1); RNA sequencing (RNA-Seq); interleukin 6 (IL-6); macrophage; tumor-necrosis factor alpha (TNF-α); type I interferon (IFN); weighted gene coexpression network analysis (WGCNA).

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

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.

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.

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.

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.

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.