Overview of Samples Sequencing
The purpose of this study was to compare the transcriptome of SbIII resistant (LiSbR) and wild-type (LiWTS) L. infantum lines. For this, cDNA libraries from these samples were constructed, sequenced and analyzed, allowing the identification of differential gene expression associated with SbIII resistance mechanisms.
Three independent biological replicates of the wild-type parasites and two of the SbIII resistant parasites were sequenced, producing ~500 million of 100 base pair reads. After quality trimming (adaptors removal and Phred quality cutoff greater than or equal to 25), approximately 2% of the reads were lost. After mapping, approximately 388 million reads were aligned to a reference genome (L. infantum JPCM5).
Differential Expression analysis
In the dataset comprising 8591 protein coding transcripts (obtained from ENA database), 933 (933/8591, 10.86%) DE transcripts were identified between wild-type and SbIII-resistant L. infantum lines considering the applied cutoffs of adjusted p-value less than 0.05 and FC greater than 2.0. Out of 933 DE transcripts, 504 (504/933, 54.01%) presented functional annotation and 429 (429/933, 45.99%) were assigned as hypothetical proteins without predicted function (Table 1). A total of 837 (837/933, 89.71%) transcripts were upregulated and 96 (96/933, 10.29%) were downregulated in the SbIII-resistant L. infantum (Table 1).
The functional enrichment analysis of 933 transcripts obtained in this study was performed on Blast2GO software [38,39]. Out of 933 transcripts, 644 (644/933, 69.02%) were associated with some GO ontology related to biological process (BP), molecular function (MF) or cellular component (CC) and 289 (298/933, 30.98%) did not present any GO term associated (Table 1). A total of 207 (207/644, 32.14%) DE transcripts presented GO ontology on biological processes, being 181 (181/207, 87.44%) functionally enriched and 26 (26/207, 12.56%) did not show enrichment (Table 1).
The distribution of 644 differentially expressed transcripts in the three different GO categories, BP, MF and CC, is represented in the Venn diagram (Fig. 1).
Many transcripts are associated with more than one GO ontology. The majority of transcripts correspond to cellular components (282/644, 47.79%), followed by biological processes (147/644, 22.83%) and molecular function (127/644, 19.72%). Fourteen transcripts present all three categories.
Gene Ontology enrichment analysis (Fisher’s exact test) was performed using as “test set” the list of upregulated (Fig. 2) and downregulated (Fig. 3) differentially expressed transcripts (DE) and as “reference set” (background) the L. infantum JPCM5 predicted proteome. A FDR less than 0.05 was set as the threshold to define the functional enrichment significance.
Gene Ontology enrichment analysis of 837 transcripts upregulated in the SbIII-resistant L. infantum showed that the overrepresented terms were related to quorum sensing (GO:0009372, GO:0052097 and GO:0052106), host-parasite interaction (GO:0044764, GO:0044114, GO:0044115 and GO:0044145), protein modifications (GO:1903320 and GO:0032446), post-translational modifications, protein phosphorylation (GO:0006468) and protein ubiquitination (GO:0016567), microtubule-based movement (GO:0007018) and regulation of membrane lipid distribution (GO:0097035) (Figs. 2, 4).
In contrast, GO enrichment analysis of 96 transcripts downregulated in the SbIII-resistant L. infantum showed overrepresentation of all GO terms linked with rRNA processing (GO:0006364), nucleosome assembly (GO:0034622, GO:0022607, GO:0065003, GO:0022618, GO:0070925, GO:0042255, GO:0000028, GO:0006334, GO:0031497, GO:0006333, GO:0065004 and GO:0034063) and maintenance of translational fidelity (GO:1990145) (Figs. 3, 4).
RNA profiling of L. infantum (MHOM/BR/74/PP75)
Genes upregulated in the LiSbR line
Out of 431 (431/834, 51.19%) upregulated enriched genes in the SbIII-resistant L. infantum line, 124 (124/431, 28.77%) presented GO ontology on biological process (111 enriched and 13 without enrichment), 289 (289/431, 67.05%) genes did not present GO ontology on biological process and 18 (18/431, 4.18%) genes had not GO ontology associated (Additional file 1: Table S1, Additional file 2: Table S2 and Additional file 3: Table S3, respectively).
According to GO enrichment analysis, some terms related to biological processes were under or overrepresented in the differentially expressed genes (Fig. 2). The most representative GO terms were protein phosphorylation, microtubule-based movement, protein ubiquitination, cellular process, quorum sensing involved in interaction with hosts and others (Fig. 4). Data from other DE genes up and downregulated in the LiSbR line that presented GO enrichment are described on Tables 2 and 3. These data are representative of complete results depicted on Additional file 1: Table S1.
Thirty seven transcripts belonging to the protein phosphorylation category were upregulated in the LiSbR line (Table 2, Additional file 1: Table S1 and Additional file 2: Table S2). This group includes five transcripts encoding phosphatidylinositol kinase (PIK) (2.52 to 3.97-fold upregulated); RAB GTPases (2-fold upregulated); dual specificity protein phosphatase (DUSP) (8.93-fold upregulated); protein phosphatase and protein phosphatase 2C (17.52 to 2.28-fold upregulated, respectively); cyclins 10, 11 and CYC2-like (2.27 to 5.56-fold upregulated) and elongation factor 2 (EF2) (3.23-fold upregulated).
In the microtubule-based movement category, 20 transcripts were upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1), including dyneins (2.04 to 5.9-fold upregulated); ten transcripts encoding kinesins (2.02 to 5.87-fold upregulated); tryptophan-aspartic acid (WD) protein (2-fold upregulated) and tetratricopeptide repeat domain (TPR) (2-fold upregulated).
GO enrichment analysis also showed that transcripts related to protein ubiquitination were upregulated in the LiSbR line. Twenty transcripts (2.03 to 9.14-fold upregulated in the LiSbR line) were assigned for this category, such as: ubiquitin, ubiquitin-transferase, cullin protein and zinc finger containing proteins. Four transcripts encoding different zinc finger proteins (C3HC4 type - RING finger and FYVE) were 2.15 to 3.77-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). Glycosomal transporter (GAT3) was 3.39-fold upregulated in the LiSbR line.
Other categories were also present among the upregulated transcripts. Serine palmitoyltransferase, included in the biosynthetic process category, was 2.99-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). Transcripts encoding ATP-dependent RNA helicase were 2.38 to 3.34-fold upregulated in the LiSbR line (Table 2, Additional file 1: Table S1 and Additional file 2: Table S2) and were included in the ribonucleoprotein complex assembly category. In the stress granule assembly category, one transcript assigned as pumilio protein was 5.59-fold upregulated in the LiSbR line. Four transcripts related to phospholipid-transporting ATPase/P-type were upregulated in the LiSbR line. In the cellular metabolic process category, a transcript encoding a 100 kDa heat shock protein was 2.86-fold upregulated in the LiSbR line. In the categories: primary metabolic process, cellular macromolecule biosynthetic process and cellular nitrogen compound metabolic process, DNAJ was found to be 2.13-fold upregulated in the LiSbR line. Many transcripts belonging to ATP-binding cassette (ABC) transporters were 2.2 to 4.6-fold upregulated in LiSbR line and were included in the categories regulation of membrane lipid distribution and phospholipid translocation. Among the transcripts implicated with quorum sensing involved in interaction with host and multi-organism cellular processes, four transcripts of RNA recognition motif (RRM) were 2.87 to 14.4-fold upregulated in the LiSbR line.
Genes downregulated in the LiSbR line
Out of 73 (73/96, 76.01%) enriched transcripts downregulated in the SbIII-resistant L. infantum line, 37 (37/73, 50.68%) presented GO enrichment on biological process, six (6/73, 8.22%) genes did not present GO ontology on this category and 30 (30/73, 41.09%) genes had no GO term associated (Additional file 1: Table S1, Additional file 2: Table S2 and Additional file 3: Table S3, respectively). According to GO enrichment analysis, some terms related to biological processes were under or overrepresented in the DE genes (Fig. 3, Table 3 and Additional file 1: Table S1). The most representative GO terms were ribonucleoprotein complex subunit organization, rRNA processing, ribosome biogenesis, translation and nucleosome assembly.
According to GO enrichment analysis, a group of terms related to ribosome were overrepresented in the DE dataset. The transcripts encoding ribosomal proteins such as: ribosomal proteins 40S and 60S, nucleolar and nuclear proteins are downregulated in the LiSbR line.
Fourteen transcripts encoding ribosomal proteins of the small ribosomal 40S subunit, 40S ribosomal S3a, S4, S15, S16, S17, S18, S19, S21, S23 and S33 were 2.01 to 3.12-fold downregulated in the LiSbR line (Table 3). In addition, 11 transcripts encoding ribosomal protein components of 60S subunit of large ribosomal subunit: 60S ribosomal L5, L7a, L11, L13, L18a, L21, L22, L31, L35 and L37 were 2.02 to 3.0-fold downregulated in the LiSbR line (Table 3).
Transcripts encoding the histones H2A, H2B, H3 and H4 were found 2.0 to 2.56-fold downregulated in antimony resistant L. infantum line (Table 3, Additional file 1: Table S1 and Additional file 3: Table S3), and were included in the nucleosome assembly category.
Proteins without GO enrichment for biological process
Some differentially expressed transcripts were not related to any category in the GO enrichment analysis (Additional file 2: Table S2), including, for example: 60S ribosomal L23a (2.07-fold upregulated in the LiSbR line); cytochrome b5 and cytochrome P450 reductase (6.37 and 4.33-fold upregulated in the LiSbR line, respectively); gamma-glutamylcysteine synthetase (2.6-fold upregulated in the LiSbR line); mannosyltransferase (2.54-fold upregulated in the LiSbR line) and two transcripts encoding protein classified as amastin (3.33 and 2.97-fold upregulated in the LiSbR).
Data from comparative transcriptomic analysis of susceptible and SbIII-resistant L. infantum lines showed that 429 differentially expressed transcripts were assigned as hypothetical protein without predicted function. From these, 406 transcripts were upregulated and 23 were downregulated in the LiSbR line (Table 1). Out of 429 DE transcripts, 46 presented GO ontology on biological process (32 functionally enriched and 13 without enrichment), 142 transcripts did not present GO ontology on biological process and 241 genes had no GO term associated (Additional file 4: Table S4). According to GO enrichment analysis, some terms related to biological processes were under or overrepresented in the transcripts differentially expressed (Figs 2 to 4). Similar to DE genes, the main terms enriched were microtubule-based movement, protein phosphorylation, protein ubiquitination, quorum sensing involved in interaction with host, ribonucleoprotein complex and others.