Leishmania infection induces a limited differential gene expression in the sand fly midgut

Background Sand flies are the vectors of Leishmania parasites. To develop in the sand fly midgut, Leishmania multiplies and undergoes various stage differentiations giving rise to the infective form, the metacyclic promastigotes. To determine the changes in sand fly midgut gene expression caused by the presence of Leishmania, we performed RNA-Seq of uninfected and Leishmania infantum-infected Lutzomyia longipalpis midguts from seven different libraries corresponding to time points which cover the various Leishmania developmental stages. Results The combined transcriptomes resulted in the de novo assembly of 13,841 sand fly midgut transcripts. Importantly, only 113 sand fly transcripts, about 1%, were differentially expressed in the presence of Leishmania parasites. Further, we observed distinct differentially expressed sand fly midgut transcripts corresponding to the presence of each of the various Leishmania stages suggesting that each parasite stage influences midgut gene expression in a specific manner. Two main patterns of sand fly gene expression modulation were noted. At early time points (days 1–4), more transcripts were down-regulated by Leishmania infection at large fold changes (> 32 fold). Among the down-regulated genes, the transcription factor Forkhead/HNF-3 and hormone degradation enzymes were differentially regulated on day 2 and appear to be the upstream regulators of nutrient transport, digestive enzymes, and peritrophic matrix proteins. Conversely, at later time points (days 6 onwards), most of the differentially expressed transcripts were up-regulated by Leishmania infection with small fold changes (< 32 fold). The molecular functions of these genes have been associated with the metabolism of lipids and detoxification of xenobiotics. Conclusion Overall, our data suggest that the presence of Leishmania produces a limited change in the midgut transcript expression profile in sand flies. Further, Leishmania modulates sand fly gene expression early on in the developmental cycle in order to overcome the barriers imposed by the midgut, yet it behaves like a commensal at later time points where a massive number of parasites in the anterior midgut results only in modest changes in midgut gene expression.


Background
Leishmania is a digenetic parasite developing in a mammalian host as well as in an insect vector. These parasites are mostly transmitted by phlebotomine sand flies (Diptera: Psychodidae) of the genera Phlebotomus and Lutzomyia in the Old and New Worlds, respectively [1].
Leishmania develops in the lumen of the sand fly midgut [2][3][4]. Once a sand fly takes up an infected blood meal, round-shaped Leishmania amastigotes (mammalian stage) are carried along within macrophages into the midgut lumen. Between 18 h and 24 h post blood meal, these parasites are released from the macrophages and start to differentiate into procyclic promastigotes within blood enveloped by the peritrophic matrix [5]. During this process, the parasites elongate their cell bodies and expose their flagella, becoming fully differentiated into procyclics by day 2. Between days 2 and 4, Leishmania multiplies and undergoes another differentiation step, acquiring an elongated form (banana-like shape) termed nectomonads [2][3][4]. Upon the breakdown of the peritrophic matrix, the nectomonads escape to the ectoperitrophic space and eventually dock on the midgut microvilli [6,7]. As the remains of the digested blood are evacuated, the parasites detach from the epithelium and further differentiate into the leptomonad stage, which exhibit a smaller cell body and a longer flagellum than nectomonads [2][3][4]. From day 6 onwards, the leptomonads undergo a differentiation process, termed metacyclogenesis, giving rise to the infective forms, the metacyclic promastigotes [8]. During metacyclogenesis, the parasites replace their glycocalyx, exhibiting different sugar side chains on their major surface glycans, reduce the size of their cell bodies, and elongate their flagella [2][3][4]. All these transformations give rise to free swimming highly motile parasites [2][3][4].
Even when developing in their natural sand fly vectors, Leishmania faces barriers imposed by the midgut; overtaking such barriers is critical for the development of mature Leishmania infections. During the transitional stages between amastigotes and procyclic promastigotes, the parasites are susceptible to the harmful action of digestive enzymes [9]. The sand fly immune system may also counteract infection with the parasites, by activation of the Imd pathway [10,11]. Escaping from the peritrophic matrix (PM) is also a crucial step for Leishmania survival [12,13]. Another critical barrier to development is attachment of Leishmania parasites to the midgut epithelium [14]. For this step, specific carbohydrate side chains are required for binding to a midgut epithelium receptor [7,15,16]. From there on, undefined parameters trigger the metacyclogenesis process in parasites leading to the development of a mature infection.
The midgut transcriptomes of three sand fly species have been described. These focused mostly on differences in gene expression triggered by blood intake and parasite infection as compared to sugar fed midguts [17][18][19]. Nonetheless, such studies took place before the advent of deep sequencing and were limited to the investigation of about 1000 transcripts due to the low dynamic range of cDNA libraries. Despite such a limited pool of genes, these studies unveiled multiple genes differentially regulated by blood and/or Leishmania infection. For the later, genes encoding digestive enzymes and components of the peritrophic matrix, the main midgut barriers to Leishmania development, were differentially regulated [17][18][19]. A broader transcriptome of Phlebotomus papatasi cDNA libraries, encompassing multiple stages and time points post-blood feeding, identified 17,120 transcripts and annotated transcripts encoding proteins participating in digestion, PM processes, and immunity [20].
In order to investigate the effects of Leishmania infection on sand fly midgut gene expression, we carried out an RNA-Seq analysis of Leishmania infantum-infected Lutzomyia longipalpis midguts at 7 time points, each corresponding to when the insect midguts are enriched with a particular Leishmania developmental stage. These encompassed early time points when blood digestion is taking place as well as late time points when the parasites are undergoing metacyclogenesis. This approach expands our breadth of knowledge by assessing the effect of Leishmania infection on over 13,000 sand fly midgut transcripts, focusing on genes encoding secreted proteins and also on genes participating in other biological processes of midgut epithelial cells.

Sand fly infection and Leishmania differentiation
Le. infantum growth in the Lu. longipalpis sand fly midgut followed a typical and expected pattern [21]. Briefly, a median of 3000 parasites detected early at 4d Pi increased to 6000 parasites by day 6 and reached about 126,000 parasites at 14d. The proportion of metacyclic stage parasites increased from 0% on 6d to 92% on 14d.
Expanding the Lu. longipalpis midgut repertoire of putative proteins We obtained a total of 53,683,499 high quality Lu. longipalpis midgut-specific reads from the de novo assembly of seven libraries prepared from RNA extracted from uninfected midguts at 1d, 2d, 4d, 6d, 8d, 12d and 14d. High quality reads were assembled in 57,016 contigs that were further filtered to 13,841 putative contigs based on the presence of an open reading frame (ORF) and on similarities to proteins deposited at Refseq invertebrate, NCBI Genbank or SwissProt databases (Additional file 1: Table S1). We also searched for putative secreted proteins where a signal peptide was predicted. The contigs or transcripts from these libraries had a mean size of 1498 bp, with the shortest comprising 150 bp and the longest 27,627 bp. Overall, 72% were categorized to a functional class after BLAST analysis (e < 10E-6) against nine distinct databases (Additional file 2: Figure S1 and Additional file 3: Table S2). Unknown contigs accounted for 28% of contigs, but only for 6.56% of transcriptome abundance.
The annotations of all the 13,841 sand fly midgut transcripts from this work are described in Table S1 (Additional file 1: Table S1), where the best match of each sequence to NCBI, KOG, and Swiss protein databases are shown in different columns. Table S1 also describes if the protein is potentially secreted by giving the term "SIG" in the SigP Result column as a result of SignalP analysis of all transcripts. Using this information, we classified the transcripts into functional categories. The most represented functional categories were secreted proteins with 25.9% of transcripts per million (TPM), protein synthesis (15.2% of TPM), metabolism (14.3% of TPM) and protein modification (11.8% of TPM; Fig. 1 and Additional file 3: Table S2). Importantly, only 0.59% of all transcripts did not have a match in any of the databases tested, indicating that more than 99% of the transcripts from this transcriptome are insect specific transcripts (Additional file 1: Table S1).
The search for transcripts encoding protein families potentially participating in biological processes important for midgut physiology and Leishmania development resulted in 740 sequences encoding proteins associated with immune responses, digestion, and chitin metabolism (Additional file 4: Table S3).
Among the immune-related genes (194 transcripts), the major components of the Toll-like receptor and Imd pathways participating in recognition (GNBPs and PGRPs), signal transduction (Spatzle and TAK1), regulation (cactus and caspar), transcription factors (dorsal and relish), as well as effector molecules (antimicrobial peptides) were identified (Additional file 4: Table S3). Members of the Reactive Oxygen Species (ROX)-producing MAP kinase pathway, such as DUOX, were also identified (Additional file 4: Table S3). Multiple transcripts encoding proteins associated with ROX metabolism (oxidative stress), cell death (JNK pathway and apoptosis) epithelium regeneration (JAK-STAT pathways) were also pinpointed (Additional file 4: Table S3).
The 198 transcripts related to chitin metabolism are possibly involved with the synthesis, scaffolding, modification, and degradation of the PM (Additional file 4: Table S3). Regarding chitin synthesis, two transcripts encoding chitin synthase were identified (Additional file 4: Table S3). Transcripts encoding PM-scaffolding proteins encompassed 28 peritrophins and 7 of 25 chitin-binding protein of the CPAP subgroup (Cuticular proteins analogous to peritrophins), which displayed relevant expression levels in the midgut (Additional file 4: Table S3). Regarding PM modification, 4 transcripts encoding chitin deacetylases were pinpointed (Additional file 4: Table  S3). For PM degradation and chitin digestion, 15 and 6 transcripts encoding chitinases and N-acetylglucosaminidases, respectively, were identified (Additional file 4: Table S3).
Compared to the annotation of the P. papatasi cDNA libraries of whole bodies and multiple stages and time points post-blood feeding [20], multiple additional homolog transcripts were identified in the Lu. longipalpis midgut RNA-Seq libraries, including transcripts encoding PGRPs, C-type lectins, and lysozymes amongst the immune genes. Also, the RNA-Seq libraries displayed more matches of sequences encoding all classes of digestive enzymes (except aminopeptidases), multiple peritrophins, and twice as many chitinase sequences (Table 1).

Sand fly midgut gene expression
The obtained midgut transcriptome dataset was used to determine the sand fly midgut differential expression caused by Leishmania infection. All transcripts used for this analysis can be found in Table S1 (Additional file 1: Table S1.).
We performed a Principal Component Analysis (PCA) to summarize the overall expression profiles of the infected and uninfected midgut transcripts from the seven time points that represent infected midguts enriched with a different Leishmania stage ( Fig. 2a and Additional file 5: Table S4) as well as amongst replicates (Additional file 6: Figure S2A-B and Additional file 5: Table S4). The PC1 axis shows a clear separation of transcripts between the midguts in which blood digestion was ongoing ( Fig. 2a left side, 1d PBM/Pi and 2d PBM/Pi) and less separation from the time points at which the blood was mostly digested (Fig. 2a right side, 4d PBM/Pi) and from the remaining time points where the midguts were clear of blood (Fig. 2a right side, 6d to 14 PBM/Pi). The PC1 accounted for 77.2% of the variance (Additional file 5: Table S4). On the other hand, the PC3, rather than PC2 (Additional file 6: Figure S2A-B), depicted the separation between infected from the uninfected samples (Fig. 2a), accounting for 4.1% of the variance (Additional file 5: Table S4).
The expression profiles of midgut transcripts were validated by assessing the expression levels of selected midgut genes (n = 28-35; Additional file 7: Table S5) using the nCounter technology (NanoString). The mean log 2 fold change (LFC) of infected over uninfected samples was compared at each time point with LFC data obtained by RNA-Seq for the same genes. Representative genes participating in chitin metabolism/peritrophic matrix scaffolding (peritrophins and chitinases), immunity (defensin, catalase, and spatzle), digestion (amylase and chymotrypsin) among others are depicted in Fig. 2bh. The regression analyses between the expression levels obtained with nCounter and RNA-Seq were statistically significant (p < 0.0001) for all seven time points ( Fig. 2bh), and the regression coefficients were greater than 0.5 for all time points, except 6d (R 2 = 0.40) and 12d (R 2 = 0.47) as shown in Fig. 2b-h.

Minimal modulation of sand fly midgut gene expression by Leishmania infection
Differences in midgut gene expression between Leishmania-infected over uninfected midguts were assessed. Overall, such differences accounted for only 113 differentially expressed transcripts (1 < LFC > 1; q-value < 0.05). The number of DE transcripts gradually increased from 2 transcripts on 1d to 53 transcripts on 4d (Fig. 3a). Thereafter, the number of DE transcripts decreased to 20 transcripts on 6d, and 15 transcripts on 8d (Fig. 3a). 4 days later, there was a strong increase in the number of DE transcripts (12d = 32 transcripts), which was reduced to 13 transcripts 2 days later at 14d (Fig. 3a).
The presence of Leishmania in the midgut led to more genes being down-regulated at d2 and up-regulated at later time points, except on 12d ( Fig. 3a and Additional file 8: Table S6). On 1d, 2d, and 4d, the early time points, 1, 11, and 30 genes were down-regulated ( Fig. 3a and Table 3 and Additional file 8: Table S6) whereas 1, 3, and 23 genes were up-regulated ( Fig. 3a and Table 3 and Additional file 8: Table S6). On 6d and 8d, on the other hand, 20 and 15 genes were up-regulated, yet none were down-regulated (Fig. 3a). Infected midguts on day 12 displayed 13 up-regulated compared to 18 down-regulated genes (Fig. 3a). The 14d time point exhibited more up-regulated (12 genes) than down-regulated (1 gene) genes in infected over uninfected midguts (Fig.  3a).
Venn diagrams show the number of DE genes at unique time points, as compared to the number of DE genes shared by multiple time points (Fig. 3b-c and Additional file 9: Table S7). In the comparisons between early time points (1d through 6d; Fig. 3b), 1 out of the 2 DE genes on 1d was only modulated at that time point (Fig. 3b). Similarly, 13 out of the 14 genes, and 43 out of 54 genes, were uniquely DE on 2d and 4d, respectively (Fig. 3b). Only the 6d DE genes exhibited as many unique as shared with 4d DE genes (10 genes; Fig. 3b). The comparisons of DE  genes between later time points (6d through 14d) showed a greater number of shared DE genes between time points (Fig. 3c). For instance, only 5 out of 15, and 5 out of 13, DE genes were unique to 8d and 14d, respectively (Fig. 3c). The 12d midguts, on the other hand, exhibited 26 uniquely expressed genes out 32, the most amongst the late time points (Fig.  3c).

Patterns of differentially expressed genes across time points
Most of the midgut genes DE by Leishmania infection were up-regulated by up to 32-fold (LFC < 5; Table 3; Additional file 10: Figure S3; Additional file 11: Table  S8). These DE genes encoded multiple digestive enzymes and allergen-related peptides at early time points and detoxification-related proteins at later time points (Table  3; Additional file 10: Figure S3). On the other hand, multiple genes were down-regulated in Leishmania-infected midguts by more than 32-fold (LFC > − 5; Table 4; Additional file 10: Figure S3; Additional file 11: Table S8).
Regarding the midgut genes downregulated by Leishmania infection (Table 4 and Additional file 10: Figure  S3), none were DE on 6d and 8d (Table 4). Such genes encode a variety of proteins of unknown function as well as proteins involved in the lipid metabolism (Table 4).
Functional profiles of the differentially expressed genes at different time points Although the midgut genes up-and down-regulated by Leishmania infection exhibited different expression patterns across time points (Additional file 10: Figure S3), such DE genes belonged to the same functional groups for the most part ( Fig. 4 and Tables 3 and 4 and Additional file 11: Table S8). Regarding the up-regulated genes, 28, 38, and 18% belonged to the detoxification (detox), metabolism (met), and secreted (s) protein molecular functions, respectively ( Fig. 4a and Table 3). In fact, the enrichment of such molecular functions amongst the up-regulated genes was consistent through time ( Fig. 4a and Table 3): between 2d through 14d for the metabolism function; and between 8d and 14d for the detoxification function. For the secreted protein category, the enrichment of up-regulated genes was more restricted to 4d and 6d ( Fig. 4a and Table 3). At earlier time points (1d and 2d), the few up-regulated genes perform different functions ranging from transporter channels (tr, 1d) to proteosome machinery (prot, 2d; Fig. 4a and Table 3). Regarding midgut genes downregulated by the Leishmania infection, 34% of these genes belonged to the metabolism (22%) and secreted protein (12%) functional groups ( Fig. 4b and Table 4). Both categories were consistently enriched on 4d, 12d, and 14d ( Fig. 4b and Table 4). At earlier time points (1d and 2d), transporter channels (tr, 1d and 2d) and signaling transduction (st, 2d) were the most enriched molecular functions amongst the down-regulated genes ( Fig. 4b and Table 4). All the molecular functions identified over all time points were matched by analogous GO terms (Additional file 12: Table S9 and Additional file 13: Table  S10).
In order to investigate in-depth the functional profiles of the DE genes, we broke down the most predominant functional classes into subclasses. For the midgut DE genes belonging to the detoxification molecular function (detox), the cytochrome P450 gene family encompassed 76% of the up-regulated genes ( Fig. 4c and Table 3 and Additional file 11: Table  S8). Such genes were consistently upregulated between 6d and 14d ( Fig. 4c and Table 3). In contrast, the down-regulated genes belonging to the detoxification molecular function were enriched in metallothioneins (4d and 12d, thio; Fig. 4d and Table 4). As far as the DE midgut genes belonging to the metabolism function, 55% of the up-regulated genes were related to the metabolism of lipids (lipd; Fig. 4e and Table 3) which was consistently the most predominant between 6d and 14d (Fig. 4e). Among the downregulated genes performing metabolic functions ( Fig.   4f and Table 4), 31% participated in the metabolism of lipids (lipd) at early time points (2d and 4d) or nucleotides (nuc) at later time point (12d and 14d, Fig.  4f and Table 4). Regarding the DE midgut genes encompassing the secreted proteins ( Fig. 4g-h), 50% of those up-regulated belonged to the 'other category' (s, multiple protein functions) that was enriched in transcripts of insect allergen proteins ( Fig. 4g; Table 3 and Additional file 11: Table S8), also known as microvilli proteins. Although the insect allergens, along with the mucins, and to a lesser extent metalloproteases (metal), were more predominant on 4d and 6d ( Fig. 4g and Table 2), up-regulated transcripts encoding proteins of unknown function were enriched at 14d, a later time point (Fig. 4g and Table 3). Among the down-regulated transcripts encoding secreted proteins, 44% belonged to the unknown function (31%, uk) and "other" (17%, s) categories ( Fig. 4h and Table 4). The "other" category (s) was consistently downregulated on 4d and 12d ( Fig. 4h and Table 4) and was enriched in transcripts encoding juvenile hormone (JH) binding proteins as well as attacin ( Table 4). Transcripts of secreted proteins related to the digestion of lipids (met-li) were down-regulated on 2d (Fig. 4h and Table 4).

Discussion
In this work, we carried out a broad RNA-Seq investigation to assess the effects of Leishmania infection in sand fly midgut gene expression. As a sand fly genome is not available for use as a reference for read mapping, all the reads obtained were assembled de novo into 13,841 putative transcripts. Additionally, we manually curated 740 transcripts potentially participating in biological processes important for midgut homeostasis, such as immune responses, digestion, and chitin metabolism. This significant number of likely unique transcripts brings to light multiple homologs to complement the gene set annotated in the cDNA transcriptome study of whole P. papatasi sand flies [20]. Altogether, these transcripts will be of great use for gene prediction and annotation of sand fly genome projects. We also used the 13,841 transcripts as a reference for gene expression quantification and comparisons between infected and uninfected samples. Out of seven time points, only about 1% of the genes were differentially expressed (113 genes) by Leishmania infection, highlighting the extent of the adaptation of Le. infantum to its natural vector, the sand fly Lu. longipalpis.
Multiple midgut genes displaying differential expression upon Leishmania infection in cDNA libraries of Le. infantum-infected Lu. longipalpis midguts [18] were also differentially expressed in our RNA-Seq libraries. For instance, all four insect allergen proteins (microvilli proteins), multiple digestive enzymes (proteases and peptidases), an astacin-metalloprotease, as well as a peritrophic matrix protein were differentially regulated by Leishmania infection in both studies [18].
The limited influence of Leishmania in midgut gene expression as observed in this study was further investigated by PC analysis. As indicated by PC1, most of the  variance (77%) in the transcriptional levels across midgut samples was caused by the presence (or lack of) blood in the midgut, sorting out the early (d1 and d2; blood engorged) from late (d4 onwards; blood passed) time points. Even though PC2 (6.4%) and PC3 (4.1%) exhibited similar levels of variance, PC3 accounted for most of the variance sorting infected from uninfected midguts, and likely represents the differential expression of the 113 genes modulated by Leishmania infection. These findings also suggest that other factors not controlled for by the experimental design accounted for the variance observed in PC2. Along these lines, it is noteworthy that Leishmania infection in the sand fly midguts also modify the microbiota composition [22], which may also have affected gene expression in the midgut samples. It is worth noting that multiple genes DE upon Leishmania infection were unique to a particular time point, being more pronounced in early time points. This phenomenon may be explained by the enrichment of different Leishmania stages at specific time points. For instance, time points 1d, 2d, 4d, and 6d were enriched in amastigotes and transitional stages, and procyclic, nectomonad, and leptomonad promastigotes, respectively. From 6d onwards, Leishmania parasites undergo metacyclogenesis. Hence, there is a gradual increase in the proportions of metacyclic compared to leptomonad promastigotes through time, which can explain the overlap of DE genes between midguts on 8d and the other late time points. Surprisingly, we observed a burst of downregulated DE genes on 12d that was not observed on 14d. At both time points the midgut infection is very similar as far as parasite stage and density, a phenomenon that needs to be further investigated.
In order to complete its life cycle in the sand fly midgut, Leishmania needs not only to develop and differentiate into the infective metacyclic stage, but also to escape the barriers imposed by the sand fly midgut early in the infection (day [1][2][3][4][5]. During this period, Leishmania needs to shield itself against the harmful actions of the proteolytic enzymes [9], avoid the immune system [10,11], escape from the peritrophic matrix [12,13], and attach to the midgut epithelium [14]. At these early time points, most of the sand fly DE genes were downregulated by large fold changes. On day 4, multiple sand fly genes encoding digestive enzyme as well as a peritrophic matrix protein were downregulated, pointing to parasite manipulation of the barriers imposed by the sand fly midgut in order to survive. Such downregulation might also lead to increased availability of nutrients for the parasites. Along the same lines, it is important to highlight that the presence of Leishmania in the sand fly midgut leads to the down regulation of genes potentially involved with the control of gene expression. For instance, the transcription factor Forkhead/HNF-3, which is involved with midgut regeneration [23], and nutrient transport and absorption [24], were down-regulated on day 2. Accordingly, we have also observed down-regulation of sand fly amino acid and trehalose transporters on 4d after Leishmania infection. Transcripts for metallothionein-2-like protein were also down-regulated at the same time point. The expression levels of these proteins are used as a proxy of heavy metals absorption [25]. Hence, their downregulation in Leishmania-infected midguts suggests that these parasites reduce nutrient uptake by the sand fly midgut epithelium. Along the same lines, genes encoding proteins associated with metabolism of hormones, such as the juvenile hormone and ecdysone, were down-regulated on days 4 and 6. Such hormone levels change during blood digestion [26], and relevantly control the expression of midgut serine proteases [27][28][29], which are also downregulated upon Leishmania infection on days 4 and 6. Together, these data suggest that the sand fly transcription factor Forkhead/HNF-3 as well as hormone metabolic enzymes might be key targets to control Leishmania infection early on. As the remains of the digested blood is flushed out and the parasites detach from the epithelium [14], the parasites undergo metacyclogenesis from day 6 onwards, migrating to the anterior midgut and differentiating into infective forms [8]. At this late period in the infection, midgut barriers to Leishmania development are unknown or negligible. The parasites seem to multiply freely, secreting a massive amount of carbohydrates (fPPG) that blocks blood intake and allows the parasites to be regurgitated into the skin [30,31]. Most of the sand fly DE genes late in infection (day 8 onwards) were upregulated by small fold change differences in response to Leishmania. Most of these genes encode proteins that participate in detoxification of xenobiotics (cytochrome P450) and metabolism of lipids. At these time points, it seems plausible that the massive amount of parasites, reaching 120,000 cells on average on day 14 [21], might be indirectly modulating sand fly gene expression by the release of cell membranes and metabolites from dead parasites and Leishmania-derived exosomes [32] throughout metacyclogenesis. Interestingly, the presence of Leishmania is undetected by the midgut immune system of the sand fly during this period. This is also noted at early time points with the exception of day 4 where the down-regulation of a gene encoding attacin, an antimicrobial peptide [33], was observed. The lack of Leishmania detection by the immune system may constitute another adaptation of Le. infantum to survive in Lu. longipalpis midguts.

Conclusion
Overall, the presence of Le. infantum in the midgut of its natural vector has direct and indirect effects on sand fly midgut gene expression. On one hand, these parasites appear to manipulate gene expression early on to weaken developmental barriers imposed by the midgut. On the other hand, Leishmania behaves like a commensal later in the infection, and changes in the sand fly gene expression caused by the parasites seem to be an indirect consequence of the massive amount of the parasites inside the anterior portion of the midgut. Altogether, our findings expose a fine-tuned relationship that has evolved over time to ensure optimal survival and transmission of the Leishmania parasite by vector sand flies.

Methods
Leishmania parasites, parasite load assessment, sand fly blood feeding and infection, and midgut dissection and storage Lutzomyia longipalpis sand flies were obtained from a sand fly colony maintained at the insectary facility, Laboratory of Malaria and Vector Research, National Institute of Allergy and Infectious diseases (LMVR/NIAID). Sand fly infection and Leishmania counts were performed as described in our companion manuscript [21].
As controls, Lu. longipalpis also fed on uninfected heparinized dog blood at the same time. Blood of beagle dogs was provided by the Division of Veterinary Resources (DVR/NIH) under animal use protocol ORS 17 and withdrawn on the same day of the experiments. After feeding, fully fed females were sorted and given 30% sucrose solution ad libitum. In order to assess how gene expression in sand fly midguts is affected by Leishmania growth and differentiation, Le. infantum infected-Lutzomyia longipalpis midguts were dissected at days one, two, four, six, eight, twelve, and fourteen after blood feeding on RNAse Free PBS (1X) for RNA-Seq library construction in triplicate and compared to midguts fed on uninfected blood at the same time points.

RNA extraction and quality control
For total RNA extraction, the PureLink RNA Mini Kit (Life Technologies, Carlsbad) was used as described in [21], following the manufacturer's recommendations. RNA amounts and purity were assessed using a Nanodrop spectrophotometer (Nano Drop Technologies Inc., Wilmingtom; ND-1000), and quality control was further evaluated using a Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA; 2100 Bioanalyzer), using the Agilent RNA 6000 Nano kit (Agilent Technologies) and following the manufacturer's recommendations. Out of the 48 samples, only 1 displayed a low RIN value (RNA integrity number = 6.7; 14d Pi, replicate 3).

RNA-Seq library preparation and deep sequencing
The RNA-Seq libraries were prepared using the NEB-Next® Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs, Ipswick MA), following manufacture's recommendation, for Single Ended sequencing by HiSeq 2500 (Illumina, San Diego, CA) of 125 nucleotides reads (SE -125), and sequenced at the NC State University Genomic Science Laboratory. All the libraries gave rise to high quality data and robust expression levels, except one replicate of the 2d PBM and another of the 12d PBM time points, which were excluded from further analyses.
Bioinformatic pipeline and de novo assembly RNA-seq data trimming and mapping were described elsewhere [21]. De novo assembly from high quality reads were a result of both Abyss (kmers from 21 to 91 in 10-fold increments) and Trinity (V2.1.1) assemblers. The output reads were assembled using an iterative blast and CAP3 pipeline as previously described [34]. Sequences sharing more than 95% identity at nucleotide level were assumed to be alleles of the same genes and clustered together as unique transcripts. Protein coding sequences were defined based on prediction of open reading frame signatures, identification of signal peptide-coding sequence, and by similarities searches in the Refseq invertebrate database from the National Center for Biotechnology Information (NCBI), sequences from Dipterans deposited at NCBI's Genbank and from Swis-sProt. Protein-coding sequences were annotated upon similarity matches to various databases, including Swissprot, Gene Ontology, KOG, Pfam, Drosophila mRNA transcripts, Virus, and SMART, Refseq-invertebrates and the Diptera subset of the GenBank sequences obtained by querying diptera [organism] and retrieving all protein sequences. Raw reads were deposited on the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI). This Transcriptome Shotgun Assembly project has been deposited at DDBJ/ENA/ GenBank under the accession GITU00000000. All sequences used in this work and their corresponding accession numbers are in Additional file 15: Table S12.
Reads were mapped to the generated dataset using the RNA-Seq by Expectation Maximization (RSEM) version 1.3.0, Bowtie version 2-2.2.5, and samtools version 1.2 [35]. Differential expression was analyzed using the Bioconductor package DeSeq2 vs 3.8 [36], using default parameters and the shrinks log2 fold-change (FC) estimation [36,37]. Genes exhibiting sum of counts inferior to 10 across all time points were removed. Differential expression was considered statistically significant under adjusted p-value (P-adj) lower than 5% (p < 0.05) and log 2 fold change higher or lower than 0.5. The whole dataset can be found in Additional file 14: Table  S11.

Filtering and annotation of specific gene families
Transcripts involved in biological processes related to immunity, digestion, and chitin metabolism were filtered upon tBLASTn searches against Drosophila melanogaster [38] and sand fly [18,19] immune-related genes, sand flies digestive enzymes [18,19], and Tribolium castaneum [39][40][41][42] and sand fly [12,13] chitin metabolism homologs. Other members of such gene families were further filtered by manual screening of similar KOG and Swiss databases matches. Only transcripts exhibiting evalues lower than 10 − 5 against most of the databases were annotated. For the annotation of mucin sequences, we relied on motif identification by the automated pipeline described above.