Tsetse flies were sampled in May and June 2015 in two active HAT foci (Campo and Bipindi), located in the Ocean Division of the southern region of Cameroon. The Campo focus (2°20′N, 9°52′E) is located on the Atlantic coast and extends along the Ntem river. The HAT National Control Program that visits Campo once per year diagnosed 61 novel HAT cases between 2001 and 2011. The passive identification of two cases in 2012 (37) indicates that HAT is still present. The Bipindi focus (3°2′N, 10°22′E) has a typical forest bioecological environment, including equatorial forest and farmland along roads and around villages. This focus has been recognized since 1920 (38) and includes several villages. Sleeping sickness is still present, since approximately 83 HAT cases were identified by the National Control Program in this focus between 1998 and 2011 (Ebo’o Eyenga, personal communication). In addition to HAT cases that involve G. palpalis and Tbg, regular global surveys have identified the presence of several other Glossina (including Gpp) and Trypanosoma species (including Tc) in both foci. Surveys have also identified a variety of domestic and wild animals that serve as reservoirs for diverse Trypanosoma species (39–42). As described below, flies were trapped in these areas in order to select non-infected and Tc-infected individuals.

The May 2015 tsetse fly trapping campaign was conducted in three Campo villages (Ipono, Mabiogo, and Campo-Beach), and the June 2015 campaign was conducted in three Bipindi villages (Lambi, Bidjouka, and Ebiminbang). The geographical positions of the sampling sites were determined by GPS. Tsetse flies were captured using pyramidal traps (43) placed in suitable tsetse fly biotopes. Each trap was installed for four consecutive days, and the flies were collected twice per day.

Prior to handling samples, work stations and dissecting instruments were cleaned with RNase away (Ambion) in order to eliminate any RNases that could degrade sample RNA. Furthermore, tsetse flies were dissected alive to prevent RNA degradation by normal post mortem degradation processes. The first step in sample processing consisted in identifying the collected tsetse flies to the species level on the basis of morphological criteria and adapted taxonomic keys (44). Next, the samples were separated into two groups of teneral and non-teneral flies. The non-teneral Gpp flies were dissected in a drop of sterile 0.9% saline solution, according to the midgut dissection protocol developed by Penchenier and Itard (45). The organs were immediately transferred to tubes containing RNAlater (Ambion) for DNA and RNA extraction. These samples were then used for parasite identification by specific PCR amplification, and ultimately for transcriptomic analysis. All tools were carefully cleaned after the dissection of each fly to prevent cross-contamination. During field manipulations, the tubes containing the organs were stored at −20°C for 5 days; subsequently, they were stored in the laboratory at −80°C until use.

To prepare for extraction, samples stored at −80°C were thawed and RNAlater was removed. The midguts were treated with the NucleoSpin TriPrep extraction kit (Macherey-Nagel) according to the manufacturer’s instructions, which allow the separate extraction of DNA and RNA. RNase free water (40 µl) was added to elute the RNA, and 100 µl of DNA elute solution was added to recover the DNA. RNA quality and the absence of any DNA contamination were checked on an Agilent RNA 6000 Bioanalyzer and quantified using the Agilent RNA 6000 Nano kit (Agilent Technologies, France).

To identify which trypanosome species had infected the sampled tsetse flies, the isolated DNA samples stored at −80°C were thawed and used as a template for PCR amplification with specific primers (Table 1). PCR amplification of parasites was performed as described by Herder et al. (46) and consisted of a denaturing step at 94°C (5 min) followed by 44 amplification cycles, each comprising a denaturing step at 94°C (30 s), annealing at 55°C (30 s), and an extension step at 72°C (1 min). A final extension was performed at 72°C for 10 min. The amplified products were separated on a 2% agarose gel containing ethidium bromide and visualized under UV illumination. Positive (2 ng of reference DNA) and negative controls were included in each PCR amplification experiment. PCR amplifications that gave a positive result were repeated once for confirmation.

A total of 1,991 tsetse flies were collected during the entomological survey (775 flies from Campo and 1,216 flies from Bipindi). The Campo fly population was composed of Gpp (95.61%), Glossina caliginea (2.06%), Glossina palicera (1.87%), and Glossina nigrofusca (0.52%). Two tsetse fly species were identified at the Bipindi focus: Gpp (99.34%) and G. palicera (0.66%). The mean apparent density was 4.24 flies per trap per day; however, this parameter was highly variable between the different villages and was higher in Bipindi (8.1) than in Campo (3.52) (Table 2). The frequency of teneral flies was typically low in both Bipindi (0.08%) and Campo (1.16%). These data are roughly in line with data recorded in 2007/2008 (13), although the rate of teneral flies was much lower in the present study. Only 1,245 of the trapped 1,991 tsetse flies were dissected, since 10 flies were teneral and 736 flies were desiccated.

