- Research article
- Open Access
Simultaneous transcriptional profiling of Leishmania major and its murine macrophage host cell reveals insights into host-pathogen interactions
BMC Genomics volume 16, Article number: 1108 (2015)
Parasites of the genus Leishmania are the causative agents of leishmaniasis, a group of diseases that range in manifestations from skin lesions to fatal visceral disease. The life cycle of Leishmania parasites is split between its insect vector and its mammalian host, where it resides primarily inside of macrophages. Once intracellular, Leishmania parasites must evade or deactivate the host's innate and adaptive immune responses in order to survive and replicate.
We performed transcriptome profiling using RNA-seq to simultaneously identify global changes in murine macrophage and L. major gene expression as the parasite entered and persisted within murine macrophages during the first 72 h of an infection. Differential gene expression, pathway, and gene ontology analyses enabled us to identify modulations in host and parasite responses during an infection. The most substantial and dynamic gene expression responses by both macrophage and parasite were observed during early infection. Murine genes related to both pro- and anti-inflammatory immune responses and glycolysis were substantially upregulated and genes related to lipid metabolism, biogenesis, and Fc gamma receptor-mediated phagocytosis were downregulated. Upregulated parasite genes included those aimed at mitigating the effects of an oxidative response by the host immune system while downregulated genes were related to translation, cell signaling, fatty acid biosynthesis, and flagellum structure.
The gene expression patterns identified in this work yield signatures that characterize multiple developmental stages of L. major parasites and the coordinated response of Leishmania-infected macrophages in the real-time setting of a dual biological system. This comprehensive dataset offers a clearer and more sensitive picture of the interplay between host and parasite during intracellular infection, providing additional insights into how pathogens are able to evade host defenses and modulate the biological functions of the cell in order to survive in the mammalian environment.
In order to establish an infection, intracellular pathogens need to adapt to their new environment to survive the innate and adaptive immune responses of the host. A number of human pathogens - including Leishmania sp., Trypanosoma cruzi, Mycobacterium tuberculosis, Toxoplasma gondii, Francisella tularensis, Legionella pneumophila, and Ehrlichia - have evolved mechanisms not only to evade the host immune system, but to infect the very immune cells that are recruited to clear an infection [1–7]. Leishmania is able to infect and to replicate within mammalian macrophages. It may serve as a model of intracellular infection of immune cells and can be used to study transcriptomic changes that take place in both the host and the pathogen over the course of an infection.
Leishmania major and related species are the causative agents of leishmaniasis, a group of diseases that vary in severity from self-healing skin lesions to disfiguring mucosal manifestations to fatal visceral disease. More than a million new cases are reported annually, mostly concentrated in the Middle East and Central and South America . The Leishmania life cycle is divided between its insect vector, the phlebotomine sand fly, and its mammalian host, where it resides primarily inside of macrophages although neutrophils, dendritic cells, and fibroblasts have also been implicated at various stages of infection [9–13]. Previous studies have shown that the parasite undergoes changes in morphology and alterations in cell surface components as it adapts to the intracellular environment [14–16]. Additionally, a CD4+ T helper (Th) type 1 response by the host leads to parasite killing while a Th2 response leads to parasite growth [17–19]. Less is known about the global changes that take place at the transcriptomic level in both the parasite and host over the course of an infection.
Once in host macrophages, Leishmania parasites rapidly transform into aflagellated amastigote forms, which are contained inside of parasitophorous vacuoles. The parasite enters host cells by receptor-mediated phagocytosis. It is thought to do so in a quiescent manner, failing to produce a significant oxidative burst and to activate the innate immune system [20–23]. Leishmania prevent their killing by altering cytokine expression (thereby influencing T cell responses), impeding antigen display by MHC class II molecules, and hindering nitric oxide production (reviewed in [24–26]).
Previous studies using microarrays have started to elucidate changes that occur within the parasite or within the host as infection occurs [27–47], but have so far not sought to look at the transcriptomes of both simultaneously and over the course of an intracellular infection. Studies on host response have identified genes that are differentially regulated upon infection with various Leishmania species, sometimes with opposing results. Some of these differences may be attributable to the parasite species and host systems used, the severity of the resulting infection, and the timepoints examined. Additionally, studies of Leishmania amastigotes have often used axenic cultures [35, 36, 38, 45, 47] or lesion-derived amastigotes [27, 28, 35, 37, 44, 48]. The former have been shown to significantly differ from the intracellular biological state [35, 45] while the latter contain a mixture of amastigotes at various timepoints post-infection, making it difficult to differentiate between changes that take place early versus later during intracellular infection.
In this study, we performed transcriptome profiling using RNA-seq to concurrently identify global changes in host and parasite gene expression that occur over the first 72 h of a synchronized intracellular infection of murine peritoneal macrophages with Leishmania major. A paired-end mRNA sequencing approach was used to allow high-confidence read mapping and transcript assembly. Collection of data from multiple biological replicates, the use of matched host control samples, careful statistical analysis of variation, and removal of batch effects provided us with a unique ability to detect biological differences between samples and timepoints with high confidence and sensitivity. Differential gene expression analysis enabled us to identify host and parasite genes that were modulated over time, and pathway and gene ontology analyses provided insights into the higher level processes activated during infection. This work builds on and improves existing expression datasets and provides insights into how Leishmania is able to evade host defenses and cause modulations in the host transcriptome in order to survive in the mammalian intracellular environment.
Results and discussion
Infection dynamics and global transcription patterns
Transcriptome profiling by RNA-seq was used to identify global changes in gene expression over the course of the first 72 h of an infection of murine macrophages with Leishmania major. Peritoneal macrophages from C57BL/6 mice were infected with L. major metacyclic promastigotes and samples collected at 4, 24, 48, and 72 h post-infection (hpi). The dynamics of the infection were monitored by counting the number of parasites per 100 macrophages (Fig. 1a). RNA sequencing was carried out for each sample and for matched uninfected controls, as well as for the L. major metacyclic promastigotes used for the infection. Over 2.4 billion sequence reads were collected across three independent experiments (biological replicates) and labeled as batches A-C (see Additional file 1).
Since there is little sequence conservation between mouse and L. major, we were able to unambiguously map reads from mouse and parasite RNAs from the mixed sample to their respective genomes. For uninfected samples, the percentage of reads mapping to the mouse genome ranged from 88.9 to 95.9 %. Greater heterogeneity was observed for the percentage of reads mapping to L. major in the infected samples, providing clues about the transcriptional activity of both parasite and host (Fig. 1b). For each batch, the percentage of reads mapping to the parasite decreased over the course of the infection (see Additional file 1). Interestingly, this decrease in the relative portion of parasite reads did not match the parasite load of the infected cells, which decreased sharply between 4 and 24 h for each batch and remained constant or slowly rose again from 24 through 72 h (Fig. 1a and Additional file 1). Possible explanations for the mismatch between the constant/increasing number of parasites and the decreasing percentage of parasite reads is that the individual parasites are less transcriptionally active as the infection progresses or that the mouse cells became more transcriptionally active over the course of the infection, thereby diluting the proportion of reads attributable to the parasite.
Samples from host and from parasite showed a similar distribution of per-gene read counts per sample, as visualized by box plots (Additional file 2: Figure S1), but striking differences were observed when comparing the two organisms. The median steady-state expression level was elevated in L. major compared to mouse (6.0 vs. 3.6 log2 counts per million) and L. major genes showed a much more narrow distribution (interquartile ranges of 5.6 to 6.5 for L. major and 0 to 6.8 for mouse). Additionally, when non- and lowly-expressed genes were removed from the datasets prior to differential expression analysis, this filter led to the removal of 10,548 genes from the mouse dataset, but only 7 genes from the parasite one.
The differences in the distribution of genes between mouse and L. major are consistent with differences in how each organism controls gene expression - mouse employs both transcriptional and post-transcriptional mechanisms to control the expression levels of individual genes, resulting in many genes that are not expressed and a wide dynamic range for those that are. This differs from Leishmania, which employs polycistronic transcription by RNA polymerase II at roughly the same rate across its genome [49, 50]. The resulting polycistronic mRNAs are split into their component mature mRNAs by coupled trans-splicing and polyadenylation events [51–53]. Steady-state mRNAs levels of individual genes are thus not determined by the amount of transcription, but by post-transcriptional processes such as RNA processing and degradation (reviewed in ). The dearth of lowly expressed L. major genes may indicate that very few mRNAs are completely degraded following polycistronic transcription.
Statistical assessment of biological replicates and batch effects
The use of multiple biological replicates necessitated the evaluation of the data to assess reproducibility and account for batch effects, i.e. experimental variation caused by sub-groups of measurements that are not related to the underlying biology of the system being studied. Previous analyses of high-throughput data, like those produced by RNA-seq, have indicated the need to assess and correct batch effects in order to prevent misinterpretation of results [55, 56]. In this study, experimental start date was used as a surrogate for batch when testing for differential expression.
We used principal component analysis (PCA) and Euclidean distance heatmap analysis to visualize the relationship between experimental datasets both prior to and after adjusting for batch effects (Fig. 2). The PCA plots revealed that samples from the same experimental condition (hpi and infection status) grouped together for the parasite and mouse macrophages demonstrating a high level of reproducibility between biological replicates. The dendrograms associated with the Euclidean distance heatmaps further illustrated this point, with like samples clustering most closely with one another. The grouping of samples by experimental condition was partially obscured when batch effects were not considered, illustrating the benefits of this approach (Additional file 2: Figure S2). Batch effects were therefore controlled for in the subsequent differential expression analyses by including experimental batch in the statistical models used.
The PCA and clustering analysis also suggested interesting biological relationships between the samples. The global profile of L. major gene expression changed over the time course of the experiment, moving from left to right across principal component (PC) 1 (Fig. 2a). The clear separation between promastigote and amastigote samples along PC1 highlights the differences between the two developmental stages. While a time progression can be seen within the intracellular amastigote samples, some intermixing of samples from different timepoints was observed beyond 4 hpi (Fig. 2c), suggesting an overlap in the gene expression profiles for these samples.
A time progression was also seen in the mouse macrophage data for both infected and uninfected samples (Fig. 2b and d). This observation underscored the necessity of collecting uninfected controls from each timepoint studied rather than relying on a control from a single timepoint. In addition, all infected macrophages from 24 to 72 hpi clustered more closely with the uninfected macrophages while the 4 hpi infected macrophages clustered apart. This is especially apparent in the heatmap dendrogram where the 4 hpi infected macrophages appear as a major outgroup.
Identification of differentially expressed genes and pathway-based enrichment analyses in murine host macrophages
While clustering analyses provided a high level overview of the behavior of the murine macrophage and parasite transcriptomes during the infection, further analyses were needed to evaluate changes in the expression levels of individual genes. A differential expression analysis was carried out to closely dissect the host murine macrophage response to infection by Leishmania major at 4, 24, 48, and 72 h after infection. Pairwise analyses were conducted within and across individual timepoints with infected macrophages evaluated against matched uninfected macrophages for each timepoint.
The most substantial response to infection by the macrophage was observed at 4 hpi, with 6897 genes differentially expressed (DE) between uninfected and infected cells with fold changes (FC) ranging from 29-fold downregulated to 56-fold upregulated. The response is reduced in later timepoints as reflected in smaller numbers of DE genes - 931, 1813, and 1460 genes at 24, 48, and 72 hpi, respectively - and in reduced fold changes with no downregulation beyond 12-fold or upregulation beyond 18-fold for these timepoints. The breakdown of up- and downregulated genes for each timepoint comparison is illustrated in Fig. 3a (top) and Additional file 2: Figure S3. Complete DE gene lists for each timepoint are provided in Additional file 3.
A direct comparison of the overlap in DE genes at each timepoint revealed that the macrophage response to infection at 4 hpi was vastly different from the response at the later timepoints (Fig. 3b and Additional file 2: Figure S4). Of the 6897 DE genes at 4 hpi, 93 % (6418 genes) were unique to 4 hpi. This is in contrast to 40 % of genes at 24 hpi (377), 33 % of genes at 48 hpi (598), and 33 % of genes at 72 hpi (484) that were DE only at that specific timepoint. Indeed, only 47 genes were consistently up- or downregulated at all 4 timepoints (38 genes upregulated at all timepoints and 9 genes downregulated at all timepoints; see Additional file 3). These genes do not appear to be functionally related, although the list does include the heavy metal transporters metallothionein 1 and 2. The latter has previously been associated with Leishmania resistance to treatment with antimonial drugs . Interestingly, two genes associated with Bcl2 (Bnip3 and Bcl2a1b), an inhibitor of apoptosis, are also on this list, suggesting that infection by L. major may induce changes in host gene expression to prevent macrophages from undergoing cell death.
In order to detect statistically significant differences in gene expression over time, we conducted differential expression analysis across timepoints. Contrasts between successive timepoints revealed a large number of DE genes during the 4 to 24 hpi transition (5674 DE genes), but almost no statistically significant genes between 24 and 48 hpi (1 gene) or between 48 and 72 hpi (0 genes) (Fig. 3a middle, Additional file 2: Figure S5B, and Additional file 3). This suggests that the large initial response of the murine macrophages to L. major infection stabilized by 24 hpi with no additional changes detected for the remainder of the time course as the infected cells mostly maintain their expression profile relative to matched uninfected cells.
The cellular processes most affected at each timepoint were characterized by KEGG pathway enrichment analysis using ConsensusPathDB . Genes that were DE >2-fold were used as input with upregulated and downregulated genes considered separately. The results of this analysis are reported in Table 1, Additional file 2: Table S1 and Table S2, and Additional file 4.
The enriched KEGG pathways were considered alongside the DE gene lists to gain insights into the cellular response to infection and how it differs between early and later infection. Early in the infection (4 hpi), many of the KEGG pathways that are most highly induced in infected macrophages are related to immune responses, specifically cytokine-cytokine receptor interactions and arginine and proline metabolism, glycolysis, and multiple signaling pathways including those for TNF, HIF-1, NF-kappa B, Jak-STAT, PI3K-Akt, and MAPK. A closer look at the DE genes in the immune response-related pathways reveals an interesting picture of the murine macrophage response to infection at 4 hpi. Infected macrophages produce a set of transcripts with paradoxical functions involved in both activating immune responses and promoting tissue regeneration and wound healing. Many of the DE genes associated with enriched pathways are anti-inflammatory in character or are involved in tissue growth and repair, including Csf1, Csf3, Il10, Il11r, Il1rn, Socs3, Hmox1, Egfr, Vegf, and fos-induced growth factor (Figf). The product of each of these genes has either been associated with reducing inflammation or promoting cell survival or differentiation. Therefore, macrophages infected with L. major appear to assume an immunoregulatory phenotype that resembles previously described macrophages treated with LPS and immune complex . However, not all of the differentially induced genes in infected macrophages at 4 hpi were anti-inflammatory. There was also an upregulation of transcripts encoding well-described inflammatory cytokines and their receptors including Il1, Il6, Tnf, and Nos2, as well as Il1rap and Il18r1. The pro-inflammatory NOD-like receptor signaling pathway is enriched among both upregulated and downregulated genes at 4 hpi with some members of the pathway upregulated (Tnf, Il6, Il1β, Cxcl2) and others downregulated (NLRP-family members, Pycard, Card9, and Il18). It is unclear to what extent the observed expression patterns are being driven by the parasite as it is taken into the host cell, or are being promoted by the host as it attempts to limit infection or to protect itself from the negative effects of a robust immune response.
Pathway analysis also showed an upregulation of glycolysis/gluconeogenesis at 4 hpi. Upregulated genes encode glycolytic enzymes, including hexokinases, enolase, phosphoglycerate kinase, glyceraldehyde-3- phosphate dehydrogenase, and lactate dehydrogenase A. An increase in anaerobic glycolysis has been associated with the stimulation of inflammatory responses in macrophages following TLR ligation . Glycolytic pathways were also enriched among upregulated genes at 24 hpi, but this effect was not seen at 48 and 72 hpi. Thus, following phagocytosis of L. major, macrophages undergo transient metabolic alterations that result in an increase in glycolysis and the generation of ATP.
KEGG pathways downregulated at 4 hpi are mostly related to lipid metabolism and biogenesis and include osteoclast differentiation, terpenoid backbone biosynthesis, steroid biosynthesis, PPAR signaling, and biosynthesis of unsaturated fatty acids. Decreases in the receptors and signaling molecules involved in the process of phagocytosis itself (“Fc gamma R-mediated phagocytosis” KEGG pathway) were also observed at 4 hpi. It has been previously shown that phagocytosis of Leishmania and other pathogens via Fc gamma receptors leads to parasite killing while phagocytosis via pathways mediated by complement receptor 3 (CR3) leads to parasite survival [61–65]. The downregulation of Fc gamma receptors may be a mechanism by which the parasite induces changes in host cells to promote alternative entry mechanisms that will support its survival. While we observed a slight downregulation of mRNA levels for C3ar1, C5ar1, and C5ar2 early in infection, upregulation of C5ar2 was detected at 48 and 72 hpi.
An examination of the enriched KEGG pathways at 24, 48, and 72 hpi reveals different enrichment patterns from those observed at 4 hpi (Additional file 2: Table S1 and Table S2 and Additional file 4). There are a greater number of pathways downregulated than upregulated during each of these later timepoints, but no clear picture otherwise emerges.
While a significant number of disease-specific KEGG pathways are enriched among up- and downregulated genes reported at various timepoints in this study (i.e., leishmaniasis, Chagas disease, tuberculosis, malaria, legionellosis, amoebiasis, Salmonella infection, Staphylococcus aureus infection, Fanconi anemia, prion diseases, and cancer subtypes; see Table 1 and Additional file 2: Table S1 and Table S2), a linkage between these disease and our dataset may not be particularly meaningful. This is because the KEGG pathways for these individual diseases appear to be constructed using a combination of heterogeneous observations from multiple studies of individual genes using data from different experiments, hosts, cell types, and timepoints, rather than using gene lists derived from a single or set of genome-wide studies, and as a result paint an incomplete picture of the global changes in gene expression. Therefore, while it is interesting to note small overlaps between the comprehensive DE profiles in this study and genes included in these pathways, the differences in scope between the lists limit the usefulness of comparisons to disease KEGG pathways as currently defined.
We compared the results of our differential expression analysis to previously published reports on the murine macrophage response to Leishmania infection. Since global transcriptome studies of L. major infection of macrophages isolated from C57BL/6 mice are not available, we compared our results to studies that used either the same mouse strain or the same Leishmania species (see Additional file 5). In a study detailing the infection of C57BL/6 peritoneal macrophages with L. amazonensis , 105 genes were identified as DE (FC > 1.5) between infected and uninfected macrophages at 18 hpi. This number compares to 3972 and 546 DE genes (FC > 1.5) at 4 and 24 hpi, respectively, in our dataset with an overlap of only 42 genes (of the 105) differentially expressed in the same direction at 4 or 24 hpi. Similarly, we compared the infection of Balb/c bone marrow-derived macrophages with L. major as reported in two recent studies. Of the 769 genes identified when comparing the transcriptomes of L. major-infected Balb/c macrophages to macrophages that had ingested latex beads (24 hpi timepoint), only 53 were in common with our 24 hpi timepoint . However, 446 genes were in common with our 4 hpi timepoint, indicating that the time progression observed between the two systems may be somewhat offset. Of the 40 genes identified as DE by more than 1.5-fold during the 24 h timecourse experiment reported by Rabhi et al. , 29 genes were also found to be differentially expressed in the same direction in our dataset at either 4 or 24 hpi with the greatest number of genes found in common with our 4 hpi dataset (26 genes). The limited agreement between the previously reported DE genes and our dataset is likely a reflection of the differences in study design (e.g., combination of mouse strain and parasite species, macrophage source, parasite opsonization), batch effects between laboratories, technical platforms, and data analysis methods.
Identification of differentially expressed genes and gene ontology-based enrichment analyses in L. major parasites
While past studies have used microarrays and other methods to study changes in Leishmania gene expression between promastigotes and amastigotes [27, 28, 31, 35, 37, 38, 44, 47, 66, 67], this study provided a unique opportunity to additionally characterize gene expression patterns within amastigotes over the course of an intracellular infection. Differential expression analysis revealed a large number of genes (2962) to be DE between metacyclic promastigotes and 4-h amastigotes at an adjusted P value cutoff of <0.05, reflecting the significant differentiation event that occurs as the parasite enters host cells. Significantly fewer changes in expression were observed in amastigotes across timepoints, with 301 DE genes found between 4- and 24-h amastigotes and no DE genes between either 24- and 48-h amastigotes or between 48- and 72-h amastigotes (Fig. 3a, bottom, Additional file 2: Figure S5A, and Additional file 6). This pattern of differential gene expression is consistent with observations of gene expression over time in the mouse macrophages and suggests that the initial reprogramming that occurs upon infection had largely been stabilized by 24 hpi.
Gene ontology (GO) analysis was used to identify cellular functions and processes that characterize the entry and survival of L. major in the murine macrophage system. These results were considered alongside the lists of DE genes to draw insights into how the parasite adapts to its new environment. Efforts were made to compare the results of this study to previous work, though our ability to conduct meaningful systematic comparisons was constrained due to the use of different parasite systems (e.g., Leishmania species, source of amastigotes, developmental stages and timepoints collected), host systems (e.g., mice strains, cell lines), experimental platforms (e.g., microarrays, number of genes interrogated), and methods for assessing and reporting differential expression (fold change vs. statistical cut-offs; details of reported gene lists).
GO analysis of the up- and downregulated genes identified during the metacyclic promastigote to 4-h amastigote transition revealed a total of 20 enriched GO categories. Enriched GO terms for genes upregulated during this transition were primarily related to mitigating the effects of an oxidative response by the host immune system and the regulation of proteins (Table 2 and Additional file 7). Heat shock proteins, especially HSP83, multiple tryparedoxin peroxidase family members, and multiple cyclophilins were upregulated upon entry of metacyclic promastigotes into host cells and contributed strongly to the GO enrichment results for the upregulated genes. HSP83, the cytoplasmic form of HSP90, is known to stabilize transcription factors and kinases and is thus largely involved in the regulation of signaling networks involved in differentiation [68, 69], and tryparedoxin peroxidase is known to contribute to Leishmania virulence and resistance to antimonial drugs . Less is known about the role of cyclophilins in Leishmania, though these peptidyl-prolyl cis/trans isomerases have been shown to reactivate cellular proteins by promoting their disaggregation . The upregulation of cyclophilin in L. major amastigotes has been previously observed (compared to procyclic promastigotes) . Many types of surface antigens were highly upregulated, as was the zinc-metalloprotease GP63 (a known virulence factor that subverts macrophage signaling) , and phosphoglycan beta 1,3 galactosyltransferase (SCG) family members, which are responsible for modifications to the Leishmania lipophosphoglycan (LPG) surface molecule side chains  and have been previously observed to be upregulated in Leishmania amastigotes . That there are significant changes taking place on the surface of the parasite is not surprising given the transformation that takes place as differentiation progresses from the promastigote to amastigote forms. Other genes upregulated during the metacyclic promastigote to 4-h amastigote transition included meta1, which is thought to play a role in virulence and was previously found to be upregulated in infective metacyclic promastigotes relative to non-infective procyclic promastigotes [55, 74, 75], and inositol-3-phosphate synthase (ino1), which synthesizes myo-inositol, a precursor molecule for the backbone of the GPI anchors used by many Leishmania surface proteins, including multiple virulence factors . Interestingly, the gene encoding the transporter for myo-inositol (mit1)  was among the most downregulated during the transition, perhaps indicating that the parasite favors synthesizing de novo rather than importing this important substrate once it enters mammalian cells.
GO terms that were enriched among downregulated genes during the metacyclic promastigote to 4-h amastigote transition were related to translation, cell signaling, fatty acid biosynthesis, and flagellum structure (Table 2 and Additional file 7). A number of genes were responsible for driving the GO results, including ribosomal proteins, casein kinase, receptor-type adenylate cyclase, fatty acid elongase, sphingolipid delta desaturase, and paraflagellar rod protein. The downregulation of ribosomal proteins during the Leishmania promastigote to amastigote transition is consistent with previous reports  and may reflect a reduction in translation taking place within the cell. Casein kinase is thought to play a role in Leishmania promastigote growth and parasite virulence via interactions with host proteins, though its exact function is unclear [78, 79]. Likewise, adenylate cyclase is a suspected regulator of cAMP signaling during Leishmania differentiation [80, 81], but the exactly signaling pathways through which this is accomplished are unknown. Fatty acids have many potential roles within the cell related to membrane structure and composition, signaling, and energy. Fatty acid elongase family members are involved in fatty acid synthesis and may be involved in GPI anchoring , while sphingolipid delta 4 desaturase is involved in the degradation of sphingolipids, a process that has been linked to enabling Leishmania to survive the acidic environment of the phagolysosome once it is taken into host cells . Genes that influence microtubule length and dynamics were systematically downregulated, including those that encode paraflagellar rod proteins, calmodulin-related proteins , multiple kinesins [85, 86], NIMA-related kinases , and mitogen-activated protein (MAP) kinase kinases [88–90]. These genes have been implicated in the regulation of flagellum length and mitotic spindle formation and may reflect changes in morphology and cell division that take place as the flagellated promastigote stage parasites transform into the aflagellated amastigotes. Consistent with previous microarray-based studies [27, 28], we also detected the downregulation of the known metacyclic markers SHERP and HASP as the parasite transformed from promastigotes to amastigotes.
Of the 14 GO categories enriched during the transition from 4 to 24-h amastigotes, those related to microtubules and protein complexes were specific to upregulated genes with many copies of β-tubulin, the primary component of microtubules, almost exclusively driving the GO analysis results. The upregulation of β-tubulin in early amastigotes compared to metacyclic promastigotes is noteworthy considering its downregulation in metacyclic promastigotes compared to procyclic promastigotes [55, 91]. Besides β-tubulin, many other top upregulated genes encoded surface antigens, including the developmentally-regulated amastigote-enriched protein amastin [27, 35, 44, 92]. GO categories related to glucose metabolism and oxidation-reduction processes were enriched among downregulated genes with many of the downregulated genes contributing to these categories belonging to the glycolysis pathway (GAPDH, pyruvate kinase, ENOL, triosephosphate isomerase, ALD, and PGI) or related to cholesterol/sterol metabolism (NSDHL, HMGR, FPPS, CYP51) or purine metabolism (XRPT). This may indicate a change in the metabolic requirements or preferences of Leishmania once the parasites have survived their initial entry into the mammalian cell. Other highly downregulated genes included heat shock proteins, biopterin transporters, and meta1 [55, 74, 75].
A very significant dataset resulting from this study includes the large proportion of the DE genes which have no known function and are annotated as hypothetical proteins. Those make up 58 % (1724 of 2962) of the DE genes in the metacyclic to 4-h amastigote transition and 49 % (147 of 301) of the DE genes in the 4- to 24-h amastigote transition. While hypothetical genes have largely been overlooked to date, they almost certainly constitute an integral part of the transcriptomic signature of the parasites as members of co-expressed gene networks which are involved in common functional pathways or regulated by shared regulatory mechanisms.
In this study, we performed transcriptomic profiling to identify genes that were differentially expressed in both parasite and host cells as L. major entered and persisted within murine macrophages during the first 72 h of an infection. The generation of RNA-seq data followed by the unambiguous mapping of reads from infected samples to the genomes of both the host and the parasite enabled us to identify genes that changed over the course of infection in the real-time context of a dynamic dual biological system, and to a much greater depth and sensitivity than has been previously reported. Collection of data from multiple biological replicates, the use of matched host control samples, careful statistical analysis of variation, and removal of batch effects provided us with a unique ability to detect biological differences between samples and timepoints with high confidence and sensitivity.
In addition to providing robust sets of markers for multiple developmental stages of L. major parasites and Leishmania-infected macrophages over several timepoints, this work contributes to a growing body of literature on the broader field of host-pathogen interactions. Indeed, the comprehensive dataset generated in this study will also serve as a reference for future studies using different Leishmania strains (or even different pathogens altogether) to examine infection of macrophages from multiple sources and in various states of activation, polarization, or rest. In conjunction with datasets that will be produced for other pathogens, a clearer picture of the signature of intracellular infection will emerge, providing additional insights into how pathogens are able to evade host defenses and modulate the biological functions of the cell in order to survive in the mammalian environment.
Leishmania major (clone V1, MHOM/IL/80/Friedlin) was isolated after passage through Balb/c mice. Promastigotes were grown in 50 % M199, 39 % Schneider’s medium along with 10 % FBS and 100 U/mL of Penicillin-Streptomycin at 25 °C. L. major promastigotes were not split for more than 5 passages to maintain virulence of the cultures. Enrichment for metacyclic promastigotes from stationary phase cultures was done by Ficoll density gradient centrifugation .
All procedures involving animals were approved by the Institutional Animal Care and Use Committee (IACUC) of the University of Maryland. Peritoneal macrophages were isolated from C57BL/6 mice (7 weeks, female) obtained from National Cancer Institute Charles River Laboratories by flushing the peritoneal cavity using 10 mL cold DPBS without calcium and magnesium. Macrophages isolated from ten mice were pooled for each experiment. Cells were plated in six well plates in 2.5 mL DMEM/F12 media supplemented with 10 % FBS, 1 % penicillin/streptomycin, and 1 % glutamine to a density of ~2.5 x 106 cells/well and incubated overnight. Two hours after washing, cells were infected with Ficoll-enriched L. major metacyclic promastigotes at a ratio of 5 parasites per macrophage along with 5 % C5-deficient serum collected from DBA2 mice. Cells were lysed using the Trizol® reagent (Invitrogen, CA) at 4, 24, 48, and 72 h following infection.
RNA isolation and cDNA library preparation
Total RNA was isolated using the Trizol® reagent (Invitrogen, CA), treated with DNase, and purified using the Qiagen RNeasy mini kit. RNA integrity was assessed using an Agilent 2100 bioanalyzer. Poly (A)-enriched cDNA libraries were generated using the Illumina TruSeq Sample Preparation kit (San Diego, CA) and checked for quality and quantity using the bioanalyzer and qPCR (KAPA Biosystems).
RNA-seq data generation, pre-processing, and quality trimming
Paired end reads (100 bp) were obtained using the Illumina HiSeq 1500 platform. Trimmomatic [ 94] was used to remove any remaining Illumina adapter sequences from reads and to trim bases off the start or the end of a read when the quality score fell below a threshold of 20. Sequence quality metrics were assessed using FastQC [http://www.bioinformatics.babraham.ac.uk/projects/fastqc/].
Mapping cDNA fragments to the reference genome, abundance estimation, and data normalization
Reads were aligned independently to the L. major genome (v. 6.0) obtained from the TriTrypDB database (www.tritrypdb.org) and to the mouse genome (v. mm10) from the UCSC genome browser (http://genome.ucsc.edu) using TopHat (v 2.0.10) . Gene model annotations were provided for the mapping (option –G) with limitations on the identification of novel splice junctions (option –no-novel-juncs). Two mismatches per read were allowed (default parameter) and reads were allowed to map only to a single locus (option –g 1). The abundance of reads mapping to each gene feature in the aligned genome was determined using HTSeq[ 96]. Each resulting count table was restricted to protein-coding genes (23,100 genes for mouse and 8486 genes for L. major).
Data quality assessment by statistical sample clustering and visualization
Multiple approaches were used to evaluate replicates and to visualize sample-sample distances. Those included Pearson correlation, median pairwise correlation analysis, box plots, Principal Component Analysis (PCA) and Euclidean distances-based hierarchical clustering. All components of the statistical pipeline, which we named cbcbSEQ, were done in R and can be accessed on GitHub (https://github.com/kokrah/cbcbSEQ/).
Differential expression analysis
Non-expressed and weakly expressed genes, defined as having less than 1 read per million in n of the samples, where n is the size of the smallest group of replicates  (here n = 3), were removed prior to differential expression analysis. This is implemented by the filterCounts function within the cbcbSEQ R package (see above). A quantile normalization scheme was applied to all samples . Following log2 transformation of the data, limma (a Bioconductor package) was used to conduct differential expression analyses. limma utilizes a standard variance moderated across all genes using a Bayesian model and produces P values with greater degrees of freedom . The voom module was used to transform the data based on observational level weights derived from the mean-variance relationship prior to statistical modeling  (Additional file 2: Figure S6). Experimental batch effects were adjusted for by including experimental batch as a covariate in the statistical model . Differentially expressed (DE) genes were defined as genes with a Benjamini-Hochberg multiple-testing adjusted P value of < 0.05.
KEGG pathway analysis
KEGG pathway analysis using ConsensusPathDB-mouse was done to identify signaling and metabolic pathways that were over-represented in the mouse DE gene lists. For each KEGG pathway, a P value was calculated using a hypergeometric test and a cutoff of 0.01 was applied to identify enriched KEGG pathways. Genes that were DE more than 2-fold in L. major-infected cells relative to uninfected controls at each timepoint were used as input with up- and downregulated genes considered separately.
Gene ontology (GO) analysis
GO categories enriched in the L. major DE gene lists were identified using the GOseq package in R . GOseq was developed specifically to account for transcript length bias in GO analyses using RNA-seq data. For each comparison, upregulated and downregulated gene sets (no fold change cut-off) were input separately into GOseq. A P value cut-off of 0.05 was used.
Availability of supporting data
The data sets supporting the results of this article are available in the NCBI Sequence Read Archive (SRA) repository [SRA: SRR1460724-SRR1460747, SRR1460767, SRR1460772, SRR2136702].
principal component analysis
Sequence Read Archive
Arango Duque G, Descoteaux A. Leishmania survival in the macrophage: where the ends justify the means. Curr Opin Microbiol. 2015;26:32–40.
Flávia Nardy A, Freire-de-Lima CG, Morrot A. Immune Evasion Strategies of Trypanosoma cruzi. J Immunol Res. 2015;2015:178947.
Jones BD, Faron M, Rasmussen JA, Fletcher JR. Uncovering the components of the Francisella tularensis virulence stealth strategy. Front Cell Infect Microbiol. 2014;4:32.
Pittman KJ, Knoll LJ. Long-Term Relationships: the Complicated Interplay between the Host and the Developmental Stages of Toxoplasma gondii during Acute and Chronic Infections. Microbiol Mol Biol Rev. 2015;79:387–401.
Rikihisa Y. Ehrlichia subversion of host innate responses. Curr Opin Microbiol. 2006;9:95–101.
Sia JK, Georgieva M, Rengarajan J. Innate Immune Defenses in Human Tuberculosis: An Overview of the Interactions between Mycobacterium tuberculosis and Innate Immune Cells. J Immunol Res. 2015;2015:747543.
Simon S, Hilbi H. Subversion of Cell-Autonomous Immunity and Cell Migration by Legionella pneumophila Effectors. Front Immunol. 2015;6:447.
Alvar J, Vélez ID, Bern C, Herrero M, Desjeux P, Cano J, et al. Leishmaniasis worldwide and global estimates of its incidence. PLoS One. 2012;7, e35671.
Bogdan C, Donhauser N, Döring R, Röllinghoff M, Diefenbach A, Rittig MG. Fibroblasts as host cells in latent leishmaniosis. J Exp Med. 2000;191:2121–30.
Laufs H, Müller K, Fleischer J, Reiling N, Jahnke N, Jensenius JC, et al. Intracellular survival of Leishmania major in neutrophil granulocytes after uptake in the absence of heat-labile serum factors. Infect Immun. 2002;70:826–35.
Moll H, Flohé S, Röllinghoff M. Dendritic cells in Leishmania major-immune mice harbor persistent parasites and mediate an antigen-specific T cell immune response. Eur J Immunol. 1995;25:693–9.
Peters NC, Egen JG, Secundino N, Debrabant A, Kimblin N, Kamhawi S, et al. In vivo imaging reveals an essential role for neutrophils in leishmaniasis transmitted by sand flies. Science. 2008;321:970–4.
Sarkar A, Aga E, Bussmeyer U, Bhattacharyya A, Möller S, Hellberg L, et al. Infection of neutrophil granulocytes with Leishmania major activates ERK 1/2 and modulates multiple apoptotic pathways to inhibit apoptosis. Med Microbiol Immunol. 2013;202:25–35.
Ambit A, Woods KL, Cull B, Coombs GH, Mottram JC. Morphological events during the cell cycle of Leishmania major. Eukaryot Cell. 2011;10:1429–38.
Beverley SM, Turco SJ. Lipophosphoglycan (LPG) and the identification of virulence genes in the protozoan parasite Leishmania. Trends Microbiol. 1998;6:35–40.
Wheeler RJ, Gluenz E, Gull K. The cell cycle of Leishmania: morphogenetic events and their implications for parasite biology. Mol Microbiol. 2011;79:647–62.
Etges R, Müller I. Progressive disease or protective immunity to Leishmania major infection: the result of a network of stimulatory and inhibitory interactions. J Mol Med. 1998;76:372–90.
Reiner SL, Locksley RM. The regulation of immunity to Leishmania major. Annu Rev Immunol. 1995;13:151–77.
Scharton-Kersten T, Scott P. The role of the innate immune response in Th1 cell development following Leishmania major infection. J Leukoc Biol. 1995;57:515–22.
Bennett CL, Misslitz A, Colledge L, Aebischer T, Blackburn CC. Silent infection of bone marrow-derived dendritic cells by Leishmania mexicana amastigotes. Eur J Immunol. 2001;31:876–83.
Laskay T, van Zandbergen G, Solbach W. Neutrophil granulocytes--Trojan horses for Leishmania major and other intracellular microbes? Trends Microbiol. 2003;11:210–4.
Locksley RM, Heinzel FP, Fankhauser JE, Nelson CS, Sadick MD. Cutaneous host defense in leishmaniasis: interaction of isolated dermal macrophages and epidermal Langerhans cells with the insect-stage promastigote. Infect Immun. 1988;56:336–42.
Zhang S, Kim CC, Batra S, McKerrow JH, Loke P. Delineation of diverse macrophage activation programs in response to intracellular parasites and cytokines. PLoS Negl Trop Dis. 2010;4, e648.
Kaye P, Scott P. Leishmaniasis: complexity at the host-pathogen interface. Nat Rev Microbiol. 2011;9:604–15.
Olivier M, Gregory DJ, Forget G. Subversion mechanisms by which Leishmania parasites can escape the host immune response: a signaling point of view. Clin Microbiol Rev. 2005;18:293–305.
Sacks D, Sher A. Evasion of innate immunity by parasitic protozoa. Nat Immunol. 2002;3:1041–7.
Akopyants NS, Matlib RS, Bukanova EN, Smeds MR, Brownstein BH, Stormo GD, et al. Expression profiling using random genomic DNA microarrays identifies differentially expressed genes associated with three major developmental stages of the protozoan parasite Leishmania major. Mol Biochem Parasitol. 2004;136:71–86.
Almeida R, Gilmartin BJ, McCann SH, Norrish A, Ivens AC, Lawson D, et al. Expression profiling of the Leishmania life cycle: cDNA arrays identify developmentally regulated genes present but not annotated in the genome. Mol Biochem Parasitol. 2004;136:87–100.
Buates S, Matlashewski G. General suppression of macrophage gene expression during Leishmania donovani infection. J Immunol. 2001;166:3416–22.
Chaussabel D, Semnani RT, McDowell MA, Sacks D, Sher A, Nutman TB. Unique gene expression profiles of human macrophages and dendritic cells to phylogenetically distinct parasites. Blood. 2003;102:672–81.
Depledge DP, Evans KJ, Ivens AC, Aziz N, Maroof A, Kaye PM, et al. Comparative expression profiling of Leishmania: modulation in gene expression between species and in different host genetic backgrounds. PLoS Negl Trop Dis. 2009;3, e476.
Ettinger NA, Wilson ME. Macrophage and T-cell gene expression in a model of early infection with the protozoan Leishmania chagasi. PLoS Negl Trop Dis. 2008;2, e252.
Giraud E, Lecoeur H, Soubigou G, Coppée J-Y, Milon G, Prina E, et al. Distinct transcriptional signatures of bone marrow-derived C57BL/6 and DBA/2 dendritic leucocytes hosting live Leishmania amazonensis amastigotes. PLoS Negl Trop Dis. 2012;6, e1980.
Gregory DJ, Sladek R, Olivier M, Matlashewski G. Comparison of the effects of Leishmania major or Leishmania donovani infection on macrophage gene expression. Infect Immun. 2008;76:1186–92.
Holzer TR, McMaster WR, Forney JD. Expression profiling by whole-genome interspecies microarray hybridization reveals differential gene expression in procyclic promastigotes, lesion-derived amastigotes, and axenic amastigotes in Leishmania mexicana. Mol Biochem Parasitol. 2006;146:198–218.
Lahav T, Sivam D, Volpin H, Ronen M, Tsigankov P, Green A, et al. Multiple levels of gene regulation mediate differentiation of the intracellular pathogen Leishmania. FASEB J. 2011;25:515–25.
Leifso K, Cohen-Freue G, Dogra N, Murray A, McMaster WR. Genomic and proteomic expression analysis of Leishmania promastigote and amastigote life stages: the Leishmania genome is constitutively expressed. Mol Biochem Parasitol. 2007;152:35–46.
McNicoll F, Drummelsmith J, Müller M, Madore E, Boilard N, Ouellette M, et al. A combined proteomic and transcriptomic approach to the study of stage differentiation in Leishmania infantum. Proteomics. 2006;6:3567–81.
Novais FO, Carvalho LP, Passos S, Roos DS, Carvalho EM, Scott P, et al. Genomic profiling of human Leishmania braziliensis lesions identifies transcriptional modules associated with cutaneous immunopathology. J Invest Dermatol. 2015;135:94–101.
Osorio y Fortéa J, de La Llave E, Regnault B, Coppée J-Y, Milon G, Lang T, et al. Transcriptional signatures of BALB/c mouse macrophages housing multiplying Leishmania amazonensis amastigotes. BMC Genomics. 2009;10:119.
Probst CM, Silva RA, Menezes JPB, Almeida TF, Gomes IN, Dallabona AC, et al. A comparison of two distinct murine macrophage gene expression profiles in response to Leishmania amazonensis infection. BMC Microbiol. 2012;12:22.
Rabhi I, Rabhi S, Ben-Othman R, Rasche A, Daskalaki A, Trentin B, et al. Transcriptomic signature of Leishmania infected mice macrophages: a metabolic point of view. PLoS Negl Trop Dis. 2012;6, e1763.
Ramírez C, Díaz-Toro Y, Tellez J, Castilho TM, Rojas R, Ettinger NA, et al. Human macrophage response to L. (Viannia) panamensis: microarray evidence for an early inflammatory response. PLoS Negl Trop Dis. 2012;6, e1866.
Rochette A, Raymond F, Ubeda J-M, Smith M, Messier N, Boisvert S, et al. Genome-wide gene expression profiling analysis of Leishmania major and Leishmania infantum developmental stages reveals substantial differences between the two species. BMC Genomics. 2008;9:255.
Rochette A, Raymond F, Corbeil J, Ouellette M, Papadopoulou B. Whole-genome comparative RNA expression profiling of axenic and intracellular amastigote forms of Leishmania infantum. Mol Biochem Parasitol. 2009;165:32–47.
Rodriguez NE, Chang HK, Wilson ME. Novel program of macrophage gene expression induced by phagocytosis of Leishmania chagasi. Infect Immun. 2004;72:2111–22.
Saxena A, Lahav T, Holland N, Aggarwal G, Anupama A, Huang Y, et al. Analysis of the Leishmania donovani transcriptome reveals an ordered progression of transient and permanent changes in gene expression during differentiation. Mol Biochem Parasitol. 2007;152:53–65.
Maretti-Mira AC, Bittner J, Oliveira-Neto MP, Liu M, Kang D, Li H, et al. Transcriptome patterns from primary cutaneous Leishmania braziliensis infections associate with eventual development of mucosal disease in humans. PLoS Negl Trop Dis. 2012;6, e1816.
Ivens AC, Peacock CS, Worthey EA, Murphy L, Aggarwal G, Berriman M, et al. The genome of the kinetoplastid parasite, Leishmania major. Science. 2005;309:436–42.
Martinez-Calvillo S, Yan S, Nguyen D, Fox M, Stuart K, Myler PJ. Transcription of Leishmania major Friedlin Chromosome 1 Initiates in Both Directions within a Single Region. Mol Cell. 2003;11:1291–9.
LeBowitz JH, Smith HQ, Rusche L, Beverley SM. Coupling of poly (A) site selection and trans-splicing in Leishmania. Genes Dev. 1993;7:996–1007.
Matthews KR, Tschudi C, Ullu E. A common pyrimidine-rich motif governs trans-splicing and polyadenylation of tubulin polycistronic pre-mRNA in trypanosomes. Genes Dev. 1994;8:491–501.
Ullu E, Matthews KR, Tschudi C. Temporal order of RNA-processing reactions in trypanosomes: rapid trans splicing precedes polyadenylation of newly synthesized tubulin transcripts. Mol Cell Biol. 1993;13:720–5.
Clayton C, Shapira M. Post-transcriptional regulation of gene expression in trypanosomes and leishmanias. Mol Biochem Parasitol. 2007;156:93–101.
Dillon LAL, Okrah K, Hughitt VK, Suresh R, Li Y, Fernandes MC, et al. Transcriptomic profiling of gene expression and RNA processing during Leishmania major differentiation. Nucleic Acids Res. 2015;43:6799-813.
Leek JT, Scharpf RB, Bravo HC, Simcha D, Langmead B, Johnson WE, et al. Tackling the widespread and critical impact of batch effects in high-throughput data. Nat Rev Genet. 2010;11:733–9.
Gómez MA, Navas A, Márquez R, Rojas LJ, Vargas DA, Blanco VM, et al. Leishmania panamensis infection and antimonial drugs modulate expression of macrophage drug transporters and metabolizing enzymes: impact on intracellular parasite survival. J Antimicrob Chemother. 2014;69:139–49.
Kamburov A, Pentchev K, Galicka H, Wierling C, Lehrach H, Herwig R. ConsensusPathDB: toward a more complete picture of cell biology. Nucleic Acids Res. 2011;39:D712–7.
Fleming BD, Chandrasekaran P, Dillon LAL, Dalby E, Suresh R, Sarkar A, et al.. The generation of macrophages with anti-inflammatory activity in the absence of STAT6 signaling. J Leukoc Biol. 2015;98:395-407.
Tannahill GM, Curtis AM, Adamik J, Palsson-McDermott EM, McGettrick AF, Goel G, et al. Succinate is an inflammatory signal that induces IL-1β through HIF-1α. Nature. 2013;496:238–42.
Babadjanova Z, Wiedinger K, Gosselin EJ, Bitsaktsis C. Targeting of a Fixed Bacterial Immunogen to Fc Receptors Reverses the Anti-Inflammatory Properties of the Gram-Negative Bacterium, Francisella tularensis, during the Early Stages of Infection. PLoS One. 2015;10, e0129981.
Da Silva RP, Hall BF, Joiner KA, Sacks DL. CR1, the C3b receptor, mediates binding of infective Leishmania major metacyclic promastigotes to human macrophages. J Immunol. 1989;143:617–22.
Fadul CE, Channon JY, Kasper LH. Survival of immunoglobulin G-opsonized Toxoplasma gondii in nonadherent human monocytes. Infect Immun. 1995;63:4290–4.
Mosser DM, Edelson PJ. The third component of complement (C3) is responsible for the intracellular survival of Leishmania major. Nature. 1987;327:329–31.
Woelbing F, Kostka SL, Moelle K, Belkaid Y, Sunderkoetter C, Verbeek S, et al. Uptake of Leishmania major by dendritic cells is mediated by Fcgamma receptors and facilitates acquisition of protective immunity. J Exp Med. 2006;203:177–88.
Quijada L, Soto M, Requena JM. Genomic DNA macroarrays as a tool for analysis of gene expression in Leishmania. Exp Parasitol. 2005;111:64–70.
Walker J, Vasquez J-J, Gomez MA, Drummelsmith J, Burchmore R, Girard I, et al. Identification of developmentally-regulated proteins in Leishmania panamensis by proteome profiling of promastigotes and axenic amastigotes. Mol Biochem Parasitol. 2006;147:64–73.
Hombach A, Clos J. No stress--Hsp90 and signal transduction in Leishmania. Parasitology. 2014;141:1156–66.
Requena JM, Montalvo AM, Fraga J. Molecular Chaperones of Leishmania: Central Players in Many Stress-Related and -Unrelated Physiological Processes. Biomed Res Int. 2015;2015:301326.
Iyer JP, Kaprakkaden A, Choudhary ML, Shaha C. Crucial role of cytosolic tryparedoxin peroxidase in Leishmania donovani survival, drug response and virulence. Mol Microbiol. 2008;68:372–91.
Chakraborty A, Das I, Datta R, Sen B, Bhattacharyya D, Mandal C, et al. A single-domain cyclophilin from Leishmania donovani reactivates soluble aggregates of adenosine kinase by isomerase-independent chaperone function. J Biol Chem. 2002;277:47451–60.
Olivier M, Atayde VD, Isnard A, Hassani K, Shio MT. Leishmania virulence factors: focus on the metalloprotease GP63. Microbes Infect. 2012;14:1377–89.
Dobson DE, Scholtes LD, Valdez KE, Sullivan DR, Mengeling BJ, Cilmi S, et al. Functional identification of galactosyltransferases (SCGs) required for species-specific modifications of the lipophosphoglycan adhesin controlling Leishmania major-sand fly interactions. J Biol Chem. 2003;278:15523–31.
Nourbakhsh F, Uliana SR, Smith DF. Characterisation and expression of a stage-regulated gene of Leishmania major. Mol Biochem Parasitol. 1996;76:201–13.
Puri V, Goyal A, Sankaranarayanan R, Enright AJ, Vaidya T. Evolutionary and functional insights into Leishmania META1: evidence for lateral gene transfer and a role for META1 in secretion. BMC Evol Biol. 2011;11:334.
Ilg T. Generation of myo-inositol-auxotrophic Leishmania mexicana mutants by targeted replacement of the myo-inositol-1-phosphate synthase gene. Mol Biochem Parasitol. 2002;120:151–6.
Drew ME, Langford CK, Klamo EM, Russell DG, Kavanaugh MP, Landfear SM. Functional expression of a myo-inositol/H+ symporter from Leishmania donovani. Mol Cell Biol. 1995;15:5508–15.
Allocco JJ, Donald R, Zhong T, Lee A, Tang YS, Hendrickson RC, et al. Inhibitors of casein kinase 1 block the growth of Leishmania major promastigotes in vitro. Int J Parasitol. 2006;36:1249–59.
Dan-Goor M, Nasereddin A, Jaber H, Jaffe CL. Identification of a secreted casein kinase 1 in Leishmania donovani: effect of protein over expression on parasite growth and virulence. PLoS One. 2013;8, e79287.
Biswas A, Bhattacharya A, Das PK. Role of cAMP Signaling in the Survival and Infectivity of the Protozoan Parasite, Leishmania donovani. Mol Biol Int. 2011;2011:782971.
Sanchez MA, Zeoli D, Klamo EM, Kavanaugh MP, Landfear SM. A family of putative receptor-adenylate cyclases from Leishmania donovani. J Biol Chem. 1995;270:17551–8.
Ramakrishnan S, Serricchio M, Striepen B, Bütikofer P. Lipid synthesis in protozoan parasites: a comparison between kinetoplastids and apicomplexans. Prog Lipid Res. 2013;52:488–512.
Xu W, Xin L, Soong L, Zhang K. Sphingolipid degradation by Leishmania major is required for its resistance to acidic pH in the mammalian host. Infect Immun. 2011;79:3377–87.
Ginger ML, Collingridge PW, Brown RWB, Sproat R, Shaw MK, Gull K. Calmodulin is required for paraflagellar rod assembly and flagellum-cell body attachment in trypanosomes. Protist. 2013;164:528–40.
Blaineau C, Tessier M, Dubessay P, Tasse L, Crobu L, Pagès M, et al. A novel microtubule-depolymerizing kinesin involved in length control of a eukaryotic flagellum. Curr Biol. 2007;17:778–82.
Vicente JJ, Wordeman L. Mitosis, microtubule dynamics and the evolution of kinesins. Exp Cell Res. 2015;334:61–9.
Cao M, Li G, Pan J. Regulation of cilia assembly, disassembly, and length by protein phosphorylation. Methods Cell Biol. 2009;94:333–46.
Bengs F, Scholz A, Kuhn D, Wiese M. LmxMPK9, a mitogen-activated protein kinase homologue affects flagellar length in Leishmania mexicana. Mol Microbiol. 2005;55:1606–15.
Erdmann M, Scholz A, Melzer IM, Schmetz C, Wiese M. Interacting protein kinases involved in the regulation of flagellar length. Mol Biol Cell. 2006;17:2035–45.
Wiese M, Kuhn D, Grünfelder CG. Protein kinase involved in flagellar-length control. Eukaryot Cell. 2003;2:769–77.
Coulson RM, Connor V, Chen JC, Ajioka JW. Differential expression of Leishmania major beta-tubulin genes during the acquisition of promastigote infectivity. Mol Biochem Parasitol. 1996;82:227–36.
Wu Y, El Fakhry Y, Sereno D, Tamar S, Papadopoulou B. A new developmentally regulated gene family in Leishmania amastigotes encoding a homolog of amastin surface proteins. Mol Biochem Parasitol. 2000;110:345–57.
Späth GF, Beverley SM. A lipophosphoglycan-independent method for isolation of infective Leishmania metacyclic promastigotes by density gradient centrifugation. Exp Parasitol. 2001;99:97–103.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25:1105–11.
Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics. 2014;31:166-9.
Anders S, McCarthy DJ, Chen Y, Okoniewski M, Smyth GK, Huber W, et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nat Protoc. 2013;8:1765–86.
Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003;19:185–93.
Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3.
Law CW, Chen Y, Shi W, Smyth GK. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29.
Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.
All the sequencing was performed at the University of Maryland Institute for Bioscience and Biotechnology Research (IBBR) sequencing core. This work was supported by a grant from NIH (AI094773) to NES and DMM.
The authors declare that they have no competing interests.
L.A.L.D., D.M.M. and N.E.S. designed research; L.A.L.D. and R.S. performed experiments; L.A.L.D. and N.E.S. analyzed data; K.O. and H.C.B. contributed statistical and analytic tools; L.A.L.D., D.M.M., and N.E.S. wrote the paper; N.E.S. and D.M.M. conceived the project. All authors have read and approved the manuscript.
Experimental design. (XLSX 14 kb)
Figure S1. Global assessment of data distribution. Figure S2. Principal Component Analysis (PCA) and hierarchical clustering analysis prior to accounting for batch effects. Figure S3. Global changes in murine macrophage gene expression upon infection by L. major. Figure S4. Expression patterns for the top 20 up- and down-regulated genes in L. major-infected murine macrophages at 4, 24, 48, and 72 hpi. Figure S5. Global changes in L. major and murine macrophage gene expression over the course of infection. Figure S6. Mean-variance curve modeling and fitting of a local regression trend line by voom. Table S1. KEGG pathways enriched among genes upregulated in L. major-infected murine macrophages at 24, 48, and 72 hpi. Table S2. KEGG pathways enriched among genes downregulated in L. major-infected murine macrophages at 24, 48, and 72 hpi. (PDF 25627 kb)
Differentially expressed gene lists for uninfected vs. L. major -infected murine macrophages within and across time using limma. (XLS 2631 kb)
Enriched KEGG pathways and corresponding murine differentially expressed genes. (XLS 228 kb)
Comparison of differentially expressed murine gene lists with published literature. (XLS 267 kb)
Differentially expressed gene lists for L. major parasite stages and timepoints using limma. (XLS 547 kb)
Enriched GO categories and corresponding L. major differentially expressed genes. (XLS 178 kb)