Comparative Genomics and Transcriptomics Analyses Reveal Divergent Lifestyle Features of Nematode Endoparasitic Fungus Hirsutella minnesotensis

Hirsutella minnesotensis [Ophiocordycipitaceae (Hypocreales, Ascomycota)] is a dominant endoparasitic fungus by using conidia that adhere to and penetrate the secondary stage juveniles of soybean cyst nematode. Its genome was de novo sequenced and compared with five entomopathogenic fungi in the Hypocreales and three nematode-trapping fungi in the Orbiliales (Ascomycota). The genome of H. minnesotensis is 51.4 Mb and encodes 12,702 genes enriched with transposable elements up to 32%. Phylogenomic analysis revealed that H. minnesotensis was diverged from entomopathogenic fungi in Hypocreales. Genome of H. minnesotensis is similar to those of entomopathogenic fungi to have fewer genes encoding lectins for adhesion and glycoside hydrolases for cellulose degradation, but is different from those of nematode-trapping fungi to possess more genes for protein degradation, signal transduction, and secondary metabolism. Those results indicate that H. minnesotensis has evolved different mechanism for nematode endoparasitism compared with nematode-trapping fungi. Transcriptomics analyses for the time-scale parasitism revealed the upregulations of lectins, secreted proteases and the genes for biosynthesis of secondary metabolites that could be putatively involved in host surface adhesion, cuticle degradation, and host manipulation. Genome and transcriptome analyses provided comprehensive understanding of the evolution and lifestyle of nematode endoparasitism.


Introduction
Nematophagous fungi are a diverse group of carnivorous fungal species that use their specific mycelia structures or conidia to capture the vermiform nematodes, or use their hypha tips to parasitize the eggs and cysts of nematodes, or produce toxins to attack nematodes ). Two main groups of nematophagous fungi are known as nematode-trapping fungi invading nematode by producing specific trapping devices, and nematode endoparasitic fungi invading nematode by using adhesive or digesting conidia. Hirsutella minnesotensis (Ophiocordycipitaceae, Hypocreales, Ascomycota) is the representative of nematode endoparasitic fungi. This species mainly invades the secondary stage juveniles of cyst nematodes by producing conidia that adhere to, penetrate cuticle, and eventually kill the nematode ( fig. 1). The fungus consumes the nematode body content and finally grows out from the cadaver to produce new conidia for the next infection cycle (Sturhan and Schneider 1980;Liu and Chen 2000).
Hirsutella minnesotensis has been mainly isolated from soybean cyst nematode (SCN) Heterodera glycines in the United States and China Ma et al. 2005); whereas, on agar plates in laboratory conditions, it can parasitize 15 different nematode species belonging to 12 genera, including plant-parasitic, entomopathogenic, fungalfeeding, and bacteria-feeding nematodes (Liu and Chen 2001). As endoparasite, H. minnesotensis is the dominant species on SCNs in China Ma et al. 2005). High percentages of nematodes parasitized in field soils and responses to the suppressiveness of SCN in soybean monoculture system indicate that H. minnesotensis is crucial factor in regulation nematode populations in nature. Field experiments and greenhouse studies also suggest that H. minnesotensis has great potential as biological control agent (Lackey et al. 1993;Liu andChen 2001, 2005;Chen and Liu 2005;Mennan et al. 2006).
Hirsutella minnesotensis and Hirsutella rhossiliensis are the two closely related species detected from cyst nematodes and are the members of the Hypocreales. Hypocreales includes many species that also tend to associated with animals (Sung et al. 2008). These species are common and/or important entomopathogenic fungi, such as Beauvaria bassiana and Metarhizium robertsii for insect biocontrol (Gao et al. 2011;Xiao et al. 2012), Cordyceps militaris and Ophiocordyceps sinensis (anamorph: Hirsutella sinensis) for medicinal application (Zheng et al. 2011;Hu et al. 2013), and Tolypocladium inflatum, which produces the immunosuppressant drug cyclosporine (Bushley et al. 2013). However, H. minnesotensis and H. rhossiliensis described to date are unique in their capacity to attack nematodes, suggesting a relatively recent, independent origin of the nematode endoparasitic lifestyle.
Unlike entomopathogenic fungi which usually attach their spores to the external body surface of insects, or trapping fungi which capture nematodes by various trap devices developed from hyphae, H. minnesotensis and H. rhossiliensis produce adhesive spores to parasitize the nematode; however, the molecular basis of this fascinating biological phenomenon remains poorly understood. Only two genes encoding a serine protease were identified and emerged as a useful factor to penetrate the host's cuticle in H. rhossiliensis (Wang et al. 2007(Wang et al. , 2009. Comparative genomic research has become an essential approach for better understanding of fungal evolution, biology, and interactions with hosts and secondary metabolism. Although trapping fungi in the Orbiliales are phylogenetically distant from entomopathogenic and endoparasitic fungi in the Hypocreales, their origin, evolution, trap formation, and trapping mechanism have been elucidated based on the genomes of Arthrobotrys oligospora forming adhesive networks (Yang et al. 2011), Monacrosporium haptotylum forming adhesive knobs (Meerupati et al. 2013), and Drechslerella stenobrocha forming constricting rings (Liu et al. 2014). Genomic studies of fungi in the Hypocreales have revealed that gene duplication has driven host range and pathogenic capabilities of Me. robertsii and B. bassiana (Gao et al. 2011;Xiao et al. 2012). Inbreeding in O. sinensis has caused massive proliferation of retrotransposons, providing a tradeoff between the advantages of increased genetic variation independent of sexual recombination and deletion of genes dispensable for its specialized pathogenic lifestyle (Xiao et al. 2012). Similarly, the current known mycotoxin genes are lacking in the medicinal fungus C. militaris (Zheng et al. 2011) and T. inflatum has developed the cyclosporine biosynthetic gene cluster (Bushley et al. 2013).
Although a great deal is known about the ecology of nematode endoparasitism and the biology of this fungus in laboratory, relatively little is known about the mechanism of endoparasitic lifestyle or how it evolved the ability to infect nematode. To address these questions, the whole genome of H. minnesotensis was sequenced and compared with trapping fungi in Orbiliales and other related entomopathogenic fungi in the Hypocreales. The transcriptome during the infection process of Caenorhabditis elegans by H. minnesotensis was also analyzed. The genome sequence and transcriptome analyses allow a comparison of genes involved in nematode endoparasitism across a range of evolutionary distances and elucidate an understanding of how the fungus evolved to associate with nematodes.

Fungal Strains
Hirsutella minnesotensis strain 3608 was selected for genome sequencing. Hirsutella minnesotensis was originally isolated from SCN (H. glycines Ichinohe) in Minnesota, USA , and strain 3608 was isolated from SCN in Jiamusi, Heilongjiang Province, China. This strain shows commercial potential for biological controls of SCNs. Twenty single spores were isolated from this fungus culture on corn meal agar (CMA; BD) plate and cultures of each spore were tested for their parasitic ability on SCNs. The single spore culture with the highest parasitic rate was selected and grown in potato dextrose agar (BD) plate for 10 days. Genomic DNA was prepared using the CTAB/SDS/Proteinase K method (Moller et al. 1992). A Quanta 200 (FEI, Netherlands) scanning electron microscope was used to observe the parasitism of H. minnesotensis on SCNs (Luo et al. 2004).

Genome Sequencing and Assembly
Whole genome shotgun sequencing of H. minnesotensis strain 3608 was performed using the Illumina next generation sequencing technology. DNA libraries with 170-, 500-bp, 2-, and 5-kb inserts were constructed and sequenced with the Illumina Genome Analyzer at the Beijing Genomics Institute at Shenzhen. Quality control of reads was realized by FastQC. Low quality reads, reads with a certain proportion of N (10%), adapter contamination (15-bp overlap between adapter and reads), and duplication contamination were removed to obtain clean reads before assembly. The Kmer value of 51 was used for the genome assembly. The genome was sequenced to 128.3-fold coverage and assembled using SOAPdenovo ).

Transposable Elements and Repeat-Induced Point Mutation Analysis
Simple repeat rate has been estimated by Basic Local Alignment Search Tool (BLAST) against the RepeatMasker library (http://www.repeatmasker.org/). The REPET pipeline (Quesneville et al. 2005) was used to predict transposable elements (TEs) de novo. The TEdenovo program was used to detect TEs and group them in families by BLASTER using BLAST with the Repbase database (Quesneville et al. 2005). Genomic regions containing more than 60% TEs were identified using sliding genomic windows of 50 kb using an inhouse script. RIP index was determined by program RIPCAL (Hane and Oliver 2008). To analyze the effect of proximity of a gene to a repeat region, the RIP index (TpA/ApT) of these genes was calculated as a function of distance to the repeat sequences (Ohm et al. 2012). RIP index (TpA/ApT) exceeding 0.89 was considered as significant (Hane and Oliver 2008). All of the TEs were taken into account.

Orthology and Phylogenomic Analysis
OrthoMCL (http://www.orthomcl.org/) was used to construct gene families for the H. minnesotensis and 14 additional genomes (supplementary table S1, Supplementary Material online). Sequences with an e value cutoff of 1e-5 and at least 40% identity over 60% coverage were considered orthologs (Liu et al. 2014). To infer phylogenomic relationships, protein alignments were generated using MUSCLE (Edgar 2004) for each of 898 single-copy orthologs. Nonconserved regions in each multiple alignment were removed using an inhouse script. The conserved region of each single copy ortholog was concatenated into one sequence and the maximum phylogenomic tree of the 16 species was created by using RAxML with the best model LG+I+G identified using ProtTest (Stamatakis 2006). The divergence time between species was estimated by the Penalized Likelihood method (Sanderson 2002) with r8s (Taylor and Berbee 2006) version 1.8 (http://loco.biosci.arizona.edu/r8s/) calibrated against the origin of Ascomycota at 500-600 Ma (Lucking et al. 2009).

Proteases
The families of proteases were classified by full-length BLASTp against the MEROPS peptidase database (E < 1 Â 10 À4 ) (http://merops.sanger.ac.uk/). False positives were eliminated following unsuccessful searches against peptidase units and peptidase domains of MEROPS (E < 1 Â 10 À4 ) and the Pfam databases, respectively. Prediction of putative secreted proteases was then identified using a combination of the SignalP 3.0 and TargetP 1.1 servers.

Secondary Metabolism
The AntiSMASH pipeline (http://antismash.secondarymetabolites.org/) with Hidden Markov Model (HMM) signatures was used to identify and annotate putative polyketide synthase (PKS), nonribosomal peptide synthetase (NRPS), and terpene synthase (TPS) genes and gene clusters, and to predict the PKS and NRPS domain architecture in all 16 genomes. Sequences were aligned with MUSCLE. RAxML (Stamatakis 2006) was used to construct a maximal-likelihood phylogenetic tree with the LG+I+G+F model identified using ProtTest. The gene order and conservation of clusters were manually inspected with the gene cluster alignment results from AntiSMASH website. A difference of 50% in the number of accessory genes between H. minnesotensis and other fungi was considered significantly different (Martinez et al. 2012).

Transcriptome Analysis
Wild-type C. elegans Bristol strain N2 was routinely cultured on Escherichia coli strain OP50 and maintained following the protocol described by Brenner (1974), with minor modifications. The gravid adult nematodes were collected on a 35-mmaperture sieve, washed with sterile water, and lysed in a sodium hypochlorite solution to isolate the eggs. The eggs were hatched overnight at 22 C to obtain the first-stage larvae (L1). The second-stage larvae (L2) of C. elegans were then obtained as described by Xie et al. (2010). Hirsutella minnesotensis strain 3608 was cultured on CMA plates for 10 days to produce a large amount of conidia and then cocultured with 6,000 L2s for each plate. After 8, 16, and 36 h at 25 C, the TRIzol method (Rio et al. 2010) was used to extract total RNA, which was then reverse transcribed into doublestranded cDNA and sequenced using the Illumina technique. Five million reads were obtained using Paired-End (PE100) sequencing for each transcript. Functions of the expressed genes were predicted based on Funcat, GO, and Pfam databases. The program maSigPro was used to enrich the genes expressed in different patterns (Conesa et al. 2006). The program NOISeq (Mortazavi et al. 2008) was used to calculate the probability value of each differently expressed gene.

Nematode Bioassays
To test the effect of lectins on parasitism, the second juveniles of SCNs (J2) were treated with different lectins at the concentration of 100 mg/ml, including mannose binding lectin, concanavalin A lectin, Ricin B lectin, and lectin from Helix pomatia, and their different mixtures for 24 h. After the treatment, 100 J2s were added to each well of a 12-well plate and cultured with H. minnesotensis strain 3608 at 20 C for 24 h. The nematodes were washed off using sterilized water. An inverted microscope was used to determine the percentage of nematodes with attached conidia or colonized by the fungi (Xiang et al. 2007). Three replicate wells were performed for each treatment. Statistical differences between each two treatments were determined using the unpaired two-tailed Student's t-test. A P value of less than 0.05 was considered to be significant.

General Features
Whole genome shotgun sequencing of H. minnesotensis strain 3608 was generated using an Illumina Hiseq 2000 sequencer with 128.3-fold coverage and assembled using SOAPdenovo ) into 967 scaffolds (N50: 382.4 kb), with an estimated total size of 51.4 Mb (table 1). The genome size is at least 25% larger than those of nematode-trapping fungi in the Orbiliales (Yang et al. 2011) and insect pathogens (Gao et al. 2011;Xiao et al. 2012), but much less than that of O. sinensis (Hu et al. 2013) in the Hypocreales.
A total of 12,902 protein-coding genes were predicted using a combination of ab initio methods and transcript alignments. The coding regions from the predicted genes constitute 40.1% of the genome with an average gene density of 247 genes per Mb (table 1), which is fewer than most of other nematode-trapping fungi and insect pathogens (Gao et al. 2011;Yang et al. 2011;Xiao et al. 2012;Bushley et al. 2013;Meerupati et al. 2013;Liu et al. 2014). The average 1,624-bp gene length includes an average 1,494 of coding region and 130 bp of noncoding region. The GC content of the genome is 52.1%. A total of 145 transfer RNA (tRNA) genes were predicted from the assembly.

Expansion of Transposons in H. minnesotensis Genome
The richness of TEs is the most remarkable genome feature of this nematode endoparasitic fungus. The massive proliferation of TEs accounts for nearly 32% of the genome (supplementary fig. S1, Supplementary Material online), which is the result of the expansion of certain families. Overall 17,154 occurrences of TEs belonging to two classes were identified. Retrotransposons are dominated by Gypsy (7.18% of the genome) and Copia (4.60%) superfamilies of the long terminal repeat retrotransposons and the I superfamily (4.50%) of the long interspersed elements (supplementary fig. S1, Supplementary Material online). The terminal inverted repeats are the dominant superfamily belonging to Class II DNA transposons (8.10% of the genome).
Like other fungal genomes enriched in repeat sequences (Haas et al. 2009

RIP Mutation
RIP mutation is a genome defense mechanism in fungi during which duplicated sequences are mutated from CpA to TpA (Galagan and Selker 2004). A significant RIP was observed in H. minnesotensis genome by the high value of TpA/ApT (1.10) and low ratio of (CpA+TpG)/(ApC+GpT) (0.99) (supplementary fig. S3, Supplementary Material online). One homolog (HIM_04934) of DNA methyltransferase RID gene in Neurospora crassa (NCU02034) required for RIP was identified in H. minnesotensis genome. This supports the active RIP defense mechanism in H. minnesotensis.
TEs are frequently a target of RIP mutations, which also preferentially occur in the flanking regions harboring the TEs. Similarly to plant pathogens such as Fusarium spp. (Rep and Kistler 2010) and Dothideomycetes fungi (Ohm et al. 2012), based on the TpA/ApT index, genes more closely located to TEs are more likely to have a RIP signature in H. minnesotensis (supplementary fig. S4, Supplementary Material online). Protein-coding genes involved in secondary metabolism, transcription factors, and members of expanded species-specific multigene families that are overrepresented in the regions flanking TEs have undergone more RIP than other genes.

Lectins
The adhesion of nematophagous fungi to their host might be mediated by the interaction between lectins on the surface of trapping devices or adhesive spores and carbohydrate ligands on the nematode cuticle (Nordbringhertz and Mattiasson 1979). Lectins are carbohydrate-binding proteins which can be highly specific to sugar molecules (Sharon and Lis 1989). Protein-coding genes containing lectin domain are identified in nematophagous fungi and insect pathogens. Totally, nematode endoparasitic fungus H. minnesotensis has fewer (40) lectin genes than trapping fungi (average 69) (supplementary table S6, Supplementary Material online). Concanavalin A-like lectin binding to a-D-glucose and a-D-mannose is the most abundant family in nematophagous fungi genomes. Compared with trapping fungi, several families including H-type lectin and fucose-specific lectin which bind to N-acetyl-D-galactosamine and fucose, respectively, are absent in H. minnesotensis genome. In addition, trapping fungi encode a bulb-type lectin which may be used as a defense against nematodes rather than to adhere to them (Liu et al. 2014). Therefore, nematode endoparasitic fungus and trapping fungi have evolved different mechanism to adhere to their host.
Following the sugar inhibition experiment procedure (Nordbringhertz and Mattiasson 1979), we pretreated the nematodes with different lectins and their mixtures. The adhesion rate of H. minnesotensis on lectin-treated nematodes significantly decreased compared with the control without lectin treatment (P < 0.05), which supports the hypothesis that lectins might involve in the process of recognition and adhesion to the host (supplementary fig. S6, Supplementary Material online).

Proteases
The nematode cuticle as the primary barrier (Cox et al. 1981) encountered by nematode parasites such as H. minnesotensis is rich of protein, and proteases such as serine proteases have been demonstrated to be involved in the infection of nematodes by fungi (Wang et al. 2007(Wang et al. , 2009). Total of 293 putative genes encoding proteases, including different exopeptidases and endopeptidases, were identified in H. minnesotensis, showing no significant differences with those of nematodetrapping fungi (216-322), but representing fewer genes than those found in three insect pathogens Me. robertsii (401), B. bassiana (359), and C. militaris (341) (supplementary table S7, Supplementary Material online). Approximately 18% of these proteases identified in H. minnesotensis contain a signal peptide, which is more likely to be involved in pathogen-host interactions (Ohm et al. 2012), such as their host preference to cyst nematodes (Xiang et al. 2007), whereas entomopathogenic fungi encode more secreted proteases (more than 20%) (supplementary tables S7 and S8, Supplementary Material online) to assist their adaptation on more complex insect cuticles compared with nematode cuticles.
The most enriched secreted protease family in H. minnesotensis and other reference fungi are the serine proteases, which comprise approximately 50% of all secreted proteases (supplementary table S8, Supplementary Material online). Serine proteases, especially the subtilases (affiliated into two families: The subtilisin-like protease S8 and the serine-carboxyl protease S53), are involved in lethal activity and the infection process of nematodes (Yang et al. 2005;Wang et al. 2007Wang et al. , 2009. They also assist in the infection processes of entomopathogenic fungi by degrading host cuticles, providing nutrition and disabling antimicrobial peptides (Bagga et al. 2004;Gao et al. 2011). Fewer secreted subtilase coding genes were identified in H. minnesotensis compared with the insect pathogen Me. robertsii and trapping fungi A. oligospora and Mo. haptotylum. A phylogenetic tree based on the conserved enzymatic S8/S53-subtilisin/kexin/sedolisin domains (PF00082) revealed that subtilase genes in the endoparasitic fungus diverged from those of insect fungi, supporting the origin of endoparasitic species from insect fungi and their similar functions in protein degradation. Compared with insect fungi, S53 serine proteases, which adapt to the specific insect hosts coupled with S8 proteases (Muszewska et al. 2011), have been lost from nematode-trapping and endoparasitic fungi ( fig. 3, Clade A). Although they both belong to the cuticle degrading proteases, the subtilases from nematode-trapping and endoparasitic fungus are not grouped together. Most of the subtilases from trapping fungi only contain the enzymatic domain, whereas subtilases from the endoparasitic fungus contain other functional domains ( fig. 3, Clades B-E), such as an N-terminal subtilisin propeptide (proteinase inhibitor I9, PF05922), the cleavage of which activates the enzyme; or the DUF1034 domain (Domain of Unknown Function 1034, PF06280) embedded with proteinase associated (PA, PF02225) domain, which may be involved in recognition of the protein by vacuolar sorting mechanisms (Luo and Hofmann 2001). Therefore, subtilases of nematode endoparasitic fungus and trapping fungi have evolved different functional domain architectures.

Carbohydrate Active Enzymes
In plant pathogens, carbohydrate-active enzymes are required to degrade complex plant cell wall composed of cellulose, hemicelluloses, and/or pectin (O'Connell et al. 2012). Hirsutella minnesotensis has fewer glycoside hydrolases (155), which resembles insect pathogens (average 155) but is significantly less than plant pathogens (average 260, P = 0.0203) (supplementary table S9, Supplementary Material online). Similar to insect pathogens, several cellulase families including GH6, GH7, GH12, GH45, and GH61 are contracted or even absent in the H. minnesotensis genome. Hirsutella minnesotensis as well as insect pathogens also lack several hemicellulase families like GH10, GH11, GH26, GH29, and GH67. With fewer cellulases and hemicellulases, carbohydrate-binding module 1 protein domains enhancing the localization of CAZy enzymes on the surface of crystalline cellulose (Liu et al. 2014) are absent in H. minnesotensis (supplementary table S10, Supplementary Material online) genome. In comparison with plant pathogens, H. minnesotensis has much fewer putative cutinases (2 vs. average 10) and pectin lyases (2 vs. average 10) (table 2) that are known to be virulence factors for plant pathogens (Xiao et al. 2012). Chitinases are involved in cuticle degradation during insect and nematode infection (Xiao et al. 2012;Liu et al. 2014). Similar to insect pathogens, H. minnesotensis encodes more GH18 chitinases (23) involved in cuticle degradation than nematode-trapping fungi (average 11) and plant pathogens (16) (supplementary tables S9 and S11, Supplementary Material online).

Signal Transduction
Fungal GPCRs are key components of signaling pathways for pheromone/nutrient sensing and host recognition (Xue et al. 2008;Gao et al. 2011). Most GPCRs encoded by H. minnesotensis resemble the rice-blast fungus Magnaporthe oryzae Pth11-like proteins which are required for sensing the host surface stimuli to produce appressorium in plant pathogens (DeZwaan et al. 1999 3.-Phylogenetic tree of nematophagous and insect fungal proteinases containing peptidase S8/S53-subtilisin/kexin/sedolisin domains annotated by Pfam. Proteinase genes from one nematode endoparasitic fungus, three nematode-trapping fungi, and five insect fungi are labeled in branch with red, green, and blue color, respectively. Phylogeny was estimated using RAxML (see Materials and Methods). Five clades were identified according to the different domain architecture of proteinase genes. Genes in Clade A contain peptidase S53 domains in combination with the peptidase S8/S53-subtilisin/kexin/ sedolisin domains. This combination is found mainly in insect fungi. Most of proteinase genes in nematode endoparasitic fungus are located in Clades B-E. (Xue et al. 1998), is also predicted in H. minnesotensis and nematode-trapping fungi, whereas most insect pathogens (except Me. robertsii) lack this GPR1-like GPCR (supplementary table S12, Supplementary Material online).
Reversible protein phosphorylation by protein kinase plays critical roles in the regulation of fungal growth, pathogenesis, stress responses, toxin production, and sexual reproduction (Wang et al. 2011). The nematode endoparasitic fungus H. minnesotensis encodes more protein kinases (139) than the trapping fungi (average 115) (supplementary table S13, Supplementary Material online). Especially, serine/argininerich protein kinases involved in the modulation of salt tolerance, membrane potential, and potassium transporter in S. cerevisiae (Forment et al. 2002) are expanded in H. minnesotensis genomes (6) compared with those of trapping fungi (average 1) (supplementary table S13, Supplementary Material online).
Physiological responses following the signal transduction are regulated by activation of different transcription factors.
Hirsutella minnesotensis encodes large number of transcription factors (979) more than nematode-trapping fungi (average 374) and insect pathogens (average 422) (supplementary table S14, Supplementary Material online). The most of H. minnesotensis gene expansion involves in Zn2/Cys6, centromere protein B, homeodomain-like, CCHC zinc finger, and bZIP transcription factors. These families especially the Zn2/ Cys6 type play important role in regulating protein and polysaccharide degradation (Pel et al. 2007). Overall, compared with nematode-trapping fungi and insect pathogens, H. minnesotensis has evolved different signal controls for response to the environment stimuli.

Expansion and Diversity of Secondary Metabolic Genes
Genes encoding secondary metabolic enzymes such as PKS, the NRPS, and the TPS are significantly expanded in H. minnesotensis genome (supplementary table S15, Supplementary Material online, P < 0.05). Totally, 101 core genes were identified, compared with 10-19 in trapping fungi and 28-70 in

FIG. 3.-Continued
They have more close relationship with those of insect fungi than with trapping fungi. Clades B and C have longer peptidase S8/S53-subtilisin/kexin/sedolisin domains compared with other clades and contain DUF1034 C-terminal and proprotein convertase P domains closest to the conserved domain, respectively. Specially nine genes in Clade B contain PA signature (protease-associated domain) embedded within the peptidase S8/S53 domain. Clade D may have arisen differently; as the majority of proteins belong to nematode-trapping fungi do not have any other domains in combination with the peptidase S8/S53 domain. Half of the genes in Clade D do not carry a signal peptide, thus they do not seem to function extracellularly. Clade E is an expanded group of genes that contain signal peptide and proteinase inhibitor I9 domain on the N-terminal of proteins. insect fungi (supplementary table S15, Supplementary Material online). The numbers of Type I PKS biosynthetic genes possessed by H. minnesotensis (26) are significantly more than trapping fungi (average 2, P = 0.0151) and insect pathogens (average 13, P = 0.0014). Hirsutella minnesotensis also encodes abundant of putative NRPS genes (23) in comparison with trapping fungi (average 2) and insect pathogens (average 12). Phylogenetic analyses indicate the PKS genes of nematode endoparasitic fungus evolved from those of insect pathogens (supplementary fig. S7, Supplementary Material online). Three, six, five, and ten PKS genes are clustered with other genes involved in the biosynthesis of pigments, mycotoxins, lovastatin, and alternapyrone, respectively (supplementary fig. S7, Supplementary Material online), indicating that H. minnesotensis might produce secondary metabolites with similar structures. Six NRPS and two terpene genes are unique to the nematode endoparasitic fungus, suggesting lineage-specific expansion of these families in H. minnesotensis genome.
The core genes for secondary metabolite production in fungal genomes are typically located in clusters with other accessory genes, such as those encoding dehydrogenases, methyltransferases, acetyl transferases, prenyltransferases, oxidoreductases, and cytochrome P450s (Gao et al. 2011). In total, 94 secondary metabolite gene clusters were identified in H. minnesotensis genome. Although most PKS and NRPS genes are orthologous across the nematophagous fungi and insect fungi, the accessory genes in the secondary metabolite gene clusters are more dynamic. Only five secondary metabolism gene clusters are shared between H. minnesotensis and other insect pathogens with limited synteny (fig. 4A). This cluster diversity seems to have resulted from gene duplication or loss related to TEs associated with the secondary metabolism gene clusters ( fig. 4B). Thirty-two of 94 gene clusters in H. minnesotensis are located in TE-rich regions (>60% TE).

Transcriptomics Analysis during Infection of H. minnesotensis
Hirsutella minnesotensis was investigated about the infection process on nematodes by the RNA-seq (RNA sequencing) approach. The relative expression of all H. minnesotensis genes was determined by Illumina sequencing across four infection stages, including i) H. minnesotensis cultures on CMA medium free of nematodes to produce large amount of spores (0 h), ii) spore adhesive stage after 8 h of culture with nematode challenge (8 h), iii) infection peg formation stage after 16 h of nematode challenge (16 h), and iv) stage of fungal mycelium growth in the body of nematode after 36 h (36 h). In total, 9,277 of 12,702 predicted genes were expressed during the infection process, representing 73.0% of the H. minnesotensis transcriptome. After infection on nematodes for 8, 16 and 36 h, 1,433, 1,448 and 973 putative protein-coding genes were significantly up-or downregulated (P < 0.05), respectively, compared with that of the culture on CMA without nematodes. Many of these differentially expressed genes encode known proteins (approximately 70%), a significantly greater proportion than for nondifferentially expressed genes.
GO terms and Pfam domains were used to assign genes to functional categories and group the genes by parasitic dynamics using the K-Means clustering algorithm. Eight main clusters ( fig. 5A, clusters 1-8) accounted for about 93% of the differentially expressed genes in the four stages. Most of the blue bins showed enrichment for particular clusters of gene expression ( fig. 5B). Genes that encode proteins involved in translation, cytoplasm, respiration activity, hydrolase activity, proteasome, transmembrane transport, lipid metabolism, and organic substance transport are greatly enriched in clusters 1 and 2, representing genes that are expressed at the highest levels at the spore adhesion stage. High expression of genes responsible for the proteasome and hydrolase activity indicated their important role in cuticle degradation. Genes that showed peak expression during the period when fungal mycelium growth within the nematode body included those that are required for phospholipid metabolism, macromolecular metabolism, amino acid metabolism, protein degradation, serine-type peptidase activity, and secondary metabolism, suggesting that the fungal transcriptome undergoes detoxification in host body, degrades nematode protein substrates, and involves nutrient uptake.
Recognition and adhesion are the first two steps of fungal invasion of nematodes. In the genome of H. minnesotensis, two concanavalin A-like lectins of 18 secreted lectin proteins were upregulated at the adhesive stages (supplementary table S16, Supplementary Material online). In addition, H. minnesotensis encodes 16 adhesive proteins containing a CFEM (eight cysteine-containing extracellular membrane domain, PF05730) that are used by pathogenic fungi to recognize and attach to their hosts (Kulkarni et al. 2003), and four of them (HIM_02815, HIM_03774, HIM_04712, and HIM_ 08135) were significantly upregulated (log2-fold change > 2, P < 0.05) at the adhesive stage (8 h) (supplementary table S16, Supplementary Material online).
After the interaction of the fungus and the nematode during the peg stage, the fungus conquers the nematode defense and consumes the nematode body content for growth. Both the endopeptidases M4 and M43, and exopeptidases M14 and M28 (zinc-metallopeptidases) were highly expressed at this stage (36 h). In addition, one general substrate transporter (HIM_03796), two oligopeptide transporters (HIM_04215, HIM_06101), and one amino acid transporter (HIM_07163) were also strongly activated to facilitate the nutrient transportation. The transcriptome profile of H. minnesotensis also indicated that the enriched genes for secondary metabolites were upregulated: 17, 6, and 2 different clusters were induced before penetration, during penetration, and at the mycelium growth stage, respectively Like the infection structures of appressoria for entomopathogenic fungi (Gao et al. 2011) and the constricting ring for nematode-trapping fungus (Liu et al. 2014), cell remodeling by the protein kinase C (PKC) signal transduction pathway is important for peg formation in H. minnesotensis. The PKC1 (HIM_00683), Bck1 (HIM_03632, MAPKKK), and Mkk1 (HIM_07373, MAPKK) genes were relatively highly expressed at 16h compared with 8h (supplementary table S17, Supplementary Material online). Adaptation to stress (osmotic, oxidative, acid, and heat) mainly occurs through the high osmolarity glycerol (HOG) pathway (Romá n et al. 2007). The MAPK Hog1 (HIM_00488) gene was also relatively highly activated at the peg formation stage (supplementary table S17, Supplementary Material online). The significant upregulation of protein kinase A (PKA) during the last stage revealed the importance of the cyclic AMP-protein kinase A pathway for filamentation (supplementary table S17, Supplementary Material online). Meanwhile, two genes (HIM_01527 and HIM_07761) resembling the STM1-like GPCR responsible for triggering adaptation to nitrogen starvation were upregulated during the whole infection process. Overall, the nematode endoparasitic fungus is able to respond to diverse host-derived stimuli using distinct or shared signal pathways involving mitogen-activated protein kinase and cAMP-PKA to carry out their infection cycle on the nematode hosts (supplementary fig. S10, Supplementary Material online).

Discussion
Hirsutella minnesotensis is one of the most important endoparasites of SCNs. Different from other animal pathogens such as nematode-trapping fungi and entomopathogenic fungi, H. minnesotensis produces adhesive spores to adhere to, penetrate, and eventually kill their nematode host. The analysis of H. minnesotensis genome compared with genomes of other insect pathogens in Hypocreales and trapping fungi in Orbiliales and the transcriptome of H. minnesotensis on C. elegans have provided insights into the evolution and life strategy of nematode endoparasites.
The most significant characteristic of H. minnesotensis genome compared with other animal pathogen genomes (except O. sinensis) is that H. minnesotensis contains a strikingly rich and diverse population of TEs, associated with its relatively larger genome size. An expanded TE-rich genome is also a feature of some biotrophic obligate plant pathogens, such as Blumeria graminis (Spanu et al. 2010) and Puccinia graminis (Duplessis et al. 2011), the ectomycorrhizal fungus Tuber melanosporum (Martin et al. 2010), and the caterpillar fungus O. sinensis (Hu et al. 2013) that are closely associated with their hosts. The disruption of gene collinearity in the TErich regions ( fig. 4) increased recombination to mediate chromosomal rearrangements, deletions, and duplications, as well as greater sequence diversity (Raffaele and Kamoun 2012). A positive correlation between recombination rates and the distribution of TEs in clusters that break synteny between Ma. oryzae and other fungi has also been noted, indicating increased rates of structural variation in TE-rich regions (Thon et al. 2006). Overrepresentation of secondary metabolism genes, zinc finger transcription factors, and species-specific genes in these TE-rich regions of H. minnesotensis indicates the proliferation and mobility of TEs contribute to the duplication, expansion, and diversity of these gene families. As in other fungi the TE-rich regions were analyzed to contain genes associated with host adaptations (Spanu et al. 2010;Duplessis et al. 2011;Klosterman et al. 2011), the proliferation of TEs in H. minnesotensis might be responsible for the close association of H. minnesotensis with its host nematodes and could explain the obligate parasitism of nematode in nature (Sturhan and Schneider 1980;Jaffee and Zehr 1985;Jaffee et al. 1991;Timper and Brodie 1994;Velvis 1995;Chen and Reese 1999). Repetitive sequences in fungal genomes appeared to be the target of RIP mutations, and RIP mutations can leak to regions flanking the repetitive sequences. This characteristic in H. minnesotensis genome is consistent with the plant pathogens such as Fusarium spp. (Rep and Kistler 2010) and Dothideomycetes fungi (Ohm et al. 2012) in which higher mutation rate associated with the RIP process accelerates the evolution of effector genes that locate in repeat-rich regions. (Rouxel et al. 2011). The potential evolutionary benefit of colonization of TEs and genes encoding secondary metabolites and transcription factors in H. minnesotensis is a higher rate of mutation caused by RIP, which in turn may lead to a higher rate of evolution. This would allow the fungus to adapt more quickly to the host nematode's defenses and ecological niches, and thereby provide an advantage in the "arms race" against their hosts (Jiang XZ, Shu C and Liu XZ, unpublished data).
Hirsutella minnesotensis is a member of the order Hypocreales which is an important group of fungi that has various lifestyles, including plant pathogens, endophytes, hyperparasites, insect pathogens, and nematode parasites. The phylogenomic analyses demonstrated that H. minnesotensis is evolutionary close to the insect fungi T. inflatum and O. sinensis in the family Ophiocordycipitaceae and evolved into nematode endoparasite independently out of nematode-trapping fungi lineage ( fig. 2). This evidence supported the viewpoint that this species has divergent after fungi adapted to insect hosts. This inference is strengthened by the consistent existence of genes for cuticle degrading enzymes and contraction of genes for cellulose/hemicellulose degrading enzymes. As both insects and nematodes are rich of protein and chitin, serine proteases and chitinases involved in cuticle degradation are also abundant in nematode endoparasitic fungus and are most closely related to those genes of insect pathogens. Both H. minnesotensis and insect pathogens lost genes involved in cellulose and hemicellulose degradation, indicating that the nematode endoparasitic fungus shares more similarities with insect pathogens than with plant pathogens or saprophytic fungi. Although nematodes are the hosts of both trapping fungi and endoparasitic fungi, trapping fungi are considered to be divergent from saprophytic fungi which is attributed to the relative abundance of enzymes for cellulose degradation (Yang et al. 2012;Liu et al. 2014).
As H. minnesotensis and insect pathogens are closely related and evolved independently of nematode-trapping fungi, we assume that the nematode endoparasitic fungus and trapping fungi have evolved different pathogenesis mechanisms. Our comparative genomics analyses found that H. minnesotensis contains fewer lectin genes than trapping fungi, especially than the adhesive networks producer A. oligospora and the adhesive knobs producer Mo. haptotylum. The interaction between lectins on trapping surface and carbohydrate ligands on nematode surface was considered to mediate the adhesion of nematodes to trapping fungi (Nordbringhertz and Mattiasson 1979). This adhesion was also indicated to exist between H. minnesotensis spores and nematodes by lectin inhibition experiment in our study. Fungi producing adhesive device types such as A. oligospora and Mo. haptotylum might need strong adhesion by lectins to prevent their prey nematode from escaping. However, for endoparasitic fungus H. minnesotensis, once its spore adheres to the host cuticle, it will fall from the conidiophores, move with the host, and eventually infect the nematode. Therefore, H. minnesotensis might not require many lectins for adhesion. Drechslerella stenobrocha with constricting-ring mechanical trapping mechanism also encodes fewer lectins (Liu et al. 2014). These differences of adhesion to host are also reflected in the contraction of several cell surface proteins including carbohydrate-binding WSC domain containing proteins and GLEYA proteins in nematode endoparasitic fungus compared with adhesive trapping fungi. Both WSC and GLEYA are carbohydrate-binding domains found in fungal adhesins (Linder and Gustafsson 2008;Andersson et al. 2013).
In H. minnesotensis, genes encoding plant cell wall degrading enzymes are reduced in number compared with plant pathogens and nematode-trapping fungi. Compared with A. oligospora and Mo. haptotylum, for example, H. minnesotensis exhibited fewer or no genes encoding enzymes to degrade celluloses, hemicelluloses, cutin or pectin. Coupled with these losses in plant degrading enzymes, there are abundance of protease families in H. minnesotensis relative to trapping fungi. Serine proteases especially the subtilase type including subtilisin-like proteinases (S8) and serine-carboxyl proteinases (S53) involve in the penetration and colonization on hosts for nematophagous fungi and entomopathogenic fungi (Li et al. 2010;Muszewska et al. 2011). Consistent with previous studies (Li et al. 2010), the higher homology is found between H. minnesotensis-predicted subtilase genes and those of entomopathogenic fungi compared with those of nematode-trapping fungi, which also reflects the differences from trapping fungi in pathogenesis mechanisms. Similar to entomopathogenic fungi, subtilases in H. minnesotensis contain the enzymatic domain combined with other functional domains, whereas most subtilases in trapping fungi only contain the core domain. Nematode endoparasitic fungus does not produce trapping devices; therefore, it is likely to rely mainly on extracellular enzymes, including the subtilases (Li et al. 2010). Thus, the variety of subtilases with different functional domains in the endoparasitic fungus may facilitate the penetration and digestion of nematode cuticles. However, serinecarboxyl proteinases containing S53 domains that can couple with S8 proteases (Muszewska et al. 2011) to degrade insectsource proteins are lost in H. minnesotensis and trapping fungi, which might limit their degradative ability on insect cuticles. Thus, H. minnesotensis adapts to nematodes as host rather than insects by protease gain and loss from its relative entomopathogenic fungi. Moreover, H. minnesotensis has greater capacity of cuticle degradation than trapping fungi by encoding more chitinase to degrade chitin, a component of nematode cuticle (Liu et al. 2014). Overall, these changes in plant cell wall degrading enzymes and serine proteases suggest that H. minnesotensis is not typical saprophytic species in that it maintains a close association with protein-rich nematodes.
The H. minnesotensis genome contains very rich repertories of secondary metabolism genes compared with entomopathogenic and trapping fungi. The gene clusters are very diverse, with only a few clusters in H. minnesotensis showing limited synteny with other fungi. As each secondary metabolism gene cluster is probably involved in the biosynthesis of a specific metabolite (Collemare et al. 2008), this species can be expected to produce large and diverse spectra of secondary metabolites, some could be novel bioactive molecules. Secondary metabolism genes in insect pathogens are responsible for the production of certain insect toxins such as destruxin (Wang et al. 2012), beauvericin (Xu et al. 2008), and bassianolide (Xu et al. 2009) and determine host range (Gao et al. 2011). Secondary metabolite genes that are induced at the adhesive stage may synthesize small molecules for signaling and host manipulation, similar to effectors produced by plant pathogenic fungi (O'Connell et al. 2012), whereas those highly expressed at the peg formation and mycelium growth stages may produce toxins that help to counter host defenses and nutrient absorption (Roberts and St Leger 2004). Nematodes are very tiny animals occurring in the complicated soil habitat. After infection, the nematodes are too small to support H. minnesotensis to conquer the soil fungistasis and other fungal and bacterial competitors (Garbeva et al. 2011). Several antibiotics are easily detected, extracted, and identified from H. minnesotensis cultures (data unpublished), including the wide-spectrum antibiotic phomalactone (5,6dihydro-5-hydroxy-6-prop-2-enyl-2H-pyran-2-one) actively against entomopathogenic fungi, plant pathogenic fungi, insects and plant (Krasnoff and Gupta 1994;Kim et al. 2001), the insect toxin musacin D (Burkhardt et al. 1996;Schneider et al. 1996), and the phytotoxin musacin E [dihydro-5-(1-hydroxy-2-butenyl)-2(3H)-furanone] (Trisuwan et al. 2009). So we hypothesize that secondary metabolites, particularly the antibiotics, help the fungus survive and colonize in soil by inhibiting microbial competitors. Confirmation of these secondary metabolites and their functions will be addressed in future researches.
The lifestyle of H. minnesotensis would require a different signaling machinery to respond to a variety of environmental stimuli. In this light, H. minnesotensis as well as trapping fungi such as A. oligospora and Mo. haptotylum contain much more numbers of Pth11-like G-proteins, which have been demonstrated to involve in physiological processes in the plant pathogen Ma. oryzae (DeZwaan et al. 1999), than entomopathogenic fungi (except Me. robertsii). Hirsutella minnesotensis and trapping fungi also encode GPR1 nutrient sensors, which are absent in entomopathogenic fungi. These similar expansion of GPCRs associated with functions necessary for nutrient sense and host recognition reflects convergent evolution between nematode endoparasitic fungus and trapping fungi. However, H. minnesotensis encodes expanded protein kinase involved in regulation of fungal growth, pathogenesis, toxin production, and sexual reproduction compared with trapping fungi especially has a complex signal transduction cascade controlling stress response. Moreover, transcription factors belonging to Zn2/Cys6 fungal type, zinc finger, and bZIP which are expressed during pathogenesis (Larriba et al. 2014) are well represented in the H. minnesotensis genome, indicating that H. minnesotensis might make different physiological responses regulated by activation of different transcription factors. These analyses highlight the importance of identification of signal and transduction pathways involved in nematode endoparasitism.
The transcriptomes during nematode infection process demonstrate that H. minnesotensis finely expresses and regulates genes to finish their infection cycle on their host. At adhesion stage, two upregulated concanavalin A-like lectins were involved in adhesion between nematode and endoparasitic fungus. Several upregulated serine proteases were responsible for protein-rich cuticle degradation during fungal penetration. Subsequently, similar to entomopathogenic fungi (Gao et al. 2011) and trapping fungus D. stenobrocha (Liu et al. 2014), cell remodeling by the PKC signal transduction pathway was activated for peg formation. The high osmolarity glycerol (HOG) pathway was induced under the osmotic and oxidative stress from the nematode host body, with upregulation of several thioredoxin. After fungal mycelium growth in the nematode body, exo-and endo-metallopeptidases with potential roles in pathogenicity (Martinez et al. 2012) were significantly expressed, indicating their important role in efficient degradation of nematode protein substrates for nutrients. The cAMP-PKA pathway may be responsible for filamentation in the host body. Thus, comparative genomics and transcriptomics revealed key genes and putative pathways active during nematode endoparasitism and provided new insights into the life strategy of this fungus.
In conclusion, we sequenced the first genome of nematode endoparasitic fungus and compared with those of other animal pathogens. Phylogenomic analysis revealed that H. minnesotensis evolved after the divergence of insect pathogens which also supported by the high similarity of serine proteases and contraction of plant cell wall degrading enzymes between H. minnesotensis and insect pathogens. Compared with trapping fungi, H. minnesotensis genome contains reduced number of cellulose/hemicellulose degrading enzymes and fewer lectins for adhesion, but possesses expanded gene inventories for protein degradation by multiple functional domains-subtilases, signal transduction for nutrient sensing and host recognition and rich gene repertories of secondary metabolism. The upregulation of lectins, serine proteases, and secondary metabolism genes during the infection of nematodes indicated that these genes are involved in adhesion, recognition, cuticle degradation, signaling, and detoxification. The comparative genomics and transcriptomics provide comprehensive understanding of the evolution and life strategy of this species and highlight the future investigation of the molecular mechanism underlying this fascinating life style of nematode endoparasitic fungi.