The number of flies carrying single or mixed trypanosome infections is presented in Table 3. Of the 337 Campo flies analyzed, 25 (7.41%) were infected by the Tc “forest” type, 16 (4.74%) by the Tc “savanah” type, and 14 (4.15%) by both parasites. In contrast, Bipindi flies only carried the Tc “forest” type (8.33%). Table S1 in Supplementary Material details the characteristics of the different samples including those used for transcriptomic analyses.

The sequencing of libraries produced a total of 400 million reads (theoretically 40 million reads per sample), which represent a satisfactory sequencing depth for subsequent differential gene expression analysis (Figure 1). Out of the 328 raw clusters generated, 77.6% were successfully filtered, with each sample producing 50–72 million clusters (mean: 61 million clusters) (Figure 1). Sequencing also revealed a total of 31,320 contigs distributed in 3,926 scaffolds with a mean size of 96,817 bp (varying in size from 545 bp to 3.6 Mb).

RNA-seq sequencing produced an average of 124 million reads per sample. From this, 111.45 million reads (83.6%) were mapped onto the genome of Gpp, 103.68 M (73.7%) of which were properly paired (Table 4).

Reads were also mapped onto the trypanosome genome in order to validate the infection status of the samples, as well as to investigate the role played by the trypanosome in this molecular dialog. Since trypanosome cells represented a small fraction of the analyzed material, only a smaller fraction of the reads could be mapped. Specifically, 300,302.5 (0.24%) reads were mapped from an average of 124 million reads per sample. This resulted in 137,305.4 (0.11%) properly paired reads and 157,712.5 (0.12%) singletons (Table 4).

We used DESeq2 (Table 5) to detect genes that were differentially expressed between the five infected and five non-infected samples. When a Benjamini–Hochberg corrected p-value lower than 0.05 was applied, 524 genes were observed to be significantly differentially expressed in infected vs. non-infected flies, among which 285 genes were downregulated and 239 were upregulated. A similar DEG analysis was performed using edgeR (Table 5), which identified only 20 downregulated genes and 53 upregulated genes. Figure 2 presents the volcano plots produced by DESeq2 (Figure 2A) and edgeR (Figure 2B); genes that were significantly (p-value <0.05) differentially expressed and with a fold change of log2 (fold change) >2 (i.e., upregulated genes) or log2 (fold change) <−2 (i.e., downregulated genes) were considered relevant.

To understand the roles of DEGs associated with tsetse fly infection by Tc, identifiers of the Drosophila orthologs of the Glossina DEGs were examined using the DAVID-functional enrichment tool. Separately, we analyzed the 290 downregulated and 213 upregulated genes reported by DESeq2 and identified 207 Drosophila best hits (121 downregulated and 86 upregulated genes). The same analysis by edgeR only resulted in 25 DEG orthologs (6 downregulated and 19 upregulated genes). This reduced number of orthologs could possibly be due to the incomplete assembly and annotation of the Glossina genome and/or the high stringency of the BBH criterion (which discards the case where several tsetse fly proteins have the same closest hit in Drosophila). These DEGs were examined using DAVID, which compares the list of input genes with a variety of functional annotations. This analysis was focused on the three primary categories of the Gene Ontology annotation: biological process (BP), molecular function (MF), and cellular component (CC). The list of different features is provided in Table 6.

The 285 downregulated genes identified by DESeq2 in flies infected with Tc mainly belonged to the BP category, in which the major functional classes were RNA processing (58 genes; 23.48%), ribosome biogenesis (24 genes; 9.71%), and translation (13 genes; 5.26%); the other genes corresponded to several poorly represented classes. The MF category included the functional RNA binding classes (20 genes; 8.09%) and catalytic activity (25 genes; 10.12%). Finally, the CC category included the intracellular lumen (48 genes; 19.43%), non-membrane-bound organelle (34 genes; 13.75%), and ribonucleoprotein complex (25 genes; 10.12%) functional classes (Figure 3).

