# Horizontally Transferred Genetic Elements in the Tsetse Fly Genome: An Alignment-Free Clustering Approach Using Batch Learning Self-Organising Map (BLSOM)

Authors: Ryo Nakao; Takashi Abe; Shunsuke Funayama; Chihiro Sugimoto
Journal: BioMed Research International
 (2016)

---

## Abstract

Tsetse flies (Glossina spp.) are the primary vectors of trypanosomes, which can cause human and animal African trypanosomiasis in Sub-Saharan African countries. The objective of this study was to explore the genome ofGlossina morsitans morsitans for evidence of horizontal gene transfer (HGT) from microorganisms. We employed an alignment-free clustering method, that is, batch learning self-organising map (BLSOM), in which sequence fragments are clustered based on the similarity of oligonucleotide frequencies independently of sequence homology. After an initial scan of HGT events using BLSOM, we identified 3.8% of the tsetse fly genome as HGT candidates. The predicted donors of these HGT candidates included known symbionts, such asWolbachia, as well as bacteria that have not previously been associated with the tsetse fly. We detected HGT candidates from diverse bacteria such asBacillus and Flavobacteria, suggesting a past association between these taxa. Functional annotation revealed that the HGT candidates encoded loci in various functional pathways, such as metabolic and antibiotic biosynthesis pathways. These findings provide a basis for understanding the coevolutionary history of the tsetse fly and its microbes and establish the effectiveness of BLSOM for the detection of HGT events.


---

## Body

## 1. Introduction

Tsetse flies (Glossina spp.) are the primary vectors of trypanosome parasites; they cause human African trypanosomiasis (or sleeping sickness) and animal African trypanosomiasis (or nagana) in Sub-Saharan African countries. The flies harbour three maternally transmitted endosymbionts,Wigglesworthia glossinidia,Sodalis glossinidius, andWolbachia pipientis, which influence host physiology. For example,Wigglesworthia provides essential nutrients, such as vitamins, to the host [1, 2] and influences host immune maturation [3]. Although the precise role ofSodalis in the tsetse fly is not clear, it appears to influence various host properties, such as longevity and susceptibility to trypanosome infections [4–6]. In many arthropod species [7, 8],Wolbachia induces strong cytoplasmic incompatibility, which was also observed in the tsetse fly [9]. In addition to these common bacteria, a recent microbial population analysis using a deep-sequencing approach revealed other facultative microorganisms from diverse bacterial families in the guts of tsetse flies, though their relative abundances were very low compared to that of the symbiontWigglesworthia [10].In addition to the parasitism of bacterial organisms themselves, partial genome sequences ofWolbachia are incorporated into the tsetse fly genome. Initially, Doudoumis et al. reported the incorporation of short fragments of threeWolbachia genes (16S rDNA,fbpA, andwsp) in the genomes of laboratory and naturalGlossina morsitans morsitans(Gmm) populations [11]. Subsequently, a whole-genome sequencing project revealed large insertions of theWolbachia genome in theGmm genome via horizontal gene transfer (HGT) events [12, 13]. These insertions were identified by extractingWolbachia-specific sequences from whole-genome Sanger sequencing reads and pyrosequencing data based on nucleotide homology with the complete genome sequences of threeWolbachia strains (wMel, wRi and wBm) [13]. Fluorescentin situhybridisation analyses further confirmed the presence of these insertions inGmm on the two sex chromosomes (X and Y) and the supernumerary B-chromosome [12, 13].HGT elements can be detected by two main methods: phylogeny-based and composition-based methods [14]. The first method relies on sequence alignments; HGT is identified when the position of a query sequence in a tree does not match that of a reference phylogeny. Although this approach is robust, the frequency of HGT events may be underestimated, especially when there is a lack of information on donor sequences [15]. The second method relies on nucleotide compositional features, such as G+C content, nucleotide frequencies, or codon usage [16–18], and theoretically does not require sequence homology. Batch learning self-organising map (BLSOM) is an alignment-free clustering method that generates a map independently of the order in which data are input via a learning process [19, 20]. This method enables the clustering of genomic sequence fragments based on the similarity of oligonucleotide frequencies, without any other taxonomical information; it has been successfully applied in genomic and metagenomic studies [20–22].The objective of this study was to characterise HGT from microorganisms in the genome ofGmmusing the alignment-free clustering method BLSOM. Using BLSOM, we detected a number of HGT candidates from diverse origins. In a comparison of the results for HGT fromWolbachia between methods, there was a high level of agreement between BLSOM and BLASTn, a homology-based approach. Based on functional annotation, these potential HGT elements encoded loci in various functional pathways.


## 2. Materials and Methods

### 2.1. Genome Sequences

