Horizontal gene transfer in silkworm, Bombyx mori

Genome

The genome of Bombyx mori has a size of 428.7 Mb and is composed of 28 chromosomes [12]. Gene candidates for HGT in B. mori were identified by screening the 14,623 coding genes downloaded from SilkDB [13].

Local database generating

In order to construct the phylogenetic trees for each sequence identified with the target taxon/similarity profile automatically and preciously, a local database (300 cores Processor BladeSystem in Zhejiang University, IBM-Biocomputing Laboratory) containing predicted protein sequences (2,385,947 coding sequences in total) from 382 species representing a wide diversity of eukaryotes and prokaryote taxa including seven insects, seven plants, 43 fungi, 30 other eukaryotes and 295 prokaryotes (Additional file 1).

Search procedure

Each candidate sequence in B. mori was compared against sequences in the local database using BLASTp [14] and the highest similarity sequences from each species were extracted for further analysis (E-value cutoff 10^-20). We used the local database first to exclude insect unique genes and the genes that, when compared against all the genomes showed a higher BLASTp similarity score to Arthropoda than to any other species (Figure 1). All sequences that had a high sequence similarity to the transposon were also excluded from further analysis as we aimed to study the effect of HGT on the putative functional genes (Figure 1). After detecting the putative transferred genes, expression sequence tags (EST) information was used to remove out the unexpressed genes. The detailed procedure was showed in Figure 1.

Phylogenetic analysis

The remaining candidate HGT genes were selected to search against GenBank non-redundant protein database (nr). Search strategy was the same as described above. These sequences were aligned using ClustalW [15], the conserved region of each alignment was trimmed using Gblocks [16] with a stringent settings described previously [17]. Maximum Likelihood (ML) phylogenies were constructed by Phyml [18] using a JTT + Γ + I substitution model (Γ + I parameters were estimated by Phyml). The proportion of invariant sites was estimated from the data. For Bayesian phylogenies generated by MrBayes [19], two independent Metropolis-coupled Monte Carlo Markov Chain (MCMC) runs, each with one cold and three heated chains (heat parameter = default), were analyzed for one million generations after a burn-in of 25,000 samples and allowed mixed models of amino acid substitution. For Maximum Likelihood phylogenies, 1,000 bootstraps were performed to gain the branch support values. To be considered as an indicative of a potential HGT event, the ML branch support about B. mori should ≥ 80% while Bayesian posterior probability should ≥ 85%.

To test the support for contentious topology, we performed nonparametric branch support tests based on a Shimodaira-Hasegawa-like procedure using Tree-Puzzle [20]. In addition, we used this software to test the statistical significance (5%) of specific topology over a collapsed version of the same branching relationship.

Functional analysis

To test whether the transferred genes were associated with the domestication of B. mori, single nucleotide polymorphism (SNP) genelist [2] produced by resequencing 40 B. mori genomes was used to search if the candidate genes are on the list. Meanwhile, expression and coexpression analysis about the transferred genes were also conducted. A Pearson correlation coefficient was calculated based on microarray gene expression profiling, which were downloaded from SilkDB [13] to investigate the function of transferred genes analyzed in this paper. Background correction and data normalization was done by RMA algorithm in Bioconductor [21] and the poor annotated probe-sets were removed. Measurements for unique gene was calculated from means of the probes belong to same gene. The gene coexpressed with one of either gene in ten putative transferred genes with the absolute correlation |r| > 0.5 was further selected to analyze the function. For each transferred gene, we obtained the top 300 genes in each list sorted by absolute correlation value r referred to as coexpressed genes in descending order. The distribution of coexpressed genes with each transferred gene was mapped onto pathways. A score reflecting the extent of coexpression in each gene was assigned to each pathway by using following formula: , where G is the number of genes in a pathway present in the top 300 of all same transferred genes. R _i is the rank of the ith gene of the specified pathway in a list, r _i is the absolute correlation value r.

The process of HGT candidate search

In this study, stringent filters were used to identify the real HGT events. Although this procedure could identify the really HGT event undoubtedly, it may neglect several potential HGT events. First, the HGT in the common ancestor of Insecta may be considered as vertical inheritance in this study if none of the other sampled taxons has homolog. As no homologs were detected outside insects, it might be assumed that this HGT was from an extinct linkages [25]. However, this assumption can only been proved by fossil evidence. Second, our filter could also remove genes transferred from other animal to B. mori as we only retained the genes showed a higher BLAST similarity score to a gene from any other species rather than Arthropoda (Material and Methods). However, fewer instances of potential HGTs between eukaryotes are known [4, 26], and most of these HGTs were plastids-eukaryote transfer. Finally, we only retained the transferred genes that have putative functions as functional genes may play more important role in physical and biochemical reaction in B. mori.