In addition, 239 DEGs were overexpressed in Tc-infected tsetse flies. These DEGs encoded proteins corresponding to the same three primary ontology categories. The BP category include the neuron morphogenesis functional class (83 genes; 38.96%), amino acid biosynthesis (24 genes; 11.26%), and carboxylic acid biosynthesis (10 genes; 4.7%). The MF category included the iron binding (64 genes, 30%) and catalytic activity (24 genes; 11.26%) functional classes. Finally, the CC category only included the mitochondrion functional class (8 genes; 3.75%).

As already shown for DEGs, edgeR provides a much lower number of functional annotations than DESeq2. Here, using edgeR, fly genes that displayed an increased expression in response to Tc infection were found to belong to the MF category, with only three functional classes: phospholipase activity (2 genes; 10.5%), lipase activity (2 genes; 10.5%), and carboxylesterase activity (2 genes; 10.5%). Finally, expression was decreased for only two genes belonging to the one-carbon metabolic process term (BP category).

To make our analysis more focused and efficient, a functional enrichment was performed to refine the list of tsetse fly DEGs in which expression was influenced by Tc infection. We therefore combined fold enrichment and the p-value at a 5% threshold, which allowed applying a Bonferroni correction to eliminate false positives. The Bonferroni correction threshold was fixed at α = 10⁻², and all functionality with a Bonferroni value below this threshold was considered to be due to trypanosome infection.

Following this correction, 16 functional classes were found to be selectively altered by trypanosome infection. These classes are mainly involved in the transcription process, including (a) RNA related processes (rRNA processing, rRNA metabolic process, ncRNA processing, ncRNA metabolic processes, and RNA processing), involving 63 DEGs; (b) monitoring processes related to the synthesis of the ribonucleoprotein complex (ribonucleoprotein complex biogenesis and ribonucleoprotein complex), involving 29 DEGs; (c) RNA binding, involving 24 DEGs; (d) nucleolus biogenesis (nucleolus and nuclear lumen), involving 23 DEGs; (e) ribosome synthesis, involving 17 DEGs; and (f) eukaryotic translation factor 3 complex, involving 4 DEGs. In contrast to the transcription process, which involved 160 DEGs, the BPs of organic acid synthesis (amine biosynthetic process, cellular amino acid biosynthetic process, carboxylic acid biosynthetic process, and organic acid biosynthetic process) that were found to be activated by trypanosome infection only involved 21 DEGs (Table 7).

The BioCyc metabolic map (Figure 4) and Table 8 illustrate the different pathways that the DEGs are involved in. Among these, the amino acid biosynthesis pathway (which includes the biosynthesis of l-glutamine, l-glutamate, l-serine, l-asparagine, l-aspartate, etc.) is controlled by genes that were shown to be overexpressed following trypanosome infection in the flies. Similarly, DEGs associated with the nucleotide biosynthesis pathway were overexpressed following trypanosome infection, especially uridine monophosphate, an RNA monomer. In contrast, genes involved in the pentose phosphate pathway, namely those implicated in the synthesis of d-ribose 5-phosphate, were downregulated. Phosphorylated pentose is converted by ribose phosphate diphosphokinase into phosphoribosylpyrophosphate, a precursor of nucleotide synthesis. Finally, regarding the carbohydrate biosynthesis pathway, we observed an overexpression of genes encoding malate dehydrogenase, which converts malate into pyruvate.

Interestingly, a large number of up- and downregulated DEGs were related to a given biosynthetic process meaning that upregulated DEGs encoding amino acids (such as tyrosine, serine, glutamine, and several others) were found in the amino acid synthesis pathway. Other overexpressed genes that were identified are involved in the biosynthesis of galactosyltransferase, N-acetylglucosamine, and beta-1,4-manosylglycolipide, which are all molecules that interact with the immune system of the tsetse fly (64–66). In contrast, genes involved in the biosynthesis of cytokines were downregulated in trypanosome-infected flies, as compared to non-infected flies. Similarly, genes involved in folate metabolism (e.g., the biosynthesis of formyltetrahydrofolate dehydrogenase), the main source of energy in flies, were downregulated.

Finally, DEGs involved in the transport of several molecules from the extracellular space toward the cytosol compartment were upregulated. This transport includes molecules with a role in cell nutrition, and nutrients such as lipids, but also ATP, succinate, and l-carnitine (which participates in the degradation of fats).