Identification of horizontally transferred genes in the genus Colletotrichum reveals a steady tempo of bacterial to fungal gene transfer

Identification of HGT candidates

The most robust method to detect HGT is phylogenetic analysis [23]. Therefore, to detect putative HGT events in three Colletotrichum species (C. graminicola M1.001 [30], C. higginsianum IMI 349063 [30] and C. gloeosporioides Cg-14 [32]), we developed a pipeline that consists of a series of BLAST searches and automated filtering designed to reduce the number of unlikely HGT candidates, followed by manual evaluation of phylogenetic trees (Additional file 1: Figure S1). Since our pipeline includes several steps of manual tree inspection, we included several filters designed to reduce the number of candidates that require manual inspection and to limit them to those with the highest sequence similarity. The first BLAST search was performed using a database of proteins from organisms with complete proteome available in UniProt (http://www.uniprot.org). Next, we selected proteins that having at least 80% of the top 120 hits (e-value e-5) with a taxonomic classification other than fungi as candidates for further analysis. The threshold of 80% was selected by evaluating previously described HGT candidates reported by Richards et al. [12], Schmitt and Lumbsch [21] and Richards [28] (see Methods for details). Next, we subjected the HGT candidates to three phylogenetic analyses using different sets of homologous sequences in each phylogeny. The first phylogeny was constructed with homologous sequences from the UniProt complete proteome database (476 proteins from the 3 species were selected). For the second phylogeny we performed a BLAST search of the GenBank nr database and included the 20 best hits from each kingdom (Archaea, Bacteria and Eukaryota) to avoid a possible under or overrepresentation for the abundance of sequences from any one kingdom in the BLAST results. This procedure was performed to observe events of inter-kingdom HGT and to observe the place into the tree of sequences with the same taxonomic label of the query but that were not the best hits in the BLAST search and therefore excluded from the first analysis that uses the best 120 hits. The third phylogeny was constructed with the best 100 BLAST hits from the nr database. The phylogenetic trees were evaluated manually, selecting only those that have well-supported topologies that are clearly incongruent with known species relationships among the taxa. Asymmetric or ladder-shape trees were excluded because such tree topologies are often a signal of long branch attraction or lack of phylogenic information [29]. Additionally, candidates with few homologues or with low sequence similarity to all of their BLAST hits were also excluded. In cases where only hits from two kingdoms were obtained (i.e. Bacteria and Eukaryota), high sequence similarity (a minimum of 30% pairwise similarity) and coverage (over 80% coverage) were required to consider them as candidates. The BLAST searches and tree evaluations were performed serially rather than in parallel to minimize the number of manual phylogenetic tree evaluations required (Additional file 1: Figure S1).

None of the 11 C. graminicola HGT candidates of bacterial origin have introns consistent with their prokaryotic origin. Only 26.4% of the genes encoded in the C. graminicola genome lack introns and the probability of selecting 11 intronless genes at random from the C. graminicola genome is 4.06e^-7 further supporting the hypothesis that the 11 candidates are of bacterial origin. To corroborate this calculate we performed a Mann–Whitney-Wilcoxon test to determine whether the two samples of genes (the whole genome and the HGT candidates) are the same, with respect to the number of introns in the genes. These two set of genes were significantly different (p ≈ 0), from which we conclude that the 11 HGT candidates have a different intron distribution than the genome.

We considered that the HGT candidates might lack introns because they are members of gene families that typically lack introns. To test this alternative hypothesis, we first identified 75 homologs to the HGT candidates by searching the proteomes for proteins with the same functional annotation as the candidates. We counted the number of introns in the 75 homologs and compared this sample to the intron content of the rest of the genes in the genome using a Mann–Whitney-Wilcoxon test. The 75 gene sample was not significantly different from the rest of the genes in the genome (p = 0.33). Thus, we conclude that the HGT candidates are not from intron-poor gene families. This evidence supports the hypothesis that the lack of introns in the HGT candidates is the consequence of their bacterial origin.

The GC content of horizontally transferred genes can be different than genes within the recipient genome [22]. We used the Mann–Whitney-Wilcoxon test to determine whether there is a difference in GC content between the 11 HGT candidates in C. graminicola and the rest of the genes encoded in the genome. We found that the GC content of the HGT candidates is not different from the rest of the genome (p = 0.337).

Functional annotation of HGT candidates

Four of the HGT candidates are involved in carbohydrate metabolism. Genes in HGT5 encode glucarate dehydratase, the enzyme that transforms D-glucarate to 5-dehydro-4-deoxy-D-glucarate + H2O in the D-glucarate degradation reaction. Genes in HGT9 encode 2-deoxy-D-gluconate 3-dehydrogenase, which is involved in pentose and glucuronate interconversion. Additionally, HGT8 genes encode acetyl-CoA synthetase, which is involved in glucose biosynthesis and in the biosynthesis of the fatty acids and in the Krebs cycle [40]. Finally, HGT4 encodes oligoxyloglucan reducing-end-specific cellobiohydrolase enzymes (glycoside hydrolases belonging to CAZy family GH74), which putatively breakdown carbohydrates in the plant cell wall [41].

Candidate HGT2 (glutathionylspermidine synthase) is the only HGT candidate that does not have vertically transmitted homologs in the Colletotrichum genomes and this may represent the acquisition of a completely new gene family by HGT. This enzyme catalyzes the synthesis of glutathionylspermidine and ADP+ orthophosphate from glutathione and spermidine. This reaction is well understood in trypanosomatid parasites and E. coli and some evidence supports the role of glutathionylspermidine synthase in detoxification of redox reactions [42, 43].

At least 4 HGT candidates have clear associations with plant interactions and virulence. Members of HGT4 encode secreted glycoside hydrolases belonging to CAZy family GH74 [39] which degrade cellulose in the plant cell wall and are important for virulence in Magnaporthe oryzae and other fungi [41, 44, 45]. Members of HGT9, encoding a short chain dehydrogenase, have significant similarity to virulence factors in the Bacterial Virulence Factors Database (Release 3) and the PHI-Base V3.4 database of proteins with roles in pathogen/host interactions. A homolog of this protein in Cochliobolous heterostrophus, OXI1, is required for biosynthesis of the secondary metabolite T-toxin [46] and null mutants of the gene show reduced virulence. Members of HGT8 have strong similarity to SidI from Aspergillus fumigatus which plays a role in siderophore biosynthesis and subsequently, virulence [47]. Finally HGT1, a family of arginosuccinate lyase genes, share homology with ARG1 of Fusarium oxysporum f. sp. melonis[48]. Null mutants in F. oxysporum have reduced virulence, linking virulence with arginine biosynthesis.

HGT candidates are expressed during plant infection

We reasoned that if the HGT candidates have roles in virulence then they should be expressed during infection of the host. Therefore, we used the transcriptome profiling experiments reported by O’Connell et al. [30] to identify differences in the level of expression of the HGT gene candidates from C. graminicola and C. higginsianum during three stages of infection of their respective hosts, maize and Arabidopsis. We also examined the expression of a selection of vertically transmitted homologs to the HGT candidates during the infection process. The transcriptional profiling experiment was conducted at three important time points during the infection process: early (in vitro or in planta appressoria, a specialized fungal structure used to penetrate into the plant; VA, PA), middle (biotrophic phase; BP) and late (necrotrophic phase; NP). Candidates GLRG_11091 (HGT2) and GLRG_11966 (HGT2) of C. graminicola and candidates CH063_02340 (HGT2), CH063_10640 (HGT2) and CH063_01625 (HGT5) of C. higginsianum are strongly upregulated at the PA and BP time points which are very early in the infection process, suggesting a role in plant penetration or establishment of infection (Additional file 2: Table S2). The candidates GLRG_11936 (HGT6), GLRG_06163 (HGT8), CH063_01794 (HGT1), CH063_13530 (HGT11) and CH063_05456 (HGT4) are upregulated at the latest stage of infection suggesting roles in nutrient uptake.

A steady tempo of HGT events in the Pezizomycotina

We hypothesized that HGT enables fungi to adopt new ecological niches by giving them access to new nutritional substrates or enabling pathogens or endophytes to jump to new hosts. If this is true then HGT might be more common after major extinction events that may open new niches for colonization by new species. To deduce the age of the HGT events, we reconstructed a species tree with all of the fungal species observed in the HGT detection pipeline. All of the HGT events occurred after the appearance of the Pezizomycotina (Figure 1) and all but three (HGT4, HGT5 and HGT6) occurred after the appearance of the Sordariomycetes. We calculated the approximate divergence times of the lineages using a Bayesian time-measured phylogeny using the fossil Paleopyrenomycites devonicus to calibrate the tree [50] (Figure 2). We estimated the minimum ages of transference from bacteria to ancestral members of the Pezizomycotina based on the lower bound of the highest posterior density (HPD) interval. The HGT events occurred over a broad range of geological periods (from the Siluric to the Tertiary). In contrast to our expectations, there is no evidence of a burst of HGT events coinciding with major geological events. In contrast, HGT appears to be a constant, albeit rare phenomenon in the Pezizomycotina, occurring with a steady tempo during their evolution.

While many of the HGT events are ancient, none of them are broadly distributed in extant species of the Pezizomycotina (Figure 1), therefore, many lineages of fungi must have lost the genes. If HGT is associated with niche adaptation then we may find that fungi with certain lifestyles have maintained horizontally transferred genes that are required for the lifestyle.

From Figure 1 we can also estimate the number of gene losses that occurred following the HGT events, based on the absence of orthologs to the HGT candidates in the whole genome sequences that are available in GenBank. We identified the species with presence/absence of each of eleven groups of candidates (Additional file 1: Table S3). 88.72% of Pezizomycotina species (with complete genome available in GenBank) had lost all members of HGT5, 91.73% of species had lost all the members of HGT4 and 93.98% of species had lost all members of HGT6. In contrast, HGT1, HGT2, HGT4, HGT5 and HGT6 are always present in species of the genus Colletotrichum. The lifestyle information for each species used in the tree reconstruction was compiled in Additional file 1: Table S4. This information was used to determine if there is a correlation between the lifestyle of the Pezizomycotina species and the presence/absence of the HGT candidates. We did not observe a direct relationship between the lifestyle and the presence of HGT candidates except in HGT4 and HGT6, which are exclusively found in pathogenic species including the entomopathogens Cordyceps militaris and Beauveria bassiana (in the case of HGT6).

Detection of candidates

The first BLAST search was performed using a database of proteins from organisms with complete proteome available in UniProt (downloaded 19-12-2011) (http://www.uniprot.org). Python scripts were used to extract the sequence, description and the two highest taxonomic levels of each sequence. The two highest taxonomic levels of UniProt correspond to kingdom and phylum in Archaea and Bacteria and to superkingdom and kingdom in Eukaryota. The two highest levels were selected to identify the taxonomic assignment at the level of kingdom of the sequences in the database. With this database a BLASTP (v2.2.29+) [58] search (maximum e-value 10^-5) was performed with all putative proteins predicted in the genomes of C. graminicola M 1001, C. higginsianum IMI 349063 [30] and C. gloeosporioides[32]. We selected proteins that having at least 80% of the top 120 hits (e-value e-5) with a taxonomic classification other than fungi as candidates for further analysis. To choose this threshold, we developed a set of true positive HGT genes by combining the HGT candidates from Richards et al. [12], Schmitt and Lumbsch [21] and Richards [28]. Using the true positive set of proteins, we performed BLAST searches and then constructed phylogenetic trees with the BLAST hits and the query sequence. We evaluated the trees for topologies consistent with HGT of the query sequence. Not all the cases showed HGT patterns after the tree evaluation. 80% was the threshold that enabled us to identify the maximum number of the true positives HGT candidates.

Next, we subjected the HGT candidates to three phylogenetic analyses using different sets of homologous sequences in each phylogeny. The first phylogeny was constructed with homologous sequences from the UniProt complete proteome database (downloaded 19-12-2011) (476 proteins from the 3 species were selected). For the second phylogeny we performed a BLAST search of the GenBank nr database (downloaded 03-02-2012) and included the 20 best hits from each kingdom (Archaea, Bacteria and Eukaryota) to avoid a possible under or overrepresentation for the abundance of sequences from any one kingdom in the BLAST results. The sixty best BLAST hits were aligned using MAFFT [59] and the alignments were edited with Gblocks [60] to remove poorly aligned regions from alignments. The alignments were used to construct phylogenetic trees using PhyML [61] using the default parameters. Each tree was evaluated to identify those with topologies consistent with HGT.

The third phylogeny was constructed with the best 100 BLAST hits from the nr database (downloaded 20-05-2012). The phylogenetic trees constructed as described in the second phylogeny and were evaluated manually, selecting only those that have well-supported topologies that are clearly concruent with HGT. Asymmetric or ladder-shape trees were excluded because such tree topologies are often a signal of long branch attraction or lack of phylogenic information [62]. Additionally, candidates with few homologues or with low sequence similarity to all of their BLAST hits were also excluded. In cases where only hits from two kingdoms were obtained (i.e. Bacteria and Eukaryota), high sequence similarity (a minimum of 30% pairwise similarity) and coverage (over 80% coverage) were required to consider them as candidates. The BLAST searches and tree evaluations were performed serially rather than in parallel to minimize the number of manual phylogenetic tree evaluations required (Additional file 1: Figure S1).

Homologs to the HGT candidates were identified in the C. graminicola, C. higginsianum and C. gloeosporioides genomes using BLAST (v2.2.29+). When no evidence of homology in one or more species was found, we additionally searched the EST database of NCBI and the RNA-seq sequences of O’Connell et al. [30] to find evidence of homologous sequences. To verify the presence of homologous sequences in C. graminicola, BLAST searches were performed in the genome of 5 sequenced strains of this species by Rech [63]. Additionally, four newly available Colletotrichum genomes were used to verify the absence/presence of each candidate in other members of the genus. This verification was performed through BLAST searches. The genomes used for this purpose were: Colletotrichum gloeosporioides 23 (teleomorph, Glomerella cingulata) (http://www.jgi.doe.gov), Colletotrichum fiorinae MH 18 (teleomorph, Glomerella acutatum) (http://www.jgi.doe.gov), Colletotrichum gloeosporioides Nara gc5 [31] and Colletotrichum orbiculare MAFF 240422 [31].

The 11 groups of candidates were subjected to three tests to ensure that there is sufficient phylogenetic signal to support the HGT hypothesis. For each candidate, we manually inspected the results of the BLAST search to the nr database (downloaded 12-11-2012), selecting appropriate percent identity and coverage thresholds to eliminate distantly related and poorly aligning homologs. Each HGT candidate and its homologs were aligned and analyzed with TREE-PUZZLE [64]. The TREE-PUZZLE algorithm uses a maximum-likelihood approach to resolve the phylogenetic tree of all possible combinations of four sequences in an alignment (quartets). If the maximum likelihood values of the three possible tree topologies of a quartet are too similar, the quartet is labeled as unresolved. Thus, the percentage of unresolved quartets for each sequence can be used as a measure of the phylogenetic signal for that sequence. Using TREE-PUZZLE, we calculated the percentage of unresolved quartets for each sequence and, following the TREE-PUZZLE documentation, discarded those with more than 10% unresolved quartets. Next, we manually edited each multiple sequence alignment, removing regions of low sequence similarity where alignment errors are more likely. A tree was reconstructed for the edited and unedited alignments and the trees were compared to ensure that there was no change in topology due to alignment editing. In no case did we observe a change in topology as a result of alignment editing.

We evaluated the stability of the tree topology when different sets of proteins are used. For each candidate we reconstructed several trees by randomly removing several proteins and comparing the resulting tree topology with that of the original. HGT2, HGT6 and HGT3 had different topologies in this analysis but were stable in the rest of the tests to evaluate the phylogenies.

The data sets were used to perform topology tests to determine whether a tree toplogy supporting vertical inheritance or a star phylogeny are equally well supported by the sequence alignments. Some of the candidates did not have high similarity with fungal sequences or sequences of any other kingdom more than bacteria and in these cases the topology tests could not be performed because the putative HGT sequences could not be constrained with the fungal branch in the tree to evaluate vertical inheritance as an alternative hypothesis. To generate the tree consistent with vertical inheritance, MrBayes [65] was used to constrain all fungal proteins in one branch and the bacterial proteins in other branch. The last topology included in the tests was the star phylogeny. TREE-PUZZLE was used to perform the Expected Likelihood Weight (ELW) test and Shimodaira and Hasegawa (SH) tests. For each candidate with sequence homology in bacteria and fungi the HGT trees were chosen as a better explanation for the data than the tree supporting vertical inheritance or the star phylogeny. Finally, for each alignment, a phylogenetic tree was constructed with PhyML [61] and with MrBayes [65] and the similarity of the two trees was evaluated with T-Coffee [66]. These analyses are summarized in Table 1 and the trees are shown in Figure 3 and Additional file 1: Figures S2-S11. In each case, the HGT candidates form a monophyletic group within the bacterial lineage, indicating a single horizontal transfer event.

To determine whether any of the HGT candidates are mitochondrial genes that were incorrectly assembled into the nuclear genome sequence, or are mitochondrion to nuclear HGT events, we performed a BLASTP (v2.2.29+) search vs all mitochondrial proteins available in the RefSeq database [67]. Proteins of HGT9 show high similarity (e-value 6.17e-22) to protein XP_001875307 a 3-hydroxyacyl-CoA-dehydrogenase (EC 1.1.1.35) from the mitochondrion of Laccaria bicolor. However, the most similar BLAST hits are from bacteria (e.g. e-value of 9.33e-116; YP_005460258 from Actinoplanes missouriensis). Therefore, we conclude that this is neither a case of mitochondrion to nucleus transfer nor an error in the genome assembly.

To detect signs of contamination of the genome sequences and identify events of HGT of clusters of genes, a BLASTp (v2.2.29+) search of the nr database of NCBI (online version 06-03-2013) was performed with the upstream and downstream proteins of each HGT candidate in the genome of C. graminicola, C. gloeosporioides and C. higginsianum. Most of the neighboring genes showed the expected distribution of BLAST hits with Colletotrichum and fungal proteins (Additional file 1: Table S1). Only proteins CH063_08061T0 and CGSP_9355, the flanking genes of members of HGT3 (CH063_08062T0 and CGSP_9355 respectively) have BLAST hits that include bacterial proteins, suggesting the possibility of contamination of the genome sequence with bacterial DNA. However, proteins CH063_08061T0 and CGSP_9355 are orthologous and may have originated from the same HGT event as the members of HGT3 but the phylogenetic evidence was insufficient to conclude this. Additionally, all the candidate proteins were analyzed in the Web server CENSOR [68] to detect repetitive or transposable elements but no evidence of repetitive sequences was found. The entire pipeline is summarized in the Additional file 1: Figure S1.

A maximum posterior tree was constructed with MrBayes, performing 2,000,000 generations of samples, using the substitution matrix and model predicted by MODELGENERATOR but allowing the program to calculate the proportion of invariable sites and the alpha parameter for gamma distribution. Two Multiple Chain Markov Chain Monte Carlo (MC³) searches were conducted with four chains each (three heated and one cold). The convergence between them was checked using a sample frequency of 1000 generations. A burn-in of 25% of generations was excluded to reconstruct the Bayesian consensus tree. The topology differences between maximum likelihood and Bayesian trees were quantified with the ratio of identical nodes performed in T-coffee [66]. The percentage of bootstrap and the posterior probability index were joined in the maximum likelihood tree using the program TreeGraph2 [69].

To determine if HGT candidates were the result of contamination a BLASTP (v2.2.29+) (vs the nr database of NCBI) search of the upstream and downstream protein of each candidate in the genome of C. graminicola, C. gloeosporioides and C. higginsianum was performed. This search was done if the neighbor gene was located in the same contig of the candidate, otherwise the analysis was not performed. Additionally, to detect putative mitochondrial to nucleus gene transfers, a BLASTP (v2.2.29+) search against all mitochondrion proteins of RefSeq database [67] was performed. Also all the candidates were submitted to CENSOR [68] to detect the presence of repetitive or transposable elements.

Functional annotation of HGT candidates

The program Blast2GO 2.8 [70] was used to annotate the HGT candidates with Gene Ontology (GO) terms. The annotations were verified with InterProScan 4 [71] and the putative biochemical functions of candidates were predicted with BRENDA (Release 2012.02) [35], KEGG (update 13-12-2012) [36] and MetaCyc 18.5 [37]. Additionally, MEROPS 9.11 [38] and CAZy [39] databases were explored to further annotate the functions of proteases and carbohydrate active enzymes. To predict the cellular localization of candidates WOLF PSORT v0.2 [72] and SignalP 4.1 server [73] were used.

BLASTP (v2.2.29+) searches of the HGT candidates were performed against PHI-base V3.4 [74], Virulence Factors Database (VFDV) (Release 3) [75] and Database of Fungal Virulence Factors (DFVF) [76] to detect proteins implied in pathogenicity.

The expression data were extracted from O’Connell et al. [30] for C. graminicola and C. higginsianum. The data is presented in the Additional file 2: Table S2.

GC content and intron content of the candidates

To determine the presence/absence of introns in HGT candidates a manual inspection of genes was made. To compare the difference of introns content among the candidates and C. graminicola a Mann–Whitney-Wilcoxon test was performed using R v2.13.1 [77]. The intron content of each C. graminicola gene was calculated from the information available on the Broad institute website (http://www.broadinstitute.org).

The program CodonW [78] was used to calculate the GC content of the whole Colletotrichum genome and HGT candidates. These indices were used to calculate the differences of GC content among candidates and whole genomes.

Species tree reconstruction

To estimate the phylogenetic relationship of the species involved in HGT events, all fungal species found in BLAST searches of the HGT candidates were used (54 species in total). The complete proteome of these species was obtained from the UniProt (http://www.uniprot.org), GenBank (http://www.ncbi.nlm.nih.gov/genbank/), Joint genome institute (http://genome.jgi-psf.org/programs/fungi/index.jsf?projectList) and Broad Institute (http://www.broadinstitute.org/) databases (all databases downloaded 03-05-2013). To reconstruct the species tree, the amino acids inferred from six nuclear genes were chosen. The proteins selected (FG533, FG570, FG832, MS277, MS413 and MS456) from FunyBase [79] demonstrate to be good phylogenetic markers for fungi species trees reconstruction [80] and for that reason these were selected as a query to perform the BLAST searches in the proteomes of the 54 fungal species. Homology was verified making a multiple sequence alignment and tree reconstruction with MAFFT and PhyML respectively to make clear the orthology relationship between taxa. The topological congruence among the protein trees was checked before accepting the protein in the analysis. The six proteins were concatenated for each taxa. The concatenated proteins were aligned with three different programs, MAFFT, MUSCLE [81] and CLUSTALW [82] to evaluate the differences in the phylogenetic reconstruction when different software is used. When the alignment was chosen TrimAl [83], GBLOCKS [60] and Guidance [84] were used to edit it. A tree was reconstructed with each edited and unedited alignment with PhyML with 100 bootstrap repetitions. The best tree was selected by the alignment that produced the tree with the highest bootstrap values. PartitionFinder [85] was used to detect accurate models for the final alignment, using each protein as a partition. The models predicted by PartitionFinder were used in RaxML to calculate the maximum likelihood tree starting with 100 random trees. Finally, a non-parametric bootstrap analysis with 100 replications was performed and the results were summarized in the maximum likelihood tree.

Molecular clock analysis

To estimate the putative age of the transferred genes a fossil calibrated molecular clock analysis was performed. To avoid problems of convergence in the calculations of the calibrated tree, we selected 15 species from the 54 taxa used in the species tree. The same matrix of 6 concatenated proteins of the species tree was used. The analysis was performed with the BEAST v1.7.5 software package [86]. For each partition (each protein) the LG + I + G model was used. To allow uncorrelated rates of evolution across the tree we use a lognormal relaxed clock model, implementing a Yule process as a tree prior. We use a maximum likelihood tree estimated in PhyML under the LG + I + G model as starting tree. To calibrate the tree we use the estimation of the Paleopyrenomycites devonicus fossil age of 400 million years ago (mya) as the lower bound for the Pezizomycotina crown [87]. A lognormal distribution with a mean of 460 (estimated from the results of Lucking et al. [87]), standard deviation of 1 and offset of 400 was used as a prior for the time to the most recent common ancestor (TMRCA) of the Pezizomycotina. Two independent BEAST runs of 15 million generations each were performed. Data was sampled every 1500 generations. The convergence of two runs was visualized with TRACER v1.5 [88] and the Log files and tree files were combined with LogCombiner v1.7.5 [86] dismissing a percentage of the sample in agree with TRACER plots of each run (23.3% for each one). With the remaining trees a maximum clade probability tree was calculated using TreeAnnotator v1.7.5 [86]. The resultant tree was visualized with the FigTree software (http://www.tree.bio.ed.ac.uk/software/figtree/).