Comparative phylogenetic analysis, which included models that account for site rate heterogeneity and, where appropriate, comparative topology tests were used in our study. We finally identified 10 genes that was the most consistent explanation to have been transferred from the genome of bacteria to the genome of B. mori. The number of transferred genes is nearly the same as recently Acyrthosiphon pisum sequencing paper [7]. The direction of transfer from prokaryote to B. mori is also highly suggested in our analysis. As in the phylogenetic trees, B. mori nested within a bacterial clade, and many bacterial species formed basel branches. Meanwhile, all of these 10 genes are intronless in the B. mori. Based on these results, the direction of transfer from prokaryote to B. mori is mainly confirmed. Recent findings suggested that substantial HGT between prokaryotic endosymbionts and their insect hosts were always happened [7, 27–32]. However, in this study, none of the putative transferred gene was acquired from a previous commonly reported endosymbiont (Wolbachia or Buchnera).

The HGT events in B. mori are rare but recent

It has long been described that many insects such as pea aphid, mosquitoes, and fruit flies have been influenced by HGT from Wolbachia[33]. Nikoh et al. have suggested that nearly 30% of a Wolbachia genome was found on the X-chromosome of the insect, probably as the result of a single HGT event [30]. In this study, we only get 10 HGTs in silkworm. Compared to previous reports, we may think that the HGT events are rare compared to other insects. The CDS of 10 putative transferred genes were used to search against wildsilkbase [34], which is a BLAST searchable catalogue of ESTs generated from three major wild silkworms, Antheraea assama, Samia cynthia ricini and Antheraea mylitta. We found that the transferred genes in B. mori were actually absent in three wild silkmoths. This phenomenon can be explained by two hypotheses. The first suggests that the aliened genome fragment was transferred to the B. mori chromosome after speciation of Bombycoidea, while the latter suggests an ancient HGT in the common ancestor with several losses. However, the "losses" hypothesis is highly unlikely for two reasons. It requires massive independent losses to explain the presence of these genes only in B. mori. Furthermore, this hypothesis has difficulty to explain the remarkable sequence identity between B. mori and its candidate donors (Table 1), while the general gene similarity between wild and domesticated silkworm is 58.6% based on the protein sequences of Antheraea assama deposited in Genbank. Taking into account, it appears most likely that the aliened genome fragment was transferred to the B. mori chromosome after speciation of Bombycoidea. The result suggested that these HGT events may be happened much recently. The recovery of only 10 potential HGT in B. mori suggests that HGT events happened in this species have played only a very minor role in its evolution (0.065% in the B. mori protein coding genes). Interestingly, none of the HGT events involved donor linkages from the common ancestor of Insecta. Instead, all of the HGTs involved transfer from bacteria or plant directly to the Bombycidae except chitinase, which has been demonstrated to be an ancient HGT event [9]. It is worth noting that half of these transfers originate from the Enterobacteriaceae species. Many Enterobacteriaceae species have been demonstrated to be as symbionts of insects [35]. It may facilitate to exchange the genetic material between symbionts and hosts.

Putative functional assignment of 10 HGTs

The potential identities and biological functions of the 10 HGT candidates were investigated based on the similarity of their putatively encoded products to known proteins. We found that the 10 putative transferred genes were all enzymes, which suggested a strong functional trend among these HGTs in B. mori. Furthermore, hypergeometric distribution analysis of the candidate HGT genes indicated catalytic activity and hydrolase activity genes are significantly overrepresented (with p-values 4.6e^-05 and 6.4e^-06, respectively). Compared with our study, none of the putative functional gene was identified in Aphid genome sequence analysis [7]. The acquisition of such genes by HGT may therefore have been selected by a competitive advantage in the recipient to degrade or take up substrates with greater efficiency.

One HGT data set (Figure 2) transferred from bacteria to B. mori in this study is a enzyme encoding NAD-dependent epimerase/dehydratase (BGIBMGA010285, pfam01370, e-value 3.2e^-25), which catalyses the last step in biosynthesis of GDP-d-rhamnose in bacteria [36]. This kind of gene has been demonstrated to be transferred in chromalveolates [37]. Expression analysis suggested that this gene was highly expressed in integument and fat body (Figure 3), whereas coexpression analysis revealed that this gene may be associated with biotin metabolism with high value support (Additional file 2). As insect store energy reserves mainly in fat body cells [38] and biotin metabolism is related to silk synthesis by nutrient affection [39]. Based on these results, we can infer that this gene may be associated with nutrient and energy metabolism in B. mori.

A putative glycerophosphoryl diester phosphodiesterase (BGIBMGA007767, pfam03009, e-value 1e^-52) was encoded by an HGT candidate from the Pseudomonas species to B. mori (Additional file 3). This enzyme display broad specificity for glycerophosphodiesters, glycerophosphocholine, glycerophosphoethanolamine, glycerophosphoglycerol and bis (glycerophosphoglycerol), all of which are hydrolysed by this enzyme. This gene also has been reported to be involved in HGT between spider and bacteria [40]. Expression analysis suggested that this gene was not only highly expressed in integument and fat body but also in midgut and hemocyte. Coexpression analysis indicated that the function of this gene may be consistent with BGIBMGA010285.