The tsetse fly (Gmm) genome (Accession number CCAG010000000) and all prokaryotic sequences identified to the species level (n = 5,600) were obtained from GenBank (http://www.ncbi.nlm.nih.gov/Genbank/). When the number of undetermined nucleotides (Ns) exceeded 10% of the window size (5 kb), the sequence was omitted from the analysis. When the number of Ns was less than 10%, the oligonucleotide frequencies were normalised to the length without Ns and included in the analysis.


### 2.2. Batch Learning Self-Organising Map

G+C% is a fundamental value for the phylogenetic classification of microbial genomes, including viral genomes, but it cannot differentiate a wide variety of genomes. Oligonucleotide composition can distinguish species, even those with the same G+C%, because it varies substantially among genomes; accordingly, it is referred to as a “genome signature” [23]. Multivariate analyses, such as factor correspondence analysis and principal component analysis (PCA), are useful to investigate variation in gene sequences [24]. However, the clustering power of conventional multivariate analyses is inadequate when massive quantities of sequence data from a wide variety of genomes are analysed collectively. Kohonen’s self-organising map (SOM) is a powerful tool for clustering and visualising high-dimensional data vectors on a two-dimensional plane [25, 26]. To handle codon and oligonucleotide composition as high-dimensional data vectors, we modified the conventional SOM to develop the BLSOM [19, 20], which is suitable for genome sequence analyses and high-performance parallel computing. The initial weight vectors were defined by PCA, instead of random values, based on the finding that PCA can classify gene sequences into groups of known biological categories. Weight vectors (wij) were arranged in the two-dimensional lattice denoted byi (=0, 1,…,I−1) andj (=0, 1,…,J−1). Weight vectors (wij) were set and updated as described previously [19, 27]. A BLSOM program suitable for PC cluster systems is available on our website (http://bioinfo.ie.niigata-u.ac.jp/?BLSOM).


### 2.3. Detection of HGT Candidates in the Tsetse Fly Genome and Prediction of Their Origins Using BLSOM

To identify HGT candidates in the tsetse fly genome derived from prokaryotes, two types of large-scale BLSOM were used, that is, Tsetse+Prokaryotes- and Genus-BLSOM, using all genome sequences deposited in DDBJ/ENA/GenBank. A Tsetse+Prokaryotes-BLSOM was constructed with a degenerate tetranucleotide composition for all 5 kb sequences derived from tsetse fly genome sequences of longer than 5 kb plus 5,600 identified prokaryotes for which at least 10 kb of sequence was available from DDBJ/ENA/GenBank. The degenerate tetranucleotide composition was the composition of degenerate sets in which a pair of complementary tetranucleotides was added (e.g., ATGC and GCAT). To obtain more detailed phylotype information for the prokaryotic sequences, Genus-BLSOM was constructed for each phylum derived from 5,600 identified prokaryotes.For tsetse fly contigs of longer than 5 kb (9,710 contigs), a 5 kb window with a 1 kb step was used to obtain 303,250 segments (Figure1, Step 1), which were mapped to Tsetse+Prokaryotes-BLSOM by identifying the lattice point with the minimum Euclidian distances in the multidimensional space (Figure 1, Step 2). For every lattice point at which tsetse fly genomic segments were mapped to prokaryotic territories, the most abundant phylum was identified, and the mapped tsetse fly genomic segments were tentatively assumed to belong to the phylum. Finally, when the most abundant phylum in more than 40% of the segments derived from a single tsetse fly contig was the same, the tsetse fly contig was assigned to this phylum by BLSOM. To identify the phylogenetic origin of the tsetse fly genomic segments that were mapped to the prokaryotic territories on Tsetse+Prokaryotes-BLSOM, they were successively mapped on Genus-BLSOM (Figure 1, Step 3). Similar stepwise mappings of tsetse fly genomic segments on BLSOMs constructed with sequences from more detailed phylogenetic categories (e.g., genera) were conducted.Figure 1
Workflow for data processing and BLSOM analysis.


### 2.4. Detection of HGT Candidates by BLASTn

To detect HGT candidates derived fromWolbachia, a local BLASTn search was conducted with the contig sequences of the tsetse fly that exceeded 5 kb against the NCBI nonredundant nucleotide database. When more than 1 kb of the sequence showed similarity withWolbachia with top-hits and anE-value threshold of 1 × 10−5, the contigs were considered HGT candidates derived fromWolbachia.


### 2.5. Functional Classification of HGT Candidates

The sequence fragments with prokaryote origins were functionally annotated using KEGG (Kyoto Encyclopedia of Genes and Genomes) mapping [28] with the KAAS web server (http://www.genome.jp/tools/kaas/) [29]. KEGG Orthology (KO) assignments were obtained using the single-directional best-hit method. The organisms included in the analysis were as follows (based on IDs): hsa, dme, ath, sce, pfa, eco, sty, hin, pae, nme, hpy, rpr, mlo, bsu, sau, lla, spn, cac, mge, mtu, ctr, bbu, syn, aae, mja, afu, pho, and ape. The organisms in the database are listed on the KAAS web server (http://www.genome.jp/kaas-bin/kaas_org).


## 2.1. Genome Sequences

The tsetse fly (Gmm) genome (Accession number CCAG010000000) and all prokaryotic sequences identified to the species level (n = 5,600) were obtained from GenBank (http://www.ncbi.nlm.nih.gov/Genbank/). When the number of undetermined nucleotides (Ns) exceeded 10% of the window size (5 kb), the sequence was omitted from the analysis. When the number of Ns was less than 10%, the oligonucleotide frequencies were normalised to the length without Ns and included in the analysis.


## 2.2. Batch Learning Self-Organising Map

G+C% is a fundamental value for the phylogenetic classification of microbial genomes, including viral genomes, but it cannot differentiate a wide variety of genomes. Oligonucleotide composition can distinguish species, even those with the same G+C%, because it varies substantially among genomes; accordingly, it is referred to as a “genome signature” [23]. Multivariate analyses, such as factor correspondence analysis and principal component analysis (PCA), are useful to investigate variation in gene sequences [24]. However, the clustering power of conventional multivariate analyses is inadequate when massive quantities of sequence data from a wide variety of genomes are analysed collectively. Kohonen’s self-organising map (SOM) is a powerful tool for clustering and visualising high-dimensional data vectors on a two-dimensional plane [25, 26]. To handle codon and oligonucleotide composition as high-dimensional data vectors, we modified the conventional SOM to develop the BLSOM [19, 20], which is suitable for genome sequence analyses and high-performance parallel computing. The initial weight vectors were defined by PCA, instead of random values, based on the finding that PCA can classify gene sequences into groups of known biological categories. Weight vectors (wij) were arranged in the two-dimensional lattice denoted byi (=0, 1,…,I−1) andj (=0, 1,…,J−1). Weight vectors (wij) were set and updated as described previously [19, 27]. A BLSOM program suitable for PC cluster systems is available on our website (http://bioinfo.ie.niigata-u.ac.jp/?BLSOM).


## 2.3. Detection of HGT Candidates in the Tsetse Fly Genome and Prediction of Their Origins Using BLSOM

To identify HGT candidates in the tsetse fly genome derived from prokaryotes, two types of large-scale BLSOM were used, that is, Tsetse+Prokaryotes- and Genus-BLSOM, using all genome sequences deposited in DDBJ/ENA/GenBank. A Tsetse+Prokaryotes-BLSOM was constructed with a degenerate tetranucleotide composition for all 5 kb sequences derived from tsetse fly genome sequences of longer than 5 kb plus 5,600 identified prokaryotes for which at least 10 kb of sequence was available from DDBJ/ENA/GenBank. The degenerate tetranucleotide composition was the composition of degenerate sets in which a pair of complementary tetranucleotides was added (e.g., ATGC and GCAT). To obtain more detailed phylotype information for the prokaryotic sequences, Genus-BLSOM was constructed for each phylum derived from 5,600 identified prokaryotes.For tsetse fly contigs of longer than 5 kb (9,710 contigs), a 5 kb window with a 1 kb step was used to obtain 303,250 segments (Figure1, Step 1), which were mapped to Tsetse+Prokaryotes-BLSOM by identifying the lattice point with the minimum Euclidian distances in the multidimensional space (Figure 1, Step 2). For every lattice point at which tsetse fly genomic segments were mapped to prokaryotic territories, the most abundant phylum was identified, and the mapped tsetse fly genomic segments were tentatively assumed to belong to the phylum. Finally, when the most abundant phylum in more than 40% of the segments derived from a single tsetse fly contig was the same, the tsetse fly contig was assigned to this phylum by BLSOM. To identify the phylogenetic origin of the tsetse fly genomic segments that were mapped to the prokaryotic territories on Tsetse+Prokaryotes-BLSOM, they were successively mapped on Genus-BLSOM (Figure 1, Step 3). Similar stepwise mappings of tsetse fly genomic segments on BLSOMs constructed with sequences from more detailed phylogenetic categories (e.g., genera) were conducted.Figure 1
Workflow for data processing and BLSOM analysis.


## 2.4. Detection of HGT Candidates by BLASTn

To detect HGT candidates derived fromWolbachia, a local BLASTn search was conducted with the contig sequences of the tsetse fly that exceeded 5 kb against the NCBI nonredundant nucleotide database. When more than 1 kb of the sequence showed similarity withWolbachia with top-hits and anE-value threshold of 1 × 10−5, the contigs were considered HGT candidates derived fromWolbachia.


## 2.5. Functional Classification of HGT Candidates

The sequence fragments with prokaryote origins were functionally annotated using KEGG (Kyoto Encyclopedia of Genes and Genomes) mapping [28] with the KAAS web server (http://www.genome.jp/tools/kaas/) [29]. KEGG Orthology (KO) assignments were obtained using the single-directional best-hit method. The organisms included in the analysis were as follows (based on IDs): hsa, dme, ath, sce, pfa, eco, sty, hin, pae, nme, hpy, rpr, mlo, bsu, sau, lla, spn, cac, mge, mtu, ctr, bbu, syn, aae, mja, afu, pho, and ape. The organisms in the database are listed on the KAAS web server (http://www.genome.jp/kaas-bin/kaas_org).


## 3. Results and Discussion

### 3.1. Detection of HGT Candidates Using BLSOM

Of 303,250 sequence segments obtained from the tsetse genome, we found that 11,524 sequences (3.8%) clustered with reads from prokaryotes and thus were HGT candidates according to Tsetse+Prokaryotes-BLSOM (Figure2). These sequences were distributed across 2,960 different contigs, corresponding to 30.48% of all contigs. We assigned the most sequences to the phylum Firmicutes (n = 758), followed by the phyla Bacteroidetes (n = 370), Alphaproteobacteria (n = 90), and Gammaproteobacteria (n = 23) (Table 1). We did not assign 1,671 contigs to phyla owing in part to the presence of HGTs from multiple phyla within the same contig. It is also possible that these candidates were introduced by ancient HGT events, and their oligonucleotide compositions drifted over time [30], limiting the use of composition-based methods for classification.Table 1
Origins of HGT candidates at the phylum level.
Phylum Number of contigs Actinobacteria 10 Alphaproteobacteria 90 Aquificae 2 Bacteroidetes 370 Betaproteobacteria 4 Crenarchaeota 1 Cyanobacteria 1 Epsilonproteobacteria 16 Euryarchaeota 12 Firmicutes 758 Fusobacteria 1 Gammaproteobacteria 23 Spirochetes 1 Unassigned 1,671  Total 2,960 Prediction was obtained using Tsetse+Prokaryotes-BLSOM.Figure 2
Tsetse+Prokaryotes-BLSOM. (a) BLSOM using the degenerate tetranucleotide set for the tsetse fly plus 5,600 identified prokaryotes. Lattice points that include the sequences from the tsetse fly are indicated in green, those that contain no genomic sequences are indicated in white, and those containing sequences from a prokaryote are indicated in blue. Lattice points that include both tsetse fly- and prokaryote-sequences are shown in black. (b) Distribution of tsetse fly genome sequences. Only green lattice points are shown.
(a) (b)We performed further characterisation of donor sequences to the genus level using Genus-BLSOM (Figure3). The results of four dominant phyla are summarised in Table 2. We assigned the most sequences to the genusBacillus (n = 239). The second most highly represented origin was the class Flavobacteria (n = 187), which we were unable to classify to the genus level owing to coclustering with genome sequences that lacked genus information, followed by the generaStaphylococcus (n = 134),Enterococcus (n = 83),Wolbachia (n = 56),Polaribacter (n = 37), andListeria (n = 31). The HGT candidates associated with the generaWigglesworthia andSodalis, which are common endosymbionts of the tsetse fly, were not detected using BLSOM. Most tsetse flies are heavily infected withWigglesworthia; its abundance reaches over 99% in naturalGmm populations [10]. The lack of genome sequences associated withWigglesworthia may support the high quality of tsetse fly genome sequences, since some genome data are contaminated by symbiont genomes, which can lead to the false-positive detection of HGT events [31]. This result also suggests that the symbiosis betweenWigglesworthia and the tsetse fly was recent, as suggested by its genome features [32]. Nonetheless, we cannot exclude the possibility of bacterial genome contaminations in the tsetse fly genome since diverse bacteria exist in tsetse fly [10] and their sequences might not have been completely removed during the genome assembly process.Table 2
Origins of HGT candidates at the genus level.
Phylum Genus1 Number of contigs Alphaproteobacteria Anaplasma 4 Bartonella 8 Ehrlichia 1 Neorickettsia 8 Rickettsia 5 Wolbachia 56 Unassigned 8  Bacteroidetes Bacteroides 4 Cytophaga 1 Dyadobacter 1 Flavobacteria2 187 Flavobacterium 12 Kordia 18 Leadbetterella 3 Mucilaginibacter 1 Paludibacter 9 Pedobacter 1 Polaribacter 37 Prevotella 10 Psychroflexus 4 Spirosoma 4 Unassigned 78  Firmicutes Clostridium 4 Bacillus 239 Enterococcus 83 Epulopiscium 4 Erysipelothrix 1 Geobacillus 3 Lactobacillus 4 Lactococcus 1 Leuconostoc 17 Listeria 31 Lysinibacillus 7 Oenococcus 1 Paenibacillus 1 Peptoniphilus 1 Staphylococcus 134 Streptococcus 9 Thermoanaerobacter 3 Turicibacter 2 Veillonella 3 Unassigned 210  Gammaproteobacteria Acinetobacter 4 Enterobacter 6 Escherichia 1 Haemophilus 1 Shewanella 1 Vibrio 1 Xylella 1 Unassigned 8 1Prediction was obtained using Genus-BLSOM. 2Classification to the class level was obtained.Figure 3
Genus-BLSOM. (a) Alphaproteobacteria. (b) Bacteroidetes. (c) Epsilonproteobacteria. (d) Euryarchaeota. (e) Firmicutes. (f) Gammaproteobacteria. Lattice points that include sequences from more than one genus are indicated in black, those including no sequences are indicated in white, and those including sequences from a single genus are indicated in individual color.
(a) (b) (c) (d) (e) (f)The high frequency of HGT candidates from the genusBacillus suggests that there was a strong association between the tsetse fly andBacillus in the past. Members of the genusBacillus are ubiquitous in nature and have been isolated from diverse environments such as water, soil, plants, animals, and air [33]. Some species, such asBacillus thuringiensis, have been well studied as agents of biological control of arthropods [34]. Kaaya and Darji infected severalBacillusspecies, includingB. thuringiensis, to the adultGmm and found that the mortality ofGmm was depending on the bacterial species [35], indicating that someBacillus species may have infectedGmm persistently without adverse effect on the hosts and served as HGT donors. In fact, in a microbiota analysis of one tsetse fly species,Glossina fuscipes fuscipes, the bacteria belonging to the genusBacillus were found dominant in a culture-dependent manner [36].In contrast, there is no report on the relationship between the tsetse fly and Flavobacteria, which was identified as a second dominant donor of HGT candidates in this study (Table2). Flavobacteria are symbionts in several arthropods [37–43], which indicates a high probability of the proliferation of this group of bacteria in arthropod hosts including tsetse fly. A comparative genome analysis of a flavobacterial symbiont (Blattabacterium strain Bge) in the omnivorous German cockroach (Blattella germanica) suggested that it plays roles in nutrient supply to the host, amino acid catabolism, and nitrogen excretion [41]. Flavobacterial symbionts in the ladybird (Coleomegilla maculata), and coccinellid beetle (Adonia variegata) induce male-killing [38, 39], in which male progeny in infected females die during embryogenesis. This phenomenon is widely recognised in other bacteria, such asWolbachia,Rickettsia,Arsenophonus,Spiroplasma, andCardinium [44]. Hurst et al. proposed that two male-killing symbionts cannot coexist at equilibrium in a single host species based on an observational study of the two-spot ladybird (Adalia bipunctata) infected with two symbionts,Rickettsia andSpiroplasma [39]. The presence ofWolbachia in the tsetse fly and HGT elements fromWolbachia in the tsetse fly genome [12, 13] may explain the absence of Flavobacteria in current tsetse fly populations.


### 3.2. HGT Candidates Derived fromWolbachia

We performed a BLASTn analysis to detect HGT candidates derived fromWolbachia. For 38 contigs, we detected sequence homology withWolbachia sequences based on the criteria described earlier. Of these 38 contigs, we identified 36 as HGT candidates fromWolbachia using BLSOM, while we assigned the remaining two contigs to the genusRickettsia using BLSOM. Accordingly, we observed a high level of agreement (36/38 = 94.7%) between the two approaches, indicating that the analytical sensitivity of BLSOM is at least comparable to that of BLASTn. These results also suggest that the genetic elements derived fromWolbachia were recently introduced to the tsetse fly genome. Based on a comparison between the BLSOM and BLASTn results, we observed 20 contigs that were only identified as HGT candidates using BLSOM. There may be as-yet unidentifiedWolbachia strains, explaining the failure to identify HGT candidates based on sequence homology. In fact, for short regions (i.e., 338, 492, 497, and 497 bp) of four contigs identified as HGT candidates fromWolbachia only using BLSOM, we detected sequence homology with aWolbachia endosymbiont ofGmm (Accession number AWUH01000121) with anE-value of <1 × 10−5.Genetic elements related toWolbachia that were presumably obtained by HGT have been detected in the genomes of multiple arthropod species, including the adzuki bean beetle (Callosobruchus chinensis) [45, 46], a fruit fly (Drosophila ananassae) [47], parasitoid wasps (Nasonia spp.) [47], the pea aphid (Acyrthosiphon pisum) [48], two mosquito species (Aedes aegypti andAedes mascarensis) [49, 50], and the longicorn beetle (Monochamus alternatus) [51]. Some HGT events could be explained by nuclear-phage recombination, as proposed previously [49], but further studies are needed to determine the specific mechanisms of transfer. Since BLSOM can be used to detectWolbachia-derived HGTs that exhibit low sequence similarity with knownWolbachia strains, it provides an alternative method with which to explore the mechanisms of HGT. The application of BLSOM to an increasing number of eukaryotic genomes will reveal the diversity and frequency ofWolbachia-derived HGTs in other arthropods, including vectors of medical and veterinary importance.


### 3.3. Functional Classification of HGT Candidates

We mapped all of the HGT candidates identified using Tsetse+Prokaryotes-BLSOM (n = 11,524) to the KEGG pathway. The KO included 317 biological pathways. The predicted pathways were mainly related to “Metabolic pathways” (193 molecules), “Biosynthesis of secondary metabolites” (75 molecules), and “Biosynthesis of antibiotics” (50 molecules). These results suggested that the HGT candidates have the potential to affect a large number of metabolic activities; however, further analyses are essential to demonstrate the active transcription of HGT-acquired genes using transcriptomics or gene-specific reverse transcription-PCR; such analyses can provide initial evidence for the functional importance of HGT-acquired genes [52]. In general, genes transferred to host genomes are pseudogenised via the acquisition of mutations, including insertions and deletions [46, 52, 53]. Unfortunately, since we could not employ RNA sequencing data into our analysis, it is not clear to what extent the detected HGT candidates have been pseudogenised. Nonetheless, active transcription of HGT-acquired genes has been detected in recipient hosts, such as aWolbachia-derived gene in theAedes albopictus C6/36 cell line [54]. Moreover, an increasing number of studies suggests that HGT-acquired genes facilitate the establishment of obligate mutualistic relationships between arthropods and their symbionts [55]. Analyses of the functional roles of HGT-acquired genes may improve our understanding of the complex interactions between the tsetse fly, microbes, and pathogens.


## 3.1. Detection of HGT Candidates Using BLSOM

Of 303,250 sequence segments obtained from the tsetse genome, we found that 11,524 sequences (3.8%) clustered with reads from prokaryotes and thus were HGT candidates according to Tsetse+Prokaryotes-BLSOM (Figure2). These sequences were distributed across 2,960 different contigs, corresponding to 30.48% of all contigs. We assigned the most sequences to the phylum Firmicutes (n = 758), followed by the phyla Bacteroidetes (n = 370), Alphaproteobacteria (n = 90), and Gammaproteobacteria (n = 23) (Table 1). We did not assign 1,671 contigs to phyla owing in part to the presence of HGTs from multiple phyla within the same contig. It is also possible that these candidates were introduced by ancient HGT events, and their oligonucleotide compositions drifted over time [30], limiting the use of composition-based methods for classification.Table 1
Origins of HGT candidates at the phylum level.
Phylum Number of contigs Actinobacteria 10 Alphaproteobacteria 90 Aquificae 2 Bacteroidetes 370 Betaproteobacteria 4 Crenarchaeota 1 Cyanobacteria 1 Epsilonproteobacteria 16 Euryarchaeota 12 Firmicutes 758 Fusobacteria 1 Gammaproteobacteria 23 Spirochetes 1 Unassigned 1,671  Total 2,960 Prediction was obtained using Tsetse+Prokaryotes-BLSOM.Figure 2
Tsetse+Prokaryotes-BLSOM. (a) BLSOM using the degenerate tetranucleotide set for the tsetse fly plus 5,600 identified prokaryotes. Lattice points that include the sequences from the tsetse fly are indicated in green, those that contain no genomic sequences are indicated in white, and those containing sequences from a prokaryote are indicated in blue. Lattice points that include both tsetse fly- and prokaryote-sequences are shown in black. (b) Distribution of tsetse fly genome sequences. Only green lattice points are shown.
(a) (b)We performed further characterisation of donor sequences to the genus level using Genus-BLSOM (Figure3). The results of four dominant phyla are summarised in Table 2. We assigned the most sequences to the genusBacillus (n = 239). The second most highly represented origin was the class Flavobacteria (n = 187), which we were unable to classify to the genus level owing to coclustering with genome sequences that lacked genus information, followed by the generaStaphylococcus (n = 134),Enterococcus (n = 83),Wolbachia (n = 56),Polaribacter (n = 37), andListeria (n = 31). The HGT candidates associated with the generaWigglesworthia andSodalis, which are common endosymbionts of the tsetse fly, were not detected using BLSOM. Most tsetse flies are heavily infected withWigglesworthia; its abundance reaches over 99% in naturalGmm populations [10]. The lack of genome sequences associated withWigglesworthia may support the high quality of tsetse fly genome sequences, since some genome data are contaminated by symbiont genomes, which can lead to the false-positive detection of HGT events [31]. This result also suggests that the symbiosis betweenWigglesworthia and the tsetse fly was recent, as suggested by its genome features [32]. Nonetheless, we cannot exclude the possibility of bacterial genome contaminations in the tsetse fly genome since diverse bacteria exist in tsetse fly [10] and their sequences might not have been completely removed during the genome assembly process.Table 2
Origins of HGT candidates at the genus level.
Phylum Genus1 Number of contigs Alphaproteobacteria Anaplasma 4 Bartonella 8 Ehrlichia 1 Neorickettsia 8 Rickettsia 5 Wolbachia 56 Unassigned 8  Bacteroidetes Bacteroides 4 Cytophaga 1 Dyadobacter 1 Flavobacteria2 187 Flavobacterium 12 Kordia 18 Leadbetterella 3 Mucilaginibacter 1 Paludibacter 9 Pedobacter 1 Polaribacter 37 Prevotella 10 Psychroflexus 4 Spirosoma 4 Unassigned 78  Firmicutes Clostridium 4 Bacillus 239 Enterococcus 83 Epulopiscium 4 Erysipelothrix 1 Geobacillus 3 Lactobacillus 4 Lactococcus 1 Leuconostoc 17 Listeria 31 Lysinibacillus 7 Oenococcus 1 Paenibacillus 1 Peptoniphilus 1 Staphylococcus 134 Streptococcus 9 Thermoanaerobacter 3 Turicibacter 2 Veillonella 3 Unassigned 210  Gammaproteobacteria Acinetobacter 4 Enterobacter 6 Escherichia 1 Haemophilus 1 Shewanella 1 Vibrio 1 Xylella 1 Unassigned 8 1Prediction was obtained using Genus-BLSOM. 2Classification to the class level was obtained.Figure 3
Genus-BLSOM. (a) Alphaproteobacteria. (b) Bacteroidetes. (c) Epsilonproteobacteria. (d) Euryarchaeota. (e) Firmicutes. (f) Gammaproteobacteria. Lattice points that include sequences from more than one genus are indicated in black, those including no sequences are indicated in white, and those including sequences from a single genus are indicated in individual color.
(a) (b) (c) (d) (e) (f)The high frequency of HGT candidates from the genusBacillus suggests that there was a strong association between the tsetse fly andBacillus in the past. Members of the genusBacillus are ubiquitous in nature and have been isolated from diverse environments such as water, soil, plants, animals, and air [33]. Some species, such asBacillus thuringiensis, have been well studied as agents of biological control of arthropods [34]. Kaaya and Darji infected severalBacillusspecies, includingB. thuringiensis, to the adultGmm and found that the mortality ofGmm was depending on the bacterial species [35], indicating that someBacillus species may have infectedGmm persistently without adverse effect on the hosts and served as HGT donors. In fact, in a microbiota analysis of one tsetse fly species,Glossina fuscipes fuscipes, the bacteria belonging to the genusBacillus were found dominant in a culture-dependent manner [36].In contrast, there is no report on the relationship between the tsetse fly and Flavobacteria, which was identified as a second dominant donor of HGT candidates in this study (Table2). Flavobacteria are symbionts in several arthropods [37–43], which indicates a high probability of the proliferation of this group of bacteria in arthropod hosts including tsetse fly. A comparative genome analysis of a flavobacterial symbiont (Blattabacterium strain Bge) in the omnivorous German cockroach (Blattella germanica) suggested that it plays roles in nutrient supply to the host, amino acid catabolism, and nitrogen excretion [41]. Flavobacterial symbionts in the ladybird (Coleomegilla maculata), and coccinellid beetle (Adonia variegata) induce male-killing [38, 39], in which male progeny in infected females die during embryogenesis. This phenomenon is widely recognised in other bacteria, such asWolbachia,Rickettsia,Arsenophonus,Spiroplasma, andCardinium [44]. Hurst et al. proposed that two male-killing symbionts cannot coexist at equilibrium in a single host species based on an observational study of the two-spot ladybird (Adalia bipunctata) infected with two symbionts,Rickettsia andSpiroplasma [39]. The presence ofWolbachia in the tsetse fly and HGT elements fromWolbachia in the tsetse fly genome [12, 13] may explain the absence of Flavobacteria in current tsetse fly populations.


## 3.2. HGT Candidates Derived fromWolbachia

We performed a BLASTn analysis to detect HGT candidates derived fromWolbachia. For 38 contigs, we detected sequence homology withWolbachia sequences based on the criteria described earlier. Of these 38 contigs, we identified 36 as HGT candidates fromWolbachia using BLSOM, while we assigned the remaining two contigs to the genusRickettsia using BLSOM. Accordingly, we observed a high level of agreement (36/38 = 94.7%) between the two approaches, indicating that the analytical sensitivity of BLSOM is at least comparable to that of BLASTn. These results also suggest that the genetic elements derived fromWolbachia were recently introduced to the tsetse fly genome. Based on a comparison between the BLSOM and BLASTn results, we observed 20 contigs that were only identified as HGT candidates using BLSOM. There may be as-yet unidentifiedWolbachia strains, explaining the failure to identify HGT candidates based on sequence homology. In fact, for short regions (i.e., 338, 492, 497, and 497 bp) of four contigs identified as HGT candidates fromWolbachia only using BLSOM, we detected sequence homology with aWolbachia endosymbiont ofGmm (Accession number AWUH01000121) with anE-value of <1 × 10−5.Genetic elements related toWolbachia that were presumably obtained by HGT have been detected in the genomes of multiple arthropod species, including the adzuki bean beetle (Callosobruchus chinensis) [45, 46], a fruit fly (Drosophila ananassae) [47], parasitoid wasps (Nasonia spp.) [47], the pea aphid (Acyrthosiphon pisum) [48], two mosquito species (Aedes aegypti andAedes mascarensis) [49, 50], and the longicorn beetle (Monochamus alternatus) [51]. Some HGT events could be explained by nuclear-phage recombination, as proposed previously [49], but further studies are needed to determine the specific mechanisms of transfer. Since BLSOM can be used to detectWolbachia-derived HGTs that exhibit low sequence similarity with knownWolbachia strains, it provides an alternative method with which to explore the mechanisms of HGT. The application of BLSOM to an increasing number of eukaryotic genomes will reveal the diversity and frequency ofWolbachia-derived HGTs in other arthropods, including vectors of medical and veterinary importance.


## 3.3. Functional Classification of HGT Candidates

We mapped all of the HGT candidates identified using Tsetse+Prokaryotes-BLSOM (n = 11,524) to the KEGG pathway. The KO included 317 biological pathways. The predicted pathways were mainly related to “Metabolic pathways” (193 molecules), “Biosynthesis of secondary metabolites” (75 molecules), and “Biosynthesis of antibiotics” (50 molecules). These results suggested that the HGT candidates have the potential to affect a large number of metabolic activities; however, further analyses are essential to demonstrate the active transcription of HGT-acquired genes using transcriptomics or gene-specific reverse transcription-PCR; such analyses can provide initial evidence for the functional importance of HGT-acquired genes [52]. In general, genes transferred to host genomes are pseudogenised via the acquisition of mutations, including insertions and deletions [46, 52, 53]. Unfortunately, since we could not employ RNA sequencing data into our analysis, it is not clear to what extent the detected HGT candidates have been pseudogenised. Nonetheless, active transcription of HGT-acquired genes has been detected in recipient hosts, such as aWolbachia-derived gene in theAedes albopictus C6/36 cell line [54]. Moreover, an increasing number of studies suggests that HGT-acquired genes facilitate the establishment of obligate mutualistic relationships between arthropods and their symbionts [55]. Analyses of the functional roles of HGT-acquired genes may improve our understanding of the complex interactions between the tsetse fly, microbes, and pathogens.


## 4. Conclusions

We investigated the use of BLSOM to detect HGT candidates in the tsetse fly genome. Using BLSOM, we successfully detected a number of HGT candidates from diverse bacterial origins. The HGT candidates represented 3.8% of the tsetse fly genome. The predicted donors of these HGT elements includedWolbachia, a well-known symbiont of the tsetse fly. In addition, using BLSOM, we identified HGT candidates from bacteria that have not previously been associated with the tsetse fly. We observed the HGT candidates from diverse bacteria such asBacillus and Flavobacteria, suggesting a strong past association between these taxa. In a comparison between BLASTn and BLSOM results for the detection of HGT candidates fromWolbachia, the analytical sensitivity of BLSOM was at least comparable to that of the sequence homology-based approach. Furthermore, BLSOM can be used to detect HGT elements from organisms with low similarity with currently available sequences. These data obtained using BLSOM provide a basis for understanding the coevolutionary history of the tsetse fly and its microbes.



---

*Source: 3164624-2016-12-15.xml*