Interspecific and host-related gene expression patterns in nematode-trapping fungi

Infection experiments

The nematodes M. hapla or H. schachtii were added to plates containing the nematode-trapping fungi A. oligospora, M. cionopagum or A. dactyloides and the infection was followed under a light microscope (Figure 1). The following five combinations of fungi and nematodes were examined: A. oligospora and M. hapla (designated Ao(Mh)), A. oligospora and H. schachtii (Ao(Hs)), A. dactyloides and M. hapla (Ad(Mh)), A. dactyloides and H. schachtii (Ad(Hs)), and M. cionopagum and H. schachtii (Mc(Hs)). There was a large difference in the infection rate, both between the nematodes and between the fungi (Table 1). Cyst nematodes (H. schachtii) showed lower susceptibility to nematode-trapping fungi than did root knot nematodes (M. hapla), which is in agreement with an earlier study by Jaffee and Muldoon [18]. This might be due to differences in the composition of the nematode cuticle [24]. Between the fungi, A. dactyloides infected M. hapla at a faster rate than A. oligospora infected M. hapla. A. dactyloides colonized H. schachtii at a faster rate than the other species, especially M. cionopagum. The observation that species with constricting rings are more efficient in capturing and killing nematodes than nematode-trapping fungi with other types of trapping structures has been made in earlier studies [16].

Characterization of the transcriptome libraries

The number of reads obtained by the 454 sequencing of the five cDNA libraries corresponding to these species combinations ranged from 70 061 to 245 264 (Table 2). Based on these sequences, three different data sets were created (Figure 2). To generate the “Highly expressed transcripts” data set, the reads were assembled into isotigs (transcripts), and low abundance reads (<5 reads), short isotigs (<100 base pairs, bp) and non-fungal sequences were removed. The number of the filtered isotigs in the five libraries varied between 1 318 and 5 140 and their average sizes varied between 1 008 to 1 237 bp. Almost all of the assembled isotigs (98.9%) had fungal matches: a few matched species from Nematoda (0.9%) and a few matched other species (0.2%). To identify the set of transcripts that was most highly expressed in each library, the reads of the isotigs were normalized using the reads per kilobase pair method. This method fitted better than the reads per kilobase per million read (RPKM) method to the expectation that most transcripts have similar relative expression abundance between samples (Additional file 1). In each library, the 500 most expressed isotigs (the “Top 500 transcripts”) were analyzed.

Librarya	Ao(Mh)	Ao(Hs)	Ad(Mh)	Ad(Hs)	Mc(Hs)
Reads
Total number of reads	114 418	70 061	226 301	191 632	245 264
Filtered readsb	97 770	58 208	183 433	141 506	206 854
Isotigs
Total number of isotigs	2 663	1 354	4 517	4 003	5 258
Fungal isotigs	2 514	1 261	4 021	3 560	4 741
Nematode isotigs	11	20	21	18	44
Others	0	5	8	4	16
Filtered isotigsc	2 634	1 318	4 428	3 926	5 140
(Average size, bp)	(1 160)	(1 008)	(1 237)	(1 128)	(1 167)
Isotigs with Pfam	1 730	856	2 317	2 504	3 071
Number of unique Pfam	1 251	707	1 555	1 592	1 847
Isotigs with UniRef50	2 485	1 239	3 916	3 472	4 601
Unique UniRef50	2 326	1 174	3 101	3 230	3 981
Mapping to Ao genome
Number of mapped reads	73 623	40 051
Number of gene modelsd	7 351	6 377

^aThe following five combinations of fungi and nematodes were characterized: A. oligospora and M. hapla (Ao(Mh)), A. oligospora and H. schachtii (Ao(Hs)), A. dactyloides and M. hapla (Ad(Mh)), A. dactyloides and H. schachtii (Ad(Hs)), and M. cionopagum and H. schachtii (Mc(Hs)).

^bNumber of reads after removal of rRNA sequences.

^cNumber of isotigs after removal of non-fungal sequences, low abundance isotigs (less than 5 reads), and short isotigs (<100 bp).

^dNumber of gene models that have ≥1 read.

Comparing transcriptomes between species with not yet sequenced genomes is challenging due to the difficulties in identifying one-to-one orthologs. To circumvent this problem, we here clustered the isotigs into UniRef50 clusters. UniRef50 clusters are based on pre-computed sequence clusters of the UniProt database that have at least 50 percent similarity and 80 percent coverage [26]. Istotigs were only grouped to UniRef50 clusters if passing a given threshold value (1e-10) and only one isotig per fungal species (displaying the highest sequence similarity) were assigned to a given UniRef50 cluster. A recent study including data from seven fungal genomes revealed that the grouping of gene sequences into UniRef50 clusters using the described procedure are in close agreement with traditional ortholog clustering methods (Canbäck et al., manuscript in preparation). Furthermore, the risk of clustering non-orthologous gene duplicates into a given UniRef50 cluster is reduced in nematode-trapping due to the rapid divergence of gene duplicates generated by repeat induced point (RIP) mutations [19, 20]. In our analyses, 15 713 of the in total 17 446 isotigs matched to 6 520 unique UniRef50 protein clusters. The reads of these putative orthologs were normalized with DESeq [27]. Based on the hypothesis that most transcripts are not differentially expressed, the analysis showed that a proper normalization was obtained using the 5% most highly expressed UniRef50 clusters (Additional file 2). This cohort (“Differentially expressed UniRef50 clusters”) contained 326 unique UniRef50 clusters.

The third data set, “Host-specific gene expression”, was generated for identifying genes that were differentially expressed due to the nematode host species. Because the genome of A. oligospora is available [19], the analysis focused on comparing the transcriptional response of this fungus when infecting M. hapla and H. schachtii. The reads from the Ao(Mh) library were mapped to 7 351 genes and those from the Ao(Hs) library to 6 377 genes (Table 2).

Divergence in gene expression

To compare the functional groups of genes that were expressed in the different libraries, the abundances and expression levels of Pfam domains in the “Top 500 transcripts” data set were analyzed. The number of Pfam domains found in the five libraries varied between 330 and 412 (Additional files 3 and 4). In total, 700 Pfam domains were found in at least one of the libraries.

A principal component analysis (PCA) of the abundances of the Pfam domains showed that the libraries mainly clustered according to the fungal species (Figure 3A). The first axis (explaining 37% of the variability) separated the two libraries of A. dactyloides (Ad(Mh) and Ad(Hs)) from those of A. oligospora (Ao(Mh) and Ao(Hs)). The second axis (27%) separated M. cionopagum (Mc(Hs)) from the two Arthrobotrys species. Clearly, in both A. dactyloides and A. oligospora, the divergence in Pfam expression patterns associated with fungal species was larger than that related to the host nematode species.

Analysis of the expression levels of the UniRef50 clusters confirmed these patterns. A PCA based on the 5% most highly expressed UniRef50 clusters showed that the first axis (49%) separated the two libraries of A. dactyloides from the two libraries of A. oligospora (Additional file 5). M. cionopagum was separated from the other libraries along the second axis (26%). Scatter plots showed that the numbers of UniRef50 clusters that were differentially expressed more than twofold between the fungal species were greater than between the host nematode species (Figure 4).

Commonly expressed transcripts

Furthermore, transcripts containing the Atg8 domain were highly expressed in all fungi. Atg8 is an essential protein in the autophagic pathway [29] and disruption of a homolog of this gene in A. oligospora leads to reduced trap formation [30]. In addition, all three fungi expressed transcripts with the Pfam domain RicinB_lectin_2 (PF14200). Ricin-B lectins are ribosome-inactivating proteins (RIPs) consisting of a catalytic A-chain and a sugar-binding B-chain [31, 32]. All fungi have a RicinB_lectin_2 transcript that match to G1X3G7 in A. oligospora. G1X3G7 is a protein with a length of 134 amino-acid residues (aa), without a secretion signal and with low sequence similarity to other RicinB lectins in the UniProt database. A RicinB_lectin_2 domain-containing protein (MOA) with nematotoxic activity against Caenorhabditis elegans has been identified in the basidiomycete Marasmius oreades[33]. The nematotoxicity was dependant on the cysteine protease activity of MOA and the binding of its lectin domain to glycosphingolipids in the worm intestine. MOA consists of 293 aa and lacks a classical secretion signal [33]. Sclerotinia sclerotiorum agglutinin (SSA) is a RicinB_lectin_2 domain-containing protein with a length more similar to the G1X3G7 protein in A. oligospora[34]. SSA has a length of 153 aa, lacks secretion signal and shows insecticidal properties when fed to the pea aphid Acyrthosiphon pisum[34].

Previous studies have shown that subtilisins (peptidase_S8) are important virulence factors in nematode-trapping fungi. In A. oligospora they have a key role in the early stages of infection, including immobilization of the captured nematode [19, 35, 36]. A. oligospora has 52 genes containing the peptidase_S8 domain [20]. However, only one transcript containing the peptidase_S8 domain was identified among the highly expressed transcripts in each library of Ao(Mh) and Ao(Hs). BLASTX searches showed that both transcripts displayed the highest sequence homology to the A. oligospora protein G1XLL2. The other cDNA libraries contained also only one transcript with the peptidase_S8 domain among the top 500 expressed genes. The three transcripts in Ad(Mh), Ad(Hs) and Mc(Hs) all displayed the highest sequence homology to H072_8474 in M. haptotylum. Interestingly, G1XLL2 and H072_8474 are orthologs (T. Meerupati, B. Canbäck, D. Ahrén, A. Tunlid, manuscript in preparation). Furthermore, G1XLL2 was the most expressed peptidase_S8 gene in A. oligospora during early infection (6 and 10 hours) of C. briggsae, and H072_8474 was the second most expressed peptidase_S8 gene in M. haptotylum during early infection (4 hours) of C. briggsae[20]. H072_8474 was also identified in the proteome of both the knob and the mycelia in M. haptotylum[22]. This shows that despite the large number of peptidase_S8 genes only a few are highly expressed during infection.

Stress proteins were highly expressed in all fungi independent of trapping structure. They included heat-shock proteins and chaperones such as DnaJ, HSP70 and HSP90; gluthatione S-transferases; and antioxidant enzymes such as thioredoxin and catalase. Antioxidants are enzymes involved in the protection of the cell from oxidative damages induced by reactive oxygen species (ROS) [37]. ROS are continuously produced in the cell as byproducts from various metabolic pathways and have an important role(s) in signaling [38]. In the plant-pathogenic fungus Magnaporthe grisea, ROS-generating NADPH oxidases (Nox1 and Nox2) are essential for pathogenicity [39]. The authors [39] suggest that the generated ROS accumulate in the appressorium to facilitate oxidative cross-linking of cell-wall proteins. This leads to a strengthening of the cell wall of the appressorium that will eventually resist high turgor pressure [39]. Transcripts with sequence similarity to the Nox proteins in M. grisea were regulated in all fungi of our study during infection. However, none of them were found among our top 500 most expressed transcripts.

Interspecific variation in gene expression

To identify the variable sets of transcripts, that is, the transcripts that were differentially regulated depending on the fungal species, the expression levels of the putative orthologs identified using the UniRef50 clusters were compared among the libraries (Figure 5; Additional file 7). Among these were transcripts encoding peptidases (peptidase_M1 and peptidase_M24), lectins (FB_lectin, B_lectin, RicinB_lectin_2), tyrosinase, transcription factors, cell-signaling components, Atg8, various stress response proteins, proteins containing the WSC domain and the DUF3129 domain. DUF3129 is a domain of unknown function that is found in the gas1 protein of M. grisea, which participates in appressorial penetration and lesion formation [40]. Interestingly, DUF3129 was highly expressed in M. cionopagum but not expressed at all in A. dactyloides. DUF3129 was identified in 12 transcripts among the top 500 most expressed transcripts in the Mc(Hs) library. During A. oligospora infections, this domain was identified in one transcript in Ao(Mh) and in four transcripts in Ao(Hs) among the top 500 most expressed transcripts. DUF3129 is an expanded gene family in nematode-trapping fungi and both M. haptotylum and A. oligospora have 33 genes encoding this domain [20]. Seventeen of these genes were previously found among the 10% most expressed genes during nematode infection by M. haptotylum, whereas only two were among the 10% most expressed genes during nematode infection by A. oligospora[20]. The DUF3129 domain is thus highly expressed during infection among the species that form adhesive branches and adhesive knobs. Further studies are needed to investigate the function of the DUF3129 domain in the nematode-trapping fungi during infection.

Transcripts encoding the fungal fruit-body lectin (FB_lectin) and the D-mannose binding lectin (B_lectin) were only highly expressed in A. oligospora and not in the other two fungi (Figure 5). In the UniRef50 cluster containing the FB_lectin domain we identified the previously studied AOL lectin (Q00233) [41]. The transcript matching to this cluster was the third most expressed of all transcripts during Ao(Mh) infection and the 25th most expressed of all transcripts during Ao(Hs) infection. Earlier studies have shown that AOL functions as a storage protein during both saprophytic and parasitic growth [42]. However, deletion of this gene did not affect the fungus’ ability to infect nematodes [43].

Previous studies have shown that proteins containing the carbohydrate-binding domain WSC comprise a large and rapidly evolving gene family in M. haptotylum[22]. Phylogenetic analysis of the 33 WSC-containing proteins in M. haptotylum revealed a clade of 15 WSC paralogs [22]. This clade contains only one (G1X6Q5) of the 16 WSC proteins identified in A. oligospora[22]. Thirteen of the 15 WSC paralogs of M. haptotylum were at least twofold upregulated during the infection of the nematode C. briggsae[20]. In this study, transcripts encoding WSC proteins were highly expressed by all fungi during infection of plant-parasitic nematodes (Table 4). The largest number of transcripts encoding WSC domain proteins was expressed by M. cionopagum. In total, 19 transcripts of WSC domain proteins were identified in the Mc(Hs) library, of which seven were found among the top 500 transcripts. Four of these transcripts displayed closest sequence homology to proteins found in the expanded clade of WSC proteins of M. haptotylum[22]. The libraries of A. oligospora and A. dactyloides contained a lower number of transcripts of WSC domain proteins. Among the top 500 transcripts, the Ao(Hs), Ao(Mh) and Ad(Mh) libraries each had two WSC domain proteins, whereas none were found in the Ad(Hs) library. Transcripts displaying highest sequence homology to the previously mentioned A. oligospora protein G1X6Q5 were identified in both the Ao(Hs) and the Ao(Mh) libraries. The deduced proteins of two highly expressed WSC transcripts in A. dactyloides did not show any sequence similarity to the proteins found in the expanded clade of paralogs in M. haptotylum[22]. Taken together, the comparative transcriptome analysis shows that the WSC domain proteins comprise a large and divergent gene family that is highly expressed during pathogenesis in nematode-trapping fungi. The specific sets of genes that are expressed depend on the fungal species and the nematodes being infected, which suggests that the function of the WSC proteins is to contribute to the specialization of the trapping mechanisms.

Virulence associated transcripts

A BLAST search of the top 500 transcripts in each library was conducted in the pathogen–host interaction protein database (PHI-base) [44]. PHI-base contains experimentally verified pathogenicity, virulence and effector genes from fungi, oomycetes and bacterial pathogens. In total, 97 unique PHI-base genes were identified. Genes with sequence similarity to seven of them were found in at least ten gene models in either M. haptotylum or A. oligospora[20]. They included RBT4 from Candida albicans, which is necessary for virulence [45]. The function of RBT4 is unknown but it contains a CAP (Cysteine-rich secretory proteins, Antigen 5 and Pathogenesis-related 1 protein) domain [45]. Among the 97 identified PHI-base genes, 15 were highly expressed by all fungal species (Additional files 8 and 9) and 82 were highly expressed by one or two of the fungal species (Additional file 10). The PHI-base genes expressed by all fungi included stress response genes and several cell signaling genes containing the Ras domain. The PHI-base genes that differ in expression between the fungi included aspartic proteases and the gas 1 and gas 2 proteins of M. grisea[40] that contains the DUF3129 domain.

Host-specific gene expression

A scatter plot of the gene expression in A. oligospora during the infection of M. hapla versus H. schachtii showed that a majority of the genes had similar expression levels (Figure 6). However, 105 genes were expressed at levels at least 5-fold higher during infection of M. hapla than during infection of H. schachtii (Additional file 11), and 65 genes were expressed at levels at least 5-fold higher in H. schachtii than in M. hapla (Additional file 12). Genes predicted to encode secreted proteins were enriched among the differentially expressed genes in both nematodes. The proportion of secreted proteins among the upregulated genes in M. hapla (Ao(Mh) or Ad(Mh)) was 12.4% (13 out of 105) and in H. schachtii (Ao(Hs) or Ad(Hs)) was 13.8% (9 out of 65). In comparison, the proportion of secreted proteins among all genes that were used for the host-specific gene expression analysis was 7.3% (304 out of 4,138).

Culture of organisms and infection experiments

Cultures of A. oligospora (ATCC 24927), M. cionopagum (CBS 220.54) and A. dactyloides (CBS 109.37) were maintained on corn meal agar 1:10. Infested soil/roots of M. hapla (Strain E 226) were obtained from Prof. Dr. Gerrit Karssen (Plant Protection Service, HC Wageningen, the Netherlands) to start the culture of this nematode. Small pieces of roots infested with M. hapla were inoculated in rhizosphere of 3-week-old tomato plants raised in the green house of Department of Biology, Lund University, for induction of root knots. After 8 weeks, the infected tomato roots with well developed knots and egg masses were gently washed under running tap water. Egg masses of M. hapla were picked by fine forceps from knots of infected roots under a stereoscopic binocular microscope and surface disinfected for 1 minute in a 0.5% NaOCl solution and rinsed three times in sterile distilled water. Egg masses were then collected in Petri dishes (30 mm) in sterile distilled water and incubated at 22 ± 1°C for 48 hours for hatching of second-stage juveniles (J₂). After incubation, freshly hatched J₂s were separated from the egg masses and collated in Eppendorf tubes, surface sterilized with 0.5% NaOCl for 2 minutes, and rinsed five times with sterilized distilled water, and used for infection experiments. J₂s of H. schachtii were obtained from HZPC in the Netherlands (http://www.hzpc.com) and used for infection experiments after sterilization and washing as described for M. hapla.

Infection experiments were performed using a dialysis membrane assay [48]. Briefly, conidia of A. oligospora, M. cionopagum and A. dactyloides were inoculated onto several pieces of sterilized dialysis membrane (spectra/por 4, Spectrumlabs). The membranes were placed over plates containing modified low-nutrient mineral salt (LNM) medium (KCl 1.0 g/l, MgSO₄ 0.2 g/l, ZnSO₄.7H₂O 0.88 mg/l, FeCl₃.6H₂O 3.0 mg/l, thiamine-HCl 0.2 mg/l, biotin 0.005 mg/l, L-phenylalanine-L-valine 0.1 g/l, agar 10 g/l, pH 6.5) [48, 49]. Infection structures (traps) were induced by adding 40–50 specimens of the nematode Panagrellus redivivus L. (Goodey) to the hyphae growing on each dialysis membrane. P. redivivus was grown axenically in a soya peptone-liver extract [50]. After several days, when substantial amount of traps have been developed and all added nematodes have been killed and digested, the infection experiments were started by adding 75–100 surface sterilized second-stage juveniles of M. hapla or H. schachtii. The following five combinations of fungi and nematodes were examined: A. oligospora and M. hapla (designated Ao(Mh)), A. oligospora and H. schachtii (Ao(Hs)), A. dactyloides and M. hapla (Ad(Mh)), A. dactyloides and H. schachtii (Ad(Hs)), and M. cionopagum and H. schachtii (Mc(Hs)). The infection was followed under a light microscope, and the number of trapped, paralyzed (i.e. with arrested movements), and colonized (hyphae growing inside the capture nematode) were counted after various time periods. For each fungal and nematode interaction, 10 replications were used. Dialysis membranes having fungal and nematode interaction of each stage were quickly transferred into liquid nitrogen and ground. Materials were collected from all infection stages (trapped, paralyzed and infected (colonized) nematodes). The ground material was stored at -80°C until use.

RNA extraction, cDNA library construction and sequencing

Total RNA was extracted from each infection stage using the RNeasy Plant Mini kit and the RLC buffer (Qiagen) and subsequently quantified using a NanoDrop 2000C spectrophotometer (Thermo Scientific). RNA integrity was inspected using a RNA 6000 Pico kit on a 2100 BioAnalyzer (Agilent). Approximately equal amounts of RNA from the three infection stages of each fungus-nematode combination were pooled. The RNA pools were concentrated by precipitation using ammonium acetate/glycogen/ethanol as described in the MicroPoly(A)Purist Kit manual (Ambion).

From total RNA, mRNA was isolated using the PolyATtract kit (Promega) according to manufacturer’s instructions. Double-stranded cDNA was synthesized using the cDNA Synthesis System (Roche Diagnostics) according to the GS FLX Titanium cDNA Rapid Library preparation protocol (454/Roche) and using adaptors with Multiplex Identifiers (MIDs) that allow for pooling of the libraries prior sequencing. Library concentration was assessed by qPCR on a Mx3005P instrument (Stratagene) and using the KAPA Library Quantification Kit - 454 Titanium (Lib-L)/Universal (Kapa Biosystems). Based on the qPCR results all libraries were pooled to contain an equal molar amount of each library. Titration and library production (aiming at 7-16% enrichment) was performed using emulsion PCR and the Lib-L kit (454/Roche). DNA-containing beads were enriched and counted using a CASY Cell Counter DT (Roche Innovatis AG), processed using aXLR70 sequencing kit (454 Life Sciences/Roche Diagnostics), and loaded onto a picotiter plate for pyrosequencing on a 454 Life Sciences Genome Sequencer FLX machine (454 Life Sciences/Roche Diagnostics). Sequencing was conducted at the Lund University Sequencing Facility (Faculty of Science).

Bioinformatic analyses

The reads obtained from the 454 sequencing were filtered, assembled and analysed according to the flowchart shown in Figure 2. Reads matching rRNA were removed using the BLASTN algorithm [25] with an E-value threshold of 1e-5 against a custom made database of rRNA sequences obtained from the 5S rRNA database [51] and the SILVA rRNA database [52]. The remaining reads for each of the five libraries were assembled separately using the GS de novo assembler 2.6 (454 Life Sciences/Roche Diagnostics) with the -cdna option. The reads were assembled into 17 785 isotigs and 10 contigs with a length > 500 bp. In the following text, both categories were referred as isotigs, i.e. transcripts. Low abundance isotigs with less than five reads and isotigs with a length shorter than 100 bp were removed. Isotigs with top hits to non-fungal species in the UniProt database [47] (the BLASTX search) [25] were also removed. The filtered dataset contained 17 446 isotigs.

The filtered isotigs from the two A. oligospora samples were mapped to the A. oligospora genome using Gmap [53] to assess the quality of the de novo isotig assemblies. In total, 3 944 out of 3 952 isotigs (including 2 634 istotigs from the Ao(Mh) library and 1 318 isotigs from the Ao(Hs) library, Table 2)) matched the genome indicating efficient filtering of non-fungal transcripts. Only 75 of the 3 944 isotigs aligned to more than one position in the genome giving a total of 4 019 genome regions aligning to the isotigs. The low number of isotigs with multiple matches indicates a low frequency of chimeric transcripts. To further investigate the quality of the transcriptome assembly, we compared the genome sequences of isotig alignments with the 11 479 predicted genes from the A. oligospora genome using the Eval software [54]. All of the 4 019 genome regions matched to the predicted A. oligospora genes. In total, 2 652 A. oligospora genes were matched giving on average 1.5 isotigs per predicted gene which suggest that some of the isotigs may represent alternative splicing forms. In total, 94.6 percent of the aligned genome regions (3 803 out of 4 019) contained both start and stop codon (i.e. considered complete by the Eval software). The mean length of the isotigs was 1 109 basepairs compared to 1 498 basepairs for the A. oligospora gene models. The difference in length may at least partly be due to alternative splicing forms where exon skipping will give shorter transcripts than the predicted gene models. The high proportion of successful matches to the A. oligospora genome as well as to its genes, the large number of complete isotigs and the long mean isotig length clearly indicate that most of the filtered isotigs have been correctly assembled into near full length transcripts.

The filtered isotigs were used to generate two different data sets (Figure 2). The first data set (“Highly expressed transcripts”) was normalized using two different approaches; reads per kilobase pair (kb), (the number of aligned reads per transcript was divided by the transcript length), and the reads per kilobase per million reads (RPKM) method [55] (Additional file 1). The isotigs were annotated based on homology using the BLASTX algorithm [25] (threshold values of 1e-10) to the UniProt sequence database [47] and proteins of M. haptotylum[20]. The isotigs were also annotated using the pfam_scan.pl tool (ftp://ftp.sanger.ac.uk/pub/databases/Pfam/Tools/) to search the Pfam-A family protein database [56] with default thresholds. Secretion signals were predicted using the SignalP 4.0 algorithm [57]. Isotigs were considered to encode putative secreted proteins if fulfilling at least one of the following three criteria: 1) Isotigs having a secretion signal in the same frame as the Pfam domain; 2) Isotigs having a secretion signal in the longest predicted open-reading frame (ORF) in the same frame as the BLASTX match (threshold value of 1e-10) to a protein in the UniProt database or protein from M. haptotylum; 3) Isotigs having a BLASTX match to a protein in the UniProt database or protein from M. haptotylum that contains a secretion signal. Orphans were identified as isotigs lacking both Pfam domains and BLASTX matches in the UniProt database and the M. haptotylum genome (threshold value of 1e-5) against a species other than itself. Orphans with a secretion signal in the longest ORF were considered to be putative secreted proteins. PCA and hierarchical clustering were performed using the Omics Explorer ver. 2.2 (Qlucore). Virulence-related genes were identified by BLASTX [25] similarity searches against the PHI-base database version 3.2 [44] using a cutoff of < 1e-10.

The second data set (“Differentially expressed UniRef50 clusters”) was obtained by matching the isotig sequences to UniRef50 clusters [26]. The UniRef50 clusters contain a representative of UniProt sequences that show 50% sequence similarity and 80% overlap with the longest sequence in the cluster. The isotigs with the highest BLASTX score to a UniRef50 sequence cluster from each library were considered as putative orthologs giving maximum one isotig from each species for any given UniRef50 cluster to take alternative splicing into account. Sometimes no isotig from one species has a significant match to a particular UniRef50 cluster (cutoff 1e-10). The transcript abundance of these putative orthologs, called “UniRef50 clusters” was normalized to correct for different library sizes using the R/Bioconductor software package DESeq 1.10.1 [27] (Additional file 2).

To identify the third data set, “Host-specific gene expression” (Figure 2), GSMapper 2.8 (454 Life Sciences/Roche Diagnostics) was used with the -cdna and -cref parameters to map the reads from the two A. oligospora libraries (Ao(Mh) and Ao(Hs)) against the coding sequences of the 11 479 genes predicted in the A. oligospora genome [19]. The read counts were normalized using DESeq [27]. A homology search of the mapped A. oligospora proteins was carried out using the BLASTP algorithm [25] to the UniProt database and proteins of M. haptotylum (threshold values of 1e-5). Secretory A. oligospora proteins were predicted using SignalP 4.0 [57].

Sequence accession numbers

Sequences can be accessed from the database http://mbio-serv2.mbioekol.lu.se/NematodeTrappingFungi/. The short read pyrosequences from Ad(Mh), Ad(Hs), Ao(Mh), Ao(Hs) and Mc(Hs) are available at NCBI SRA database with the Bioproject IDs PRJNA230433, PRJNA230458, PRJNA230459, PRJNA230446 and PRJNA230448, respectively.