Two genes, encoded chitinase and sucrose-6-phosphate hydrolase, are the only two genes which have been reported to be transferred from bacteria to B. mori previously [8, 9] were also identified in this study. B. mori acquired the first gene from a bacterium or an ancestral baculovirus and the latter gene from a bacterium [8, 9]. The results suggested that both genes bestow a clear selective advantage to B. mori. The discovery of both genes indicated that our strategy in searching candidate HGTs in B. mori is reliable.

Another candidate HGT was a gene transferred from Edwardsiella tarda (Additional file 4). Homology analysis suggested that this gene was a glucose-1-phosphatase/inositol phosphatase (BGIBMGA011204, pfam00328, e-value 1.4e^-17). This enzyme has been reported to enhance the phytase activity in bacteria [41]. Expression pattern revealed that this gene is extremely highly expressed in head (Figure 2). Coexpression analysis suggested that the function of this gene may be most associated with Butirosin and neomycin biosynthesis (Additional file 2). As the gene in this pathway is involved in aminoglycoside antibiotics biosynthesis, which is effective against several types of bacteria [42]. The result indicated that this gene may be associated with the resistance of B. mori to bacterial pathogen.

The sixth HGT candidate (Additional file 5) showed sequence similarity to alginate lyase (BGIBMGA005615, pfam08240, e-value 2.1e^-08), which can catalyze the degradation of alginate by a β-elimination mechanism [43]. In bacteria, this enzyme has been reported to destroy the cell detachment, which is the key process in biofilm formation [44]. Unfortunately, there is no corresponding probe to detect this gene in microarray data.

The seventh, eighth and ninth HGT are putative kynureninase (BGIBMGA007146, pfam00266, e-value 2.6e^-18), N-methyltryptophan oxidase (BGIBMGA008215, pfam01266, e-value 4.2e^-54) and gamma- glutamyltranspeptidase (pfam01019, e-value 4.3e^-125), respectively. The seventh gene encoded by an HGT candidate from Listeria grayi to B. mori (Additional file 6). The mutation of this gene has been reported to be associated with abnormal body coloration in B. mori[45]. The eighth and ninth genes were both encoded by an HGT candidate from Enterobacteriaceae to B. mori (Additional file 7 and 8). In both genes, the eighth one has been reported to be involved in substrate specificity determination [46] while the ninth gene plays a key role in the gamma-glutamyl cycle, a pathway for the synthesis and degradation of glutathione and drug and xenobiotic detoxification [47]. All of the three genes are highly expressed in malpighian tubule (Figure 2), and the coexpression patterns are similar (Additional file 2). The results indicated that these genes may be associated with valine, leucine and isoleucine degradation pathway (Additional file 2). It has been reported that after B. mori was infected by its bacterial pathogen Bacillus bombyseptieus, the expression level of the genes which are associated with valine, leucine and isoleucine degradation pathway were more than three times the baseline levels in the malpighian tubule. Huang et al. (2009) speculated that these genes might be involved in the detoxification of malpighian tubules [48]. The results suggested that these genes may be associated with the resistance of B. mori to bacterial pathogen.

The tenth gene which encoded as aromatic ring-opening dioxygenase LigB subunit (BGIBMGA003842, pfam02900, e-value 1.2e^-51) was transferred from a plant donor (Additional file 9). No tissue specific expression pattern was observed in this gene (Figure 2). However, coexpression analysis suggested that this gene may be most associated with Butirosin and neomycin biosynthesis (Additional file 2). It has been reported that Butirosin and neomycin belong to a family of clinically valuable 2-deoxystreptamine (2DOS)-containing aminoglycoside antibiotics [49]. The result indicated that this gene may be associated with the resistance of B. mori to bacterial pathogen.

In this study, the coexpressed genes involved in HGT are mapped onto biological pathway to investigate the potential function in B. mori. This scheme provides a comprehensive way to identify previously unknown functional patterns in sets of genes with known functions. It has been well known that the domesticated B. mori has gained or lost some representative behavioral characteristics compared to its wild ancestor. A lot of researches have focused on the search for 'domestication genes' in mammals by the polygenic nature of most behavioral traits [50]. Indeed, the linking of coat-color genes in mammals to brain biochemistry and behaviors has been found to be able to facilitate domestication and ease of handling [50]. In our study, SNP data searching indicated that all of the 10 genes have no direct relationship with human domestication although all of these genes are suggested to be transferred recently after the divergence from the wild silkworm. Therefore, we could suggest that HGT in B. mori maybe mainly depend on its nature requirement, but not on human activity. However, our study strongly suggested that these transferred genes increase the ability of B. mori to natural survival and adaption.