The Genome of the Fungal Pathogen Verticillium dahliae Reveals Extensive Bacterial to Fungal Gene Transfer.

Horizontal gene transfer (HGT) involves the transmission of genetic material between distinct evolutionary lineages and can be an important source of biological innovation. Reports of interkingdom HGT to eukaryotic microbial pathogens have accumulated over recent years. Verticillium dahliae is a notorious plant pathogen that causes vascular wilt disease on hundreds of plant species, resulting in high economic losses every year. Previously, the effector gene Ave1 and a glucosyltransferase-encoding gene were identified as virulence factor-encoding genes that were proposed to be horizontally acquired from a plant and a bacterial donor, respectively. However, to what extent HGT contributed to the overall genome composition of V. dahliae remained elusive. Here, we systematically searched for evidence of interkingdom HGT events in the genome of V. dahliae and provide evidence for extensive horizontal gene acquisition from bacterial origin.

Introduction

Genetic information is generally vertically transferred from parents to their offspring. However, genetic information can also be transmitted laterally between reproductively isolated species, often referred to as horizontal gene transfer (HGT). It has been well established that HGT plays a significant role in the adaptive evolution of prokaryotic species (Eisen 2000; Koonin et al. 2001; Bapteste et al. 2009). Three well-characterized mechanisms contribute to DNA uptake by prokaryotes, namely transformation, conjugation, and transduction. With transformation, a DNA fragment from a dead, degraded bacterium or another donor enters a competent recipient bacterium (Johnston et al. 2014), whereas conjugation is the active DNA transfer between prokaryotic cells by direct cell-to-cell contact or by a bridgelike connection between two cells (Norman et al. 2009). Finally, transduction involves the transfer of a DNA fragment from one prokaryotic cell to another by a virus or viral vector. HGT in prokaryotes takes place at all taxonomic levels: from individuals of the same population up to interkingdom transfers.

HGT also contributes to the evolution of eukaryotes, although it occurs much less frequent than in prokaryotes (Kurland et al. 2003; Bock 2010). HGT has played an important role in the evolution of pathogenic traits in plant pathogens (Soanes and Richards 2014). For instance, the wheat pathogenic fungi Pyrenophora tritici-repentis and Phaeosphaeria nodorum both contain the near-identical effector gene ToxA that encodes a host-specific toxin that acts as pathogenicity factor (Friesen et al. 2006). Compelling evidence revealed that ToxA was acquired by Pyrenophoratritici-repentis from Phaeosphaerianodorum via HGT, giving rise to pathogenicity of the former species on wheat plants (Friesen et al. 2006). Interestingly, ToxA was similarly found in the genome of yet another wheat pathogen, Bipolaris sorokiniana (McDonald et al. 2018).

Although knowledge on mechanisms of HGT in eukaryotes remains limited, HGT in filamentous plant pathogens is thought to occur more often between closely related species because donor and recipient species share similar genomic architectures (Mehrabi et al. 2011). For the filamentous plant pathogenic fungus Fusarium oxysporum, an entire chromosome can be transferred from one strain to another through coincubation under laboratory conditions in vitro (Ma et al. 2010; van Dam et al. 2017).

Intriguingly, interkingdom HGT also contributed to the genome evolution of filamentous plant pathogens despite the evolutionary distant relation with the donor. For example, phylogenetic analyses of plant pathogenic oomycete species revealed 34 gene families that have undergone HGT between fungi and oomycetes (Richards et al. 2011). The repertoire of HGT candidates includes genes encoding proteins with the capacity to break down plant cell walls and effector proteins to support interaction with host plants (Richards et al. 2011). Furthermore, genomic and phylogenetic analyses of the fungus Fusarium pseudograminearum, the major cause of crown and root rot of barley and wheat in Australia, have revealed that a novel virulence gene was horizontally acquired from a bacterial species (Gardiner et al. 2012).

The fungal genus Verticillium contains nine haploid species plus the allodiploid V. longisporum (Inderbitzin et al. 2011; Depotter et al. 2017). These ten species are phylogenetically subdivided into two clades: Flavexudans and Flavnonexudans (Inderbitzin et al. 2011; Shi-Kunne et al. 2018). The Flavnonexudans clade comprises Verticilliumnubilum, Verticilliumalfalfae, Verticilliumnonalfalfae, Verticilliumdahliae, and V. longisporum, whereas the Flavexudans clade comprises Verticilliumalbo-atrum, Verticilliumisaacii, Verticilliumtricorpus, Verticilliumklebahnii, and Verticilliumzaregamsianum (Inderbitzin et al. 2011). Among these Verticillium spp., V. dahliae is the most notorious plant pathogen that is able to cause disease in hundreds of plant species (Fradin and Thomma 2006; Inderbitzin and Subbarao 2014). Furthermore, Verticilliumalbo-atrum, V. alfalfae, V. nonalfalfae, and V. longisporum are pathogenic, albeit with narrower host ranges (Inderbitzin and Subbarao 2014). Although the remaining species Verticilliumtricorpus, Verticilliumzaregamsianum, Verticilliumnubilum, Verticilliumisaacii, and Verticilliumklebahnii have incidentally been reported as plant pathogens too, they are mostly considered saprophytes that thrive on dead organic material and their incidental infections should be seen as opportunistic (Ebihara et al. 2003; Inderbitzin et al. 2011; Gurung et al. 2015). Verticillium spp. are considered to be strictly asexual, yet various mechanisms contributing to the genomic diversity of V. dahliae have been reported (Seidl and Thomma 2014, 2017; Faino et al. 2016), including HGT (Klosterman et al. 2011; de Jonge et al. 2012). Previously, Ave1 and a glucosyltransferase-encoding gene were found to be acquired from a plant and a bacterial donor, respectively, both of which were found to contribute to V. dahliae virulence during plant infection (Klosterman et al. 2011; de Jonge et al. 2012). These two interkingdom HGT events inspired us to study the extent and potential impact of interkingdom HGT to V.dahliae.

Materials and Methods

Construction of Local Protein Database and Protein BLAST

The protein database used for BLAST analyses was generated using the reference proteomes from UniProtKB excluding viruses (downloaded July 7, 2016). Entries from UniProt were renamed to include their UniProt ID and NCBI taxonomy ID. The proteome database contains proteomes of 5,074 species, including 4,305 prokaryotes and 769 eukaryotes. The prokaryotic species of this database comprise 4,123 bacteria and 182 archaea. The bacterial species comprise 1,737 Terrabacteria, 1,383 Proteobacteria, 368 Sphingobacteria, and 57 Plancterbacteria. The eukaryotic species of the database comprise 373 fungi and 14 oomycete species. Among the fungi species are 257 ascomycete and 84 basidiomycete species, respectively. Three Verticillium species are included in this database, namely V. dahliae (strain VdLs17), V. alfalfae (strain VaMs102), and V. longisporum (strain VL1) (Klosterman et al. 2011). Sequences shorter than ten amino acids were removed from the database. All predicted proteins sequences from V. dahliae strain JR2 were queried to search homologs in the local protein database (UniProtKB) using BlastP (E-value = 0.001, max target seqs = 1,000). We searched for the presence of the HGT candidates in one strain per Verticillium species that were studied previously (Depotter et al. 2017; Shi-Kunne et al. 2018). Protein sequences of each HGT candidate were queried against the proteome of each strain using BlastP (with default settings).

Identification of HGT Candidates

Alien Index (AI) scores were calculated using a custom-made Python script which applied the following formula: , where bbhG and bbhO are the E-values (<10−3) of the best BLAST hit from the ingroup and outgroup, respectively. Best BLAST hits were determined by the highest bit score. Genes with AI score >1 were classified as positive AI genes. For each AI positive gene, homologs were selected based on the BLAST criteria that at least 70% of query length is covered by at least 60% of hit sequence using a custom-made Python script. These homologs were then aligned using MAFFT with default setting (Katoh and Standley 2013). Alignments were subsequently curated using Gblocks (Castresana 2000) with nonstringent parameters. Phylogenetic trees of the curated alignments were constructed using FastTree (Price et al. 2010) with gamma likelihood and WAG amino acid substitution matrix. Phylogenetic trees were automatically evaluated using a custom-made Python script that removes phylogenetic trees with candidate genes that are directly adjacent to non-Verticillium fungal species (ingroup species) rather than nonfungal species (outgroup species). Subsequently, new phylogenetic trees (for manual inspection) were reconstructed using RAxML (Stamatakis 2014) with automatic substitution model determination for branching values.

Gene Expression Analysis

To obtain RNA-Seq data for V. dahliae grown in culture medium, isolate JR2 was grown for 3 days in potato dextrose broth with three biological replicates. To obtain RNA-Seq data from V. dahliae grown in planta, the roots of 3-week-old uprooted Arabidopsisthaliana (Col-0) plants were dipped in a suspension of 106 conidiospores per ml of water for 5 min and replanted in soil. After root inoculation, plants were grown in individual pots in a greenhouse under a cycle of 16 h of light and 8 h of darkness, with temperatures maintained between 20 and 22 °C during the day and a minimum of 15 °C overnight. At 21 days postinoculation, three samples per treatment containing ten stem fragments of 3 cm length, harvested above the soil line, and all the flowering stems from ten plants were pooled. Plant material was freeze dried for 24 h and total RNA was extracted based on TRIzol RNA extraction (Simms et al. 1993). cDNA synthesis, library preparation (TruSeq RNA-Seq short insert library), and Illumina sequencing (single-end 50 bp) was performed at the Beijing Genome Institute (Hong Kong, China). In total, ∼2 and ∼1.5 Gb of filtered reads were obtained for the V. dahliae samples grown in culture medium and in planta, respectively. RNA-Seq data were submitted to the Sequence Read Archive database under the accession number: SRP149060 (accession: PRJNA473305).

The RNA sequencing reads were mapped to their previously assembled genomes using the Rsubread package in R (Liao et al. 2013). The comparative transcriptomic analysis was performed with the package edgeR in R (v3.4.3) (Robinson et al. 2010; McCarthy et al. 2012). Genes are considered differently expressed when P value <0.05 with a log 2-fold-change ≥1. P values were corrected for multiple comparisons according to Benjamini and Hochberg (1995).

Results

Interkingdom HGT to V. dahliae and Other Ascomycete Fungi

To systematically search for genes in the V. dahliae genome that are derived from interkingdom HGT, we downloaded a recent UniProtKB proteome database that contains proteomes of 5,074 species, including 4,123 bacteria, 182 archaea, and 769 eukaryotes (fig. 1). The database contains proteomes of 373 fungal species, including V. dahliae (strain VdLs17), V. alfalfae (strain VaMs102), and V. longisporum (strain VL1) (Klosterman et al. 2011). In our analyses, we focused on the complete telomere-to-telomere genome assembly of V. dahliae strain JR2 (Faino et al. 2015). We first queried each of the 11,430 protein sequences of V. dahliae strain JR2 using BLAST against the aforementioned UniProtKB protein database. Subsequently, we utilized the AI method (Gladyshev et al. 2008; Alexander et al. 2016) that generates AI scores for each V. dahliae gene based on the comparison of best BLAST hits with ingroup and outgroup species; in this case non-Verticillium fungal species and nonfungal species, respectively (fig. 2). In this manner, we queried for HGT candidates that were acquired only by Verticillium spp., or that were acquired by ancestral fungal species but retained only by Verticillium spp. during evolution. When the best BLAST hit of a gene within the outgroup is better than the best hit within the ingroup, the AI score is positive (i.e., AI > 1) and this gene is considered a potential HGT candidate. This AI analysis yielded 42 HGT candidates with AI positive scores (AI > 1) for V. dahliae JR2 (fig. 2), which were further analyzed. To this end, we aligned protein sequences of all the homologs for each candidate (from ingroups and outgroups) and constructed phylogenetic trees that were automatically evaluated to remove phylogenetic trees with candidate genes that are directly adjacent to non-Verticillium fungal species (ingroup species) rather than nonfungal species (outgroup species). However, candidate genes that are closely related to the homologs of outgroup species in a phylogenetic tree can also result from duplication followed by multiple gene losses in the ingroup species. Typically, in such a phylogenetic tree, paralogs of the candidate gene cluster with genes of non-Verticillium fungal species (ingroup species) and form separate branches from the branch of the candidate gene. Thus, we removed the candidates with phylogenetic trees that can be explained by gene losses in multiple species. In the end, we obtained two HGT candidates (HGT-1 and HGT-2) that were acquired only by Verticillium spp., or that were acquired by ancestral fungal species but kept only by Verticillium spp. during evolution (table 1). To check whether these two V. dahliae candidates are also present in other Verticillium spp., we queried them using BLAST against the proteomes of previously published Verticillium spp. (Depotter et al. 2017; Shi-Kunne et al. 2018). We found that candidate VDAG_JR2_Chr4g10560 (HGT-1) is present in all Verticillium species and that VDAG_JR2_Chr8g03340 (HGT-2) is only present in V. alfalfae, V. nonalfalfae, and V. longisporum (figs. 3 and 4).

. 1.

—Phylogenetic tree of the species in the selected local database. Different colors represent different taxonomical groups in the UniProtKB protein database. The number of species (spp.) on each branch is indicated. The phylogenetic position of Verticillium spp. is indicated on the right.

. 2.

—AI-based HGT detection method for detecting HGT candidates. (A) Simplified phylogenies of ingroup and outgroup. In AI-A setting, ingroup and outgroup species are non-Verticillium fungal species and nonfungal species, respectively. In AI-B setting, ingroup and outgroup species are nonascomycete fungal species and nonfungal species, respectively. (B) Distribution of AI sores. AI scores from two settings were calculated for every gene in the genome of Verticillium dahliae strain JR2 and they were ordered by decreasing AI score. (C) Numbers of genes with AI positive scores from two different settings.

Table 1

Information on HGT Candidates

HGT ID	Gene ID	Putative Donor	Putative Gene Product	Functional Categorya	Secreted	In Planta Expression
HGT-1	Chr4g10560	Bacterial	Unknown	Unknown	No	No
HGT-2	Chr8g03340	Bacterial	Pyridoxal-dependent decarboxylase	Pyridoxal phosphate binding	No	No
HGT-3	Chr1g00380	Bacterial	Glycosyl hydrolase	Carbohydrate metabolism	Yes	No
HGT-4	Chr1g00710	Bacterial	Unknown	Unknown	Yes	No
HGT-5	Chr1g17450	Bacterial	Glycosyl hydrolase (GH27)	Carbohydrate metabolism	Yes	No
HGT-6	Chr1g24620	Bacterial	Luciferase-like monooxygenase	Oxidoreductase activity	No	No
HGT-7	Chr1g24680	Bacterial	Glutathione-dependent formaldehyde-activating enzyme	Detoxification of formaldehyde	No	Yes
HGT-8	Chr1g25610	Bacterial	Glutathione S-transferase	Detoxification of reactive Electrophillic compounds	No	No
HGT-9	Chr1g28650	Bacterial	Glycosyl hydrolases	Carbohydrate metabolism	No	No
HGT-10	Chr2g00370	Bacterial	Unknown	Unknown	No	No
HGT-11	Chr2g02230	Bacterial	Carbohydrate esterase (CE1)	Carbohydrate metabolism	No	No
HGT-12	Chr2g02370	Bacterial	3-Demethylubiquinone-9 3-methyltransferase	Unknown	No	Yes
HGT-13	Chr2g02510	Bacterial	Methyltransferase	Unknown	No	No
HGT-14b	Chr2g04700	Bacterial	Glycosyl transferase family group 2	Carbohydrate metabolism	No	No
HGT-15	Chr2g07220	Bacterial	Pyridine nucleotide-disulfide oxidoreductase	Oxidoreductase activity	No	No
HGT-16	Chr2g08740	Bacterial	Enolase	Carbohydrate metabolism	No	Yes
HGT-17	Chr3g05430	Bacterial	Chondroitin AC/alginate lyase	Glycosaminoglycans degradation	Yes	Yes
HGT-18	Chr3g06830	Bacterial	Acetyltransferase (GNAT) family	Transcriptional regulation	No	No
HGT-19	Chr3g12010	Bacterial	Glycosyl hydrolase (GH74)	Hydrolase activity	Yes	Yes
HGT-20	Chr3g13080	Bacterial	Phosphoinositide phospholipase	Hydrolase activity	Yes	Yes
HGT-21	Chr4g04370	Bacterial	Membrane protein	Unknown	No	No
HGT-22	Chr4g11370	Bacterial	Glycosyl hydrolase (GH74)	Carbohydrate metabolism	No	No
HGT-23	Chr4g11740	Bacterial	Amidohydrolase	Hydrolase activity	Yes	No
HGT-24	Chr4g11940	Bacterial	Carbohydrate esterases (CE1)	Carbohydrate metabolism	Yes	No
HGT-25	Chr5g01710	Bacterial	l-Asparaginase	Hydrolase activity	No	No
HGT-26c	Chr5g02170	Plant	Plant natriuretic peptide	Virulence factor	Yes	Yes
HGT-27	Chr5g05890	Bacterial	Glycosyl hydrolases (GH43/29)	Carbohydrate metabolism	Yes	No
HGT-28	Chr5g09290	Bacterial	Protein of unknown function	Unknown	No	Yes
HGT-29	Chr6g02080	Bacterial	Bacteriocin-protection, YdeI or OmpD-associated protein	Oxidase activity	No	No
HGT-30	Chr6g03000	Bacterial	Pyridoxal phosphate–dependent enzyme	Metabolic activity	No	No
HGT-31	Chr6g05320	Bacterial	OsmC-like protein	Oxidative stress regulation	No	No
HGT-32	Chr6g05650	Bacterial	Luciferase-like monooxygenase	Oxidoreductase activity	No	Yes
HGT-33	Chr6g10680	Bacterial	Acetyltransferase (GNAT) family	N-Acetyltransferase activity	No	Yes
HGT-34	Chr7g01340	Bacterial	M6 family metalloprotease	Peptidase activity	No	No
HGT-35	Chr7g01630	Bacterial	Glycosyl hydrolases (GH136)	Carbohydrate metabolism	No	No
HGT-36	Chr7g03210	Bacterial	Aminotransferases	Catalytic activity	No	No
HGT-37	Chr7g03340	Bacterial	Alpha/beta hydrolase	Hydrolytic activity	No	Yes
HGT-38	Chr7g03590	Bacterial	Insecticide toxin TcdB	Unknown	No	No
HGT-39	Chr8g00620	Bacterial	Glycerate kinase	Glycerate kinase activity	No	No
HGT-40	Chr8g06130	Bacterial	Unknown	Unknown	No	No
HGT-41	Chr8g08880	Bacterial	Nucleotidyltransferase	DNA repair	No	No
HGT-42	Chr8g11000	Bacterial	Glycosyl hydrolase (GH67)	Carbohydrate metabolism	Yes	Yes
HGT-43	Chr8g11270	Bacterial	Glycosyl hydrolases (GH43/26), (CBM42)	Carbohydrate metabolism	Yes	No
HGT-44	Chr8g01080	Bacterial	Peptidase m20	Peptidase activity	No	No

Functional annotation was performed by searching for conserved domains in each of these proteins using Interproscan (Jones et al. 2014).

Previously identified candidate (Klosterman et al. 2011).

Previously identified candidate (de Jonge et al. 2012).

. 3.

—Evolutionary relationship of HGT-1 homologs. Protein sequences of HGT-1 homologs were aligned (using MAFFT) and the resulting alignment was used to infer a maximum-likelihood phylogeny (using RAxML). The phylogeny suggests that HGT-1 is transferred from a bacterial species. Different colors depict different groups or species. A more detailed part of the tree that contains Verticillium species is shown on the right. Orange squares indicate branches with bootstrap values ≥60.

—Evolutionary relationship of HGT-2 homologs. Protein sequences of HGT-2 homologs were aligned (using MAFFT) and the resulting alignment was used to infer a maximum-likelihood phylogeny (using RAxML). Different colors depict different groups or species. The phylogeny suggests that HGT-2 is transferred from a bacterial species. A more detailed part of the tree that contains Verticillium species is shown on the right. Orange squares indicate branches with bootstrap values ≥60.

In order to search for HGT candidates that are not only present in Verticillium spp. but also may be present in other ascomycete species, we again used the AI method, albeit with adjusted AI group settings. In this case, we set nonfungal species and nonascomycete fungal species as outgroups and ingroups, respectively. AI scanning resulted in 603 genes with AI positive scores (AI > 1), which were selected for further assessment through phylogenetic tree analysis. After assessing phylogenetic trees, we identified 43 additional HGT candidates with nonfungal origins that are also present in other fungal ascomycete species (table 1). Among these candidates, two appeared as paralogs in the phylogenetic tree. The clusters of these two genes are close to each other in the tree, which indicates that the two paralogs are likely the result of a duplication and not of independent acquisitions. Therefore, we only considered 42 HGT candidates for further analyses. Of these, one is of plant origin and 41 are of bacterial origin (table 1). The candidate of plant origin is the previously identified HGT gene, Ave1 (fig. 5) (de Jonge et al. 2012). The previously reported glucosyltransferase gene (Klosterman et al. 2011) likewise was identified as HGT event in our study (fig. 6). Among these 42 candidates, 32 are present in all of the Verticillium species. Eight of the remaining ten are only found in the Verticillium species of the Flavnonexudans clade (fig. 7). Subsequently, we queried the protein sequence of each of the HGT candidates against the NCBI database using TBlastN (web version). We found that the top best hits for each of the candidates, except for Ave1, are from bacterial species, further confirming that the HGT candidates are likely acquired from bacteria (table 2).

. 5.

—Evolutionary relationship of Ave1 homologs. Protein sequences of Ave1 homologs were aligned (using MAFFT) and the resulting alignment was used to infer a maximum-likelihood phylogeny (using RAxML). Different colors depict different groups or species. A more detailed part of the tree that contains Verticillium species is shown on the right. Orange squares indicate branches with bootstrap values ≥60 or above.

. 6.

—Evolutionary relationship of Verticillium dahliae glucosyltransferase homologs. Homologs of glucosyltransferase were aligned (using MAFFT) and the resulting alignment was used to infer a maximum-likelihood phylogeny (using RAxML). Different colors depict different groups or species. A more detailed part of the tree that contains Verticillium species is shown on the right. Orange squares indicate branches with bootstrap values ≥60.

. 7.

—Presence and absence of each of the HGT candidates in each of the Verticillium species. Presence of the HGT candidates was assessed in one strain for each of the species that was studied previously (Depotter et al., 2017; Shi-Kunne et al. 2018) using BlastP (with default settings). Black dots indicate the presence of an HGT candidate in the corresponding Verticillium species.

Table 2

Overview of the Most Significant (Nonascomycete) BLAST Hits

HGT ID	Gene ID	Best (Nonascomycete) BLAST Hita
		Species	E-Value	Accession No.
HGT-1	Chr4g10560	Pseudomonas fluorescens	8.00E-31	CP012830.1
HGT-2	Chr8g03340	Solitalea canadensis	9.00E-27	CP003349.1
HGT-3	Chr1g00380	Cellulosimicrobium cellulans	3.00E-105	CP020857.1
HGT-4	Chr1g00710	Pseudarthrobacter equi	6.00E-101	XM_007833612.1
HGT-5	Chr1g17450	Streptosporangium sp.	0	CP034463.1
HGT-6	Chr1g24620	Bacillus sp.	2.00E-142	CP005586.1
HGT-7	Chr1g24680	Azospirillum sp.	3.00E-71	KX488120.1
HGT-8	Chr1g25610	Pectobacterium atrosepticum	1.00E-54	CP024956.1
HGT-9	Chr1g28650	Spirosoma aerolatum	7.00E-66	CP020104.1
HGT-10	Chr2g00370	Amycolatopsis albispora	1.00E-108	CP015163.1
HGT-11	Chr2g02230	Gemmata sp.	3.00E-41	CP011271.1
HGT-12	Chr2g02370	Janthinobacterium agaricidamnosum	5.00E-59	HG322949.1
HGT-13	Chr2g02510	Paenibacillus ihbetae	4.00E-28	AP018227.1
HGT-14	Chr2g04700	Neorhizobium sp.	2.00E-128	CP030830.1
HGT-15	Chr2g07220	Asticcacaulis excentricus	1.00E-104	CP002395.1
HGT-16	Chr2g08740	Actinoalloteichus sp.	0	CP016077.1
HGT-17	Chr3g05430	Streptomyces bingchenggensis	2.00E-135	CP002047.1
HGT-18	Chr3g06830	Thielavia terrestris	3.00E-24	XM_003652118.1
HGT-19	Chr3g12010	Actinoplanes friuliensis	6.00E-99	CP006272.1
HGT-20	Chr3g13080	Kribbella flavida	4.00E-88	CP001736.1
HGT-21	Chr4g04370	Streptomyces hygroscopicus	2.00E-38	CP003720.1
HGT-22	Chr4g11370	Paraburkholderia phenoliruptrix	3.00E-104	CP003863.1
HGT-23	Chr4g11740	Sorangium cellulosum	5.00E-84	AM746676.1
HGT-24	Chr4g11940	Pseudomonas salegens	1.00E-106	LT629787.1
HGT-25	Chr5g01710	Rhizobium phaseoli	2.00E-111	CP013525.1
HGT-26	Chr5g02170	Nicotiana attenuata	2.00E-35	XM_019410592.1
HGT-27	Chr5g05890	Sphingopyxis sp.	8.00E-151	LT598653.1
HGT-28	Chr5g09290	Comamonas testosteroni	2.00E-52	CP006704.1
HGT-29	Chr6g02080	Shinella sp.	1.00E-93	CP015740.1
HGT-30	Chr6g03000	Desulfitobacterium metallireducens	1.00E-28	CP007032.1
HGT-31	Chr6g05320	Marmoricola scoriae	4.00E-38	LT629757.1
HGT-32	Chr6g05650	Rhizobium tropici	0	CP004018.1
HGT-33	Chr6g10680	Actinoplanes teichomyceticus	2.00E-21	CP023865.1
HGT-34	Chr7g01340	Streptomyces sp.	5.00E-52	CP029541.1
HGT-35	Chr7g01630	Streptomyces sp.	5.00E-139	KF386871.1
HGT-36	Chr7g03210	Halioglobus japonicus	8.00E-66	CP019450.1
HGT-37	Chr7g03340	Gluconobacter oxydans	2.00E-24	CP004373.1
HGT-38	Chr7g03590	Psychrobacter sp.	0	CP012533.1
HGT-39	Chr8g00620	Bacillus cereus	4.00E-117	CP001176.1
HGT-40	Chr8g06130	Mesorhizobium sp.	2.00E-38	LN649230.1
HGT-41	Chr8g08880	Streptomyces albidoflavus	1.00E-72	CP004370.1
HGT-42	Chr8g11000	Candidatus koribacter	0	CP000360.1
HGT-43	Chr8g11270	Streptomyces qaidamensis	1.00E-146	CP015098.1
HGT-44	Chr8g01080	Polaromonas sp.	0	CP031013.1

HGT candidates were queried against NCBI database (web version) using TBlastN (E-value = E-06).

A High Proportion of V. dahliae HGT Candidates Encodes Secreted Proteins

Filamentous plant pathogens secrete large numbers of proteins, including plant cell wall-degrading enzymes and effector proteins, to interact with host plants (Cook et al. 2015). We used a combination of SignalP, TMHMM, and TargetP (Shi-Kunne et al. 2018) to identify HGT-derived genes that encode such secreted proteins. This analysis showed that 12 of the 44 candidates encode secreted proteins. Compared with the total number of the secreted proteins in the genome of V. dahliae strain JR2 (858 out of 11,430), HGT candidates are significantly enriched for genes that are predicted to encode secreted proteins (Fisher exact test, P < 0.0001). To know more about the putative functions of these proteins, we searched for conserved domains in each of these proteins using Interproscan (Jones et al. 2014). Among the 12 secreted proteins, 5 are glycoside hydrolases (GHs) and two are carbohydrate esterases (CEs), which are all carbohydrate-active enzymes (CAZymes). CAZymes are responsible for the synthesis and breakdown of glycoconjugates, oligo- and polysaccharides. CAZymes of plant-associated fungi usually comprise a large number of plant cell wall-degrading enzymes (Zhao et al. 2013). The remaining five are Ave1, a chondroitin AC/alginate lyase, a phosphoinositide phospholipase, an amidohydrolase, and a protein without functional annotation. Subsequently, we also searched for the functional annotation of the 30 nonsecreted proteins. Besides the previously identified glucosyltransferase and seven proteins without functional annotation, we obtained diverse predicted functions including transcription regulation, DNA repair, hydrolase, and metabolic activities. Overall, we observed that the full set of HGT candidate proteins (secreted and nonsecreted) is also enriched for GHs (10 out of 44) (Fisher exact test, P < 0.0001) when considering to the total amount of the GH genes in the genome of V. dahliae strain JR2 (265 out of 11,430).

To investigate the putative role of the HGT-derived genes in plant pathogen interactions, we compared expression patterns in planta and in vitro. It was previously shown that Ave1 is expressed in planta during interaction of V. dahliae with Nicotiana benthamiana plants (de Jonge et al. 2012). We assessed transcription (RNA-Seq) data of V. dahliae during colonization of A. thaliana and found that 12 HGT candidates showed in planta expression, one of which is Ave1. Four of these 12 were found to be induced when compared with the expression in in vitro–cultured mycelium (fig. 8). This suggests that besides Ave1, three additional HGT candidates might play a role during host colonization. Moreover, one of the three induced candidates encodes a secreted protein.

. 8.

—Pair-wise comparison of HGT candidates with differential expression in vitro and in planta. Gene expressions are depicted for Verticillium dahliae strain JR2 cultured in liquid medium and upon Arabidopsis thaliana colonization, respectively. Bars represent the mean gene expression with standard deviations. The significance of difference in gene expression was calculated using t-tests relative to a threshold (TREAT) of log 2-fold-change ≥1 (McCarthy and Smyth 2009).

Verticillium dahliae HGT Candidates Preferentially Occur at Repetitive Genomic Regions

In principle, foreign DNA can be incorporated anywhere in a recipient genome as long as it does not disrupt an essential genomic element (Husnik and Mccutcheon 2017). Several studies suggest that horizontally transferred sequences are often integrated into genomic regions that are enriched in transposable elements (Gladyshev et al. 2008; Mcnulty et al. 2010; Acuña et al. 2012; Pauchet and Heckel 2013; Husnik and Mccutcheon 2017). Previously, we observed that Ave1 is located in a lineage-specific (LS) region of V. dahliae that is enriched in repeats (de Jonge et al. 2013; Faino et al. 2016). In order to assess whether the other HGT candidates are LS and associated with repeats as well, we plotted genomic repeat densities and the genomic location of all HGT candidates. Although only Ave1 is located at a genuine LS region (Depotter et al. 2018), most of the HGT candidates (34 out of 44) reside in proximity (within 1 kb) to a repetitive sequence (fig. 9). Subsequently, we compared repeat counts within 1-kb genomic windows around each of the 44 HGT candidates to the repeat counts in 1-kb genomic windows of 44 randomly picked genes (1,000 permutations), which confirmed that the HGT candidates are significantly more often adjacent to repeats than expected by chance (P = 0.0331).

. 9.

—Genomic location of the HGT candidates. The outer lane represents the chromosomes of Verticillium dahliae strain JR2. The middle lane shows the relative position of the HGT candidates in each chromosome. The inner lane shows the repeat density. The red lines indicate the locations of the LS regions (Depotter et al. 2018).

Discussion

In this study, we systematically searched for evidence of interkingdom HGT events in the genome of V. dahliae using an AI based method. Subsequently, we verified HGT candidates using a phylogenetic approach. In general, there are two classical ways to identify HGT candidates: intrinsic and extrinsic methods. Intrinsic methods focus on patterns in the primary structure of genes and genome sequences and aim to find genes or genomic regions with composition patterns that differ significantly from the rest of the genome. For example, deviation in GC content and codon usage can be considered a sign for horizontal acquisition (Lawrence and Ochman 1997). In contrast, extrinsic methods focus on comparing similarity metrics between closely related and distant taxa. For example, when a gene from a species of interest (recipient) shows higher similarity (lower E-value, higher bit score in BLAST) to sequences from distantly related (donor) species than to genes of close relatives, this gene may have been horizontally acquired. Another frequently applied extrinsic approach is to assess the phylogenetic distribution of genes across a panel of diverse organisms. These phylogeny-based methods are generally thought to be more accurate (Poptsova and Gogarten 2007; Poptsova 2009). In principle, an HGT event will create a discrepancy between the gene and species trees. However, phylogeny-based methods are much more computationally intensive, because they require sequence alignments and phylogeny reconstructions (Dupont and Cox 2017). Thus, we chose to utilize an AI based method to pre-select candidates that can be further analyzed through phylogenetic analyses.

Although the AI method provided a rapid identification of HGT candidates, some limitations of the method should be kept in mind. The AI method relies on the quality of the reference database, the accuracy of the taxonomic assignment and the diversity of taxa covered by the reference database. The two HGT candidates (HGT-1 and HGT-2) that appear to be Verticillium-specific could also be explained by the absence of the homolog-containing species. In other words, the more closely related species are included in the database, the less Verticillium-specific HGT candidates might be identified. Furthermore, AI based methods rely on predefining ingroups and outgroups and can only identify HGT candidates that are absent in the ingroup species. In our case, the method is unable to identify candidates that were acquired before the formation of ascomycetes, and candidates that were also acquired by fungal nonascomycete species independently at more recent events from the same donors.

Nevertheless, the 42 HGT candidates that were also found in other ascomycete species were most likely acquired during ancient transfer events. The extent of gene transfer from prokaryotes to fungi is highly different across diverse ascomycete fungal lineages. A whole-genome study of HGT in Aspergillus fumigatus revealed that 40% of 189 transferred regions (containing 214 genes) were of bacterial origin (Mallet et al. 2010). In contrast, a study with a similar objective found only 11 genes with a convincing signature of transfer from bacterial origin in the genomes of three Colletotrichum species (Jaramillo et al. 2015). Three of the 11 Colletotrichum HGT candidates were also found in V. dahliae (strain VdLs17) and V. alfalfae (strain VaMs102) (Klosterman et al. 2011). In our study, we also found these three candidates (HGT-14, HGT-15, and HGT-44) in V. dahliae strain JR2.

We found that HGT candidates are enriched for genes encoding secreted proteins, many of which belong to a family of carbohydrate-active enzymes (CAZymes) or GHs. GHs play a fundamental role in the decomposition of plant biomass (Murphy et al. 2011), which suggests that bacterial genes might contribute to the plant-associated life styles of Verticillium species as well as of many other ascomycetes. Although these V. dahliae GH genes are not induced during infection of A. thaliana when compared with in vitro expression data, it is possible that the GH genes are induced during interaction with other host plants as V. dahliae is able to cause disease in hundreds of plant species (Fradin and Thomma 2006; Inderbitzin and Subbarao 2014). Alternatively, these proteins might play roles during life stages outside the host. Overall, when assessing transcription (RNA-Seq) data of V. dahliae-infected A. thaliana, we identified four HGT candidates (including Ave1) that were induced when compared with expression in in vitro–cultured mycelium. One of the induced candidates encodes a secreted protein that is predicted to be a phosphoinositide phospholipase, a key metabolic enzyme that is needed by all living organisms to hydrolyze phospholipids into fatty acids and lipophilic substances (Ghannoum 2000). Phospholipases are also signaling molecules that elicit stress tolerance and host immune responses in fungi (Ghannoum 2000; Soragni et al. 2001; Köhler et al. 2006). Phospholipases have been studied in some plant pathogenic fungi, such as Fusarium graminearum (Zhu et al. 2016) and Magnaporthe oryzae (Ramanujam and Naqvi 2010; Zhang et al. 2011; Yin et al. 2016), and the phospholipase (FgPLC1) regulates development, stress responses, and pathogenicity of Fusariumgraminearum (Zhu et al. 2016). Collectively, this may suggest that the V. dahliae phosphoinositide phospholipase plays an import role during interaction with host plants.

Conclusions

In this study, we applied a conservative approach combining an AI method and phylogenetic analysis to analyze interkingdom HGT to V.dahliae. Besides the previously identified effector gene Ave1 and a glucosyltransferase-encoding gene, we revealed additional HGT candidates, all of which are of bacterial origin, indicating a high number of interkingdom gene acquisitions.