Sand fly infection and Leishmania differentiation
In order to assess how gene expression in sand fly midguts is affected by Leishmania growth and differentiation, Le. infantum infected-Lutzomyia longipalpis midguts, 1 day (1d) through 14 days (14d) post-infection (Pi) were dissected for RNA-Seq library construction in triplicate and compared to midguts fed on uninfected blood at the same time points,1d through 14d, post blood meal (PBM). 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. For infected midguts, Le. infantum growth in the Lu. longipalpis sand fly midgut followed a typical and expected pattern [20]. Briefly, a median of 3,000 parasites detected early at 4d Pi increased to 6,000 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
A Lu. longipalpis midgut-specific de novo assembly was made from libraries prepared from RNA extracted from uninfected midguts at each of the seven study timepoints (total of 53,683,499 high quality reads). High quality reads were assembled in 57,016 contigs that were further down-selected to 13,841 putative contigs based on the presence of an ORF and similarities to proteins deposited at Refseq invertebrate, NCBI Genbank or SwissProt. Putative proteins where a signal peptide was predicted were also considered. Selected contigs varied in size with the shortest comprising 150 bp, the longest at 27,627 bp, and had a mean size of 1,498bp. Overall, 72% could be categorized to a functional class after BLAST analysis (e<10E-6) to nine distinct databases (Additional file 1: Fig. S1 and Additional file 2: Table S1). Figure 1 shows an overview of the transcriptome repertoire displaying the overall percentage of contigs (% of contigs) or abundance as transcripts per million (%TPM) for all time points and conditions combined, highlighting the distribution of the mapped reads to functional classifications Unknown contigs accounted for 28% of contigs, but only for 6.56 % of transcriptome abundance. The most represented functional categories were secreted proteins with 25.9 % of TPM, protein synthesis (15.2 % of TMP), metabolism (14.3 % of TMP) and protein modification (11.8 % of TMP) (Fig. 1 and Additional file 3: Table S2). This dataset was used to map the individual samples and determine the sand fly midgut differential expression caused by Leishmania infection (Additional file 3: Table S2.)
Sand fly midgut gene expression
The overall expression profiles of the infected and uninfected midguts obtained at seven time points each representing infected midguts enriched with a different Leishmania stage is summarized by PCA analyses of the average expression for each time point (Fig 2A and Additional file 4: Table S3) as well as amongst replicates (Additional file 5: Fig. S2A-B and Additional file 4: Table S3). The PC1 axis showed a clear separation between the midguts in which blood digestion is 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 4: Table S3). On the other hand, the PC3, rather than PC2 (Additional file 5: Fig. S2A-B), depicts the separation between infected from the uninfected samples (Fig. 2A), accounting for 4.1% of the variance (Additional file 4: Table S3).
The expression profiles of midguts were validated by assessing the expression levels of selected midgut genes (n = 28-35; Additional file 6: Table S4) using the nCounter technology (NanoString). The mean log2 fold change (LFC) of infected over uninfected samples was compared at each time point with LFC data obtained with the RNA-Seq technique for the same genes. Representative genes participate in chitin metabolism/peritrophic matrix scaffolding (peritrophins and chitinases), immunity (defensin, catalase, and spatzle), digestion (amylase and chymotrypsin) among others are depicted in Fig. 2B-H. 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. 2B-H), and the regression coefficients were greater than 0.5 for all time points, except 6d (R2 = 0.40) and 12d (R2 = 0.47) as shown in Fig. 2B-H.
Modulation of sand fly midgut gene expression by Leishmania infection
Differences in gene expression between Leishmania-infected over uninfected midguts at the seven time points were assessed. Overall, such differences accounted for only 113 differentially expressed transcripts (1 < LFC > 1; q-value < 0.05; Additional file 7: Table S5). The number of DE genes gradually increased from 2 genes on 1d to 53 genes on 4d (Fig. 3A). Thereafter, the number of DE genes decreased to 20 genes on 6d, and 15 genes on 8d (Fig. 3A). Four days later, there was a strong increase in the number of DE genes (12d = 32 genes), which was reduced to 13 genes two days later at 14d (Fig. 3A).
Table 1 Selected midgut transcripts differentially regulated upon Leishmania infection.
Transcript name
|
Best match
|
E-value
|
Time-Point(s)
|
Up/Down Regulated
|
lulogut44569
|
Forkhead/HNF-3-related transcription factor
|
0
|
2d
|
Down
|
lulogut32574
|
17-beta-hydroxysteroid dehydrogenase 13-like isoform X2
|
8E-66
|
2d
|
Down
|
lulogut40195
|
juvenile hormone esterase
|
6.00E-29
|
2d
|
Down
|
lulogutSigP-24104
|
JAV08889.1 juvenile hormone binding protein
|
0
|
4d
|
Down
|
lulogutSigP-40401
|
Chitin binding Peritrophin-A
|
4.00E-12
|
4d
|
Down
|
lulogutSigP-8812
|
attacin precursor
|
5.00E-64
|
4d
|
Down
|
lulogut16004
|
Amino acid transporters
|
0
|
4d
|
Down
|
lulogutSigP-40100
|
Facilitated trehalose transporter Tret1
|
5.00E-93
|
4d
|
Down
|
lulogutSigP-25516
|
chymotrypsin-2
|
8.00E-80
|
4d
|
Up
|
lulogutSigP-33169
|
Trypsin-like serine protease
|
4.00E-67
|
4d
|
Up
|
lulogutSigP-12857
|
carboxypeptidase A
|
0
|
4d/6d
|
Up
|
lulogutSigP-53922
|
Secreted metalloprotease
|
0
|
6d
|
Up
|
lulogutSigP-646
|
Insect allergen related repeat
|
5.00E-28
|
4d
|
Up
|
lulogutSigP-16736
|
Insect allergen related repeat
|
4.00E-30
|
4d/6d
|
Up
|
lulogutSigP-13949
|
Insect allergen related repeat
|
2.00E-42
|
4d/6d
|
Up
|
lulogutSigP-13652
|
Insect allergen related repeat
|
2.00E-32
|
4d/6d
|
Up
|
lulogutSigP-54492
|
Insect allergen related repeat
|
5.00E-42
|
6d
|
Up
|
lulogutSigP-8474
|
probable cytochrome P450 6a14
|
0
|
8d/12d/14d
|
Up
|
lulogut46050
|
cytochrome P450 4C1
|
0
|
8d
|
Up
|
lulogut33084
|
Cytochrome P450 CYP3/CYP5/CYP6/CYP9 subfamilies
|
0
|
12d
|
Up
|
lulogut34615
|
probable cytochrome P450 6d5
|
0
|
8d/14d
|
Up
|
lulogut41307
|
JAV11511.1 ecdysteroid kinase
|
0
|
12d
|
Down
|
Amongst the midgut genes differentially expressed upon Leishmania infection, some appear to play a role in specific biological processes (Table 1; Additional File 7: Table S5). A gene encoding the transcription factor Forkhead/HNF-3 (lulogut44569) was down-regulated on 2d. Genes encoding proteins potentially involved with metabolism of steroid hormones, such as 17-beta-hydroxysteroid dehydrogenase 13-like (lulogut32574) and juvenile hormone esterase (lulogut40195) were down-regulated on 2d; a putative juvenile hormone binding protein (lulogutSigP-24104) was down-regulated on 4d; and an ecdysteroid kinase (lulogut41307) was down-regulated on 12d. Also, genes encoding a peritrophic matrix protein (lulogutSigP-40401), involved with the peritrophic matrix scaffolding, the antimicrobial peptide attacin (lulogutSigP-8812), and amino acid (lulogut16004) and trehalose (lulogutSigP-40100) transporters, were down-regulated on 4d. Amongst the up-regulated genes, multiple peptidases and proteases were up-regulated on 4d and 6d. Likewise, multiple insect allergen proteins (microvilli proteins) of unknown function were up-regulated on 4d and 6d upon Leishmania infection. From 8d onwards, multiple cytochrome p450 transcripts were up-regulated.
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 7: Table S5 and Additional File 8: Fig. S3). On 1d, 2d, and 4d, the early time points, 1, 11, and 30 genes were down-regulated (Fig. 3A and Table 2 and Additional File 7: Table S5 and Additional File 8: Fig. S3) whereas 1, 3, and 23 genes were up-regulated (Fig. 3A and Table 3 and Additional File 7: Table S5 and Additional File 8: Fig. S3), respectively. On 6d and 8d, on the other hand, 20 and 15 genes were up-regulated, yet none were down-regulated (Fig. 3A and Additional File 8: Fig. S3). Infected midguts on day 12 displayed 13 up-regulated compared to 18 down-regulated genes (Fig. 3A and Additional File 8: Fig. S3). The 14d time point exhibited more up-regulated (12 genes) than down-regulated (1 gene) genes in infected over uninfected midguts (Fig. 3A and Additional File 8: Fig. S3).
Venn diagrams show that most genes were differentially expressed at specific time points (Fig. 3B-C and Additional File 9: Table S6). 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).
The expression patterns of all DE genes across time points were assessed through PCA analysis (Fig. 3D) and Additional file 10: Table S7). The 113 DE genes were mapped onto a two-dimensional space (expression space), whereby DE genes located close together displayed similar expression profiles through time than those that mapped farther away (Fig 3D). In fact, the DE genes located in the first quadrant of the expression space exhibited about 25-fold greater overall expression levels than those that mapped onto the second and third quadrants (Fig. 3E; Mann Whitney U test, p < 0.0001). Likewise, the DE genes located on the fourth quadrant of the expression space exhibited about 177-fold higher overall expression levels than those that mapped onto the second and third quadrants (Fig. 3E; Mann Whitney U test, p < 0.0001). Looking through time, the location of the DE genes in different quadrants further highlighted temporal expression differences in both early (blood-engorged) infected midguts and late time point (blood-passed) infected midguts (Fig. 3F; Additional file 11: Fig. S4). For example, the DE genes mapping onto the first and second quadrants were either down-regulated at early time points (1d and/or 2d) and up-regulated at later time points (d4 onwards; Fig. 3F and Additional file 11: Fig. S4). On the other hand, the DE genes located on the third and fourth quadrants were up-regulated at 1d and 2d and down-regulated from 4d onwards (Fig. 3F and Additional file 11: Fig. S4). Hence, DE genes located on the first quadrant were expressed at high abundance and were late up-regulated; DE genes mapping onto the second quadrant were expressed transcripts at low abundance and were up-regulated at late time points; the third quadrant housed the DE genes expressed at low abundance and up-regulated at early time points; and the DE genes transcribed at high abundance and up-regulated at early time points were localized on the fourth quadrant of the transcriptional space (Fig. 3D-E).
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 (Table 2, LFC < 5). Interestingly, most of the DE midgut genes up-regulated by Leishmania infection encompassed highly expressed genes, yet the up-regulated genes were more predominant at early time points (4d and 6d) and the late expressed genes were more predominant at late time points, 8d to 14d (Table 2).
The up-regulated genes (Fig. 4A-G and Table 2 and Additional file 12: Figure S5 and Additional file 13: Table S8) at each time point were plotted onto the transcriptional space in order to assess whether or not the expression of the genes modulated by Leishmania across time points followed a specific or a random expression pattern by mapping onto specific quadrants or randomly. All the 113 DE genes were distributed throughout the four quadrants in different proportions: 20%, 37%, 22%, and 21% of the DE genes mapped onto the first through fourth quadrants, respectively (Fig. 4 and Additional file 12: Figure S5 and Additional file 13: Table S8). The DE genes at 4d through 8d followed specific expression patterns (Chi-square test, p < 0.01; Fig. 4C-E). At such time points, 74% (4d), 85% (6d), 87% (8d) of genes up-regulated by Leishmania infection mapped onto either the first or fourth quadrant (Additional file 12: Figure S5 and Additional file 13: Table S8), which housed genes transcribed at high abundance. These DE genes encoded multiple digestive enzymes and allergen-related peptides at early time points and detoxification-related proteins at later time points (Table 2).
Many of the genes were down-regulated in Leishmania-infected midguts by more than 32-fold (LFC > -5; Table 3). Regarding the midgut genes down-regulated by Leishmania infection (Fig. 5A-E and Table 3 and Additional file 14: Figure S6 and Additional file 15: Table S9), none were DE on 6d and 8d (Table 3). Contrasting to the midgut up-regulated genes, which exhibited similar expression profiles and were fine-tuned through time (Fig. 4C-G), for the most part the midgut down-regulated genes displayed more diverse expression patterns, highlighted by the random distribution of such genes in the transcriptional space on 1d and 2d (Fig. 5A-B; Table 3), and 12d and 14d (Fig. 5F-G; Table 3). On the other hand, the 4d midguts displayed most of the down-regulated genes on the second quadrant (73%, p < 0.05; Fig. 5C and Table 3), belonging to the group transcribed at low abundance and up-regulated late in infection (Fig. 5C). Such genes encode a variety of proteins of unknown function as well as proteins involved in the lipid metabolism (Table 3).
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 (Figs. 4 and 5), such DE genes belonged to the same functional groups for the most part (Fig. 6 and Tables 2 and 3 and Additional file 16: Table S10). Regarding the up-regulated genes, 28%, 38%, and 18% belonged to the detoxification (detox), metabolism (met), and secreted (s) protein molecular functions, respectively (Fig. 6A and Table 2). In fact, the enrichment of such molecular functions amongst the up-regulated genes was consistent through time (Fig. 6A and Table 2): 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. 6A and Table 2). 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. 6A and Table 2). Regarding midgut genes down-regulated by the Leishmania infection, 34% of these genes belonged to the metabolism (22%) and secreted protein (12%) functional groups (Fig. 6B and Table 3). Both categories were consistently enriched on 4d, 12d, and 14d (Fig. 6B and Table 3). 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. 6B and Table 3). All the molecular functions identified over all time points were matched by analogous GO terms (Additional file 17: Table S11 and Additional file 18: Table S12).
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. 6C and Table 2 and Additional file 16: Table S10). Such genes were consistently up-regulated between 6d and 14d (Fig. 6C and Table 2). In contrast, the down-regulated genes belonging to the detoxification molecular function were enriched in metallothioneins (4d and 12d, thio; Fig. 6D and Table 3). 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. 6E and Table 2) which was consistently the most predominant between 6d and 14d (Fig. 6E). Among the down-regulated genes performing metabolic functions (Fig. 6F and Table 3), 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. 6F and Table 3). Regarding the DE midgut genes encompassing the secreted proteins (Fig. 6G-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. 6G; Table 2 and Additional file 15: Table S9), 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. 6G and Table 2), up-regulated transcripts encoding proteins of unknown function were enriched at 14d, a later time point (Fig. 6G and Table 2). Among the down-regulated transcripts encoding secreted proteins, 44% belonged to the unknown function (31%, uk) and “other” (17%, s) categories (Fig. 6H and Table 3). The “other” category (s) was consistently down-regulated on 4d and 12d (Fig. 6H and Table 3) and was enriched in transcripts encoding juvenile hormone (JH) binding proteins as well as attacin (Table 3). Transcripts of secreted proteins related to the digestion of lipids (met-li) were down-regulated on 2d (Fig. 6H and Table 3).