Introduction Saprolegnia species are watermolds or oomycetes that are endemic to probably all fresh water ecosystems. These understudied pathogens can cause destructive diseases of amphibians, crustaceans, fish and insects in aquaculture and in natural environments worldwide [1], [2]. With the rise of fish as a principal source of animal protein, and the decline of wild fish stocks, aquaculture production has increased on average by 11% per year worldwide over the past ten years (FAO Fishery Information). Intensive aquatic farming practices have resulted in explosive growth in pathogen populations, which has been exacerbated by the ban of malachite green as a pesticide. Losses due to microbial, parasitic and viral infections are the largest problem in fish farms nowadays, and have a significant effect on animal welfare and sustainability of the industry. The salmon farming industry is particularly affected by Saprolegnia parasitica. This pathogen causes Saprolegniosis (also known as Saprolegniasis), a disease characterized by visible grey or white patches of mycelium on skin and fins, and subsequent penetration of mycelium into muscles and blood vessels [1], [3]. It is estimated that 10% of all hatched salmon succumb to Saprolegnia infections and losses are estimated at tens of millions of dollars annually [2]. In addition to the damage to the aquaculture industry, the declines of natural salmonid populations have also been attributed to Saprolegnia infections [1]. More in-depth knowledge of the epidemiology, biology and pathology of the pathogen is urgently needed. A draft genome sequence of S. parasitica provides an excellent starting point to investigate the disease process at the molecular and cellular level and may lead to novel avenues for sustainable control of Saprolegniosis. Animal pathogens have evolved independently multiple times in lineages such as Stramenopila, Alveolata, Amebozoa, Euglenozoa and Mycota, as well as in numerous bacterial lineages. Oomycetes such as Saprolegnia belong to the kingdom Stramenopila (Patterson, 1989), syn. Straminipila (Dick, 2001), that includes photosynthetic algae such as kelp and diatoms, ubiquitous saprotrophic flagellates such as Cafeteria roenbergensis, and obligate mammalian parasites such as Blastocystis [4], [5]. Although many Saprolegnia and related species are capable of causing diseases on a wide range of animal hosts including humans, relatively little is known about their mechanisms of pathogenicity. Among the oomycetes, most animal pathogens including S. parasitica belong to the class Saprolegniomycetidae (Figure 1A). The oomycetes also include many plant pathogens and these are mainly concentrated within the class Peronosporomycetidae. There are a small number of interesting exceptions to this otherwise sharp dichotomy, including the mammalian pathogen Pythium insidiosum (Peronosporomycetidae) and the plant pathogens Aphanomyces euteiches and Aphanomyces cochlioides (Saprolegniomycetidae) [1], [6], [7]. 10.1371/journal.pgen.1003272.g001 Figure 1 Taxonomy and ancestral genomic features in S. parasitica. (A) Animal pathogenic and plant pathogenic oomycetes reside in different taxonomic units. (B) Comparison of intron number in phytopathogenic oomycetes (the average count from the total genes of P. infestans, P. ramorum, P. sojae, Py. ultimum and H. arabidopsidis) and S. parasitica among all genes. (C) Significant difference in intron number in 4008 orthologous genes shared by S. parasitica and Phytophthora species (average intron count of P. infestans, P. sojae and P. ramorum). (Wilcoxon test, p 50%). 20% of the core set is not detectable in the S. parasitica proteome (Figure S2B). Although the genomes of the peronosporomycetes show substantial conservation of gene order (synteny), little of this synteny is preserved in S. parasitica, as was observed for A. candida [16]. Interestingly, compared to other oomycetes, S. parasitica genes contain a larger number of introns (Figure 1B, Table 1). More than 73% of the S. parasitica genes contain at least one intron, compared to 50–60% in other oomycete species (Table 1). Among 4008 orthologs shared between S. parasitica and three Phytophthora species, the majority of the genes have different numbers of introns. For example, more than half of the S. parasitica genes have more introns than their orthologs in Phytophthora, and 15% of the S. parasitica genes have 5 or more additional exons compared to their Phytophthora orthologs (Figure 1C). The intron abundance in S. parasitica potentially more closely matches the ancestral state, assuming a trend of intron reduction as found in animal and fungal lineages [34]. The S. parasitica genome has very few known mobile elements, which is consistent with its smaller size compared to the transposon-rich Phytophthora genomes. Of the 160 repeat families identified among all sequenced Phytophthora species, only one LTR retrotransposon family was found in the S. parasitica genome (Figure S3A). This group of LTR elements, which occur at low copy numbers ( 30%) is indicated by +. c Lectin and lectin-like families. S. parasitica proliferates in host tissue rich in proteins and ammonium. Concomitantly, its pathways involved in inorganic nitrogen and sulfur assimilation have degenerated (Figure 2A). The loss of these metabolic capabilities has occurred independently in the obligate oomycete plant pathogen H. arabidopsidis [15], as well as within several lineages of obligate fungal plant pathogens, presumably due to the high level of parasitic adaptation in these organisms. Strikingly the same physical clusters of genes have been lost in each lineage, namely the genes encoding nitrate reductase, nitrite reductase, sulfite reductase and nitrate transporters (Table S5). Also in line with a protein-rich environment that is a major source of both carbon and nitrogen, the S. parasitica genome contains 56 genes predicted to encode amino acid transporters. Most of the S. parasitica transporters appear to be novel because less than 20 of the predicted amino acid transporters have orthologs or closely related paralogs in other oomycete genomes. Phylogenetic analysis shows there are lineage-specific expansions of amino acid transporter genes in the different oomycete genomes, with recently duplicated S. parasitica genes forming the largest group (Figure 2B). 10.1371/journal.pgen.1003272.g002 Figure 2 Metabolic adaptations to animal pathogenesis. (A) Independent degeneration of nitrite and sulfite metabolic pathways in animal pathogens and obligate biotrophic plant pathogens. Red cross indicates the gene encoding the enzyme is absent in the genome. (B) Lineage specific expansion of amino acid transporters. Members from Pythium (black), Hyaloperonospora (green), Albugo (blue) and S. parasitica (red) are included. - The S. parasitica-specific clade is marked with red dots. (C) Secreted peptidase families in S. parasitica and phytopathogenic oomycetes (the average count from the total peptidase genes of P. infestans, P. ramorum, P. sojae, Py. ultimum and H. peronospora) . Peptidase_C1, Peptidase_S8 and Peptidase_S10 are the largest families in S. parasitica. (D) Lineage-specific expansion of peptidase_C1 family. Members from P. sojae, P. ramorum and P. infestans (black) and S. parasitica (red) are included. The S. parasitica-specific clade is marked with red dots. Other metabolic differences with oomycete plant pathogens The gene for phospholipase C (PLC) is absent in all of the sequenced peronosporomycete plant pathogens, but is present in S. parasitica (SPRG_04373). Phylogenetic analysis groups the S. parasitica PLC gene with that of other heterokont species (Figure S5, Table S6). This shows that the S. parasitica PLC is most likely to be ancestral and that the absence of PLC in other oomycetes is due to gene loss. Peronosporomycete plant pathogens are sterol auxotrophs and their genomes are missing most genes encoding enzymes involved in sterol biosynthesis [44]. In contrast, analysis of the EST collection from A. euteiches and the S. parasitica genome predicts the existence of enzymes that function in a novel sterol biosynthetic pathway [26] which has been shown to lead to the synthesis of fucosterol in A. euteiches [45]. Importantly, one of the genes SPRG_09493 encodes a CYP51 sterol-demethylase (Figure S6), a major target of antifungal chemicals that could perhaps also be used to combat Saprolegniomycetes. Candidate virulence proteins Like plant pathogens, S. parasitica presumably secretes a battery of virulence proteins to promote infection. Due to co-evolution with the host, virulence proteins are typically rapidly evolving and may appear to be unique to the species, or encoded by recently expanded gene families [17], [19]. The S. parasitica genome contains a large number of genes (11,825) that are not orthologous to any known genes in other species (Figure S2A and S2B), and many recently expanded gene families. There are at least 87 pfam domains that are either unique or show recent expansions in S. parasitica as compared to other oomycete species (Table S7). An estimated 970 proteins (Table S8) were predicted to be extracellular based on previously established bioinformatics criteria [11], [12], such as the presence of a eukaryotic signal peptide, and lack of targeting signals to organelles or membranes. Many of the expanded families appear to function at the exterior or cell surface of the pathogens, such as proteins containing CBM1 (Carbohydrate Binding Module Family I according to the CAZy database (http://www.cazy.org/; [46]), ricin B lectin, Notch domains, and also numerous peptidases. Among the proteins that are unique to S. parasitica compared to plant pathogenic oomycetes, the largest families have similarities to animal-pathogenesis-associated proteins, such as disintegrins, ricin-like galactose-binding lectins and bacterial toxin-like proteins (haemolysin E). Oomycetes contain an unusually large number of proteins with novel domain combinations, recruited from common metabolic, regulatory and signaling domains [47], [48]. S. parasitica contains in total 169 novel domain combinations that are specific to this pathogen (Table S9). As described above, some of the lineage-expanded domains such as CBM and ricin are used for novel combinations to form composite proteins. Additional domains used for novel combinations are the cytochrome p450 and tyrosinase domains. Proteins carrying S. parasitica-specific domain combinations are significantly enriched (hypergeometric test, p 50). A caveat for our analysis is that S. parasitica is the only available genome so far outside of the plant pathogenic oomycetes. There are a great variety of basal oomycetes [56], [57] that do not have genome sequence information and have not been investigated. Therefore, it could be that these HGT events could have occurred in some ancestral oomycetes . 10.1371/journal.pgen.1003272.t003 Table 3 Predicted horizontally transferred genes that may be associated with pathogenesis in Saprolegnia parasitica. Pfam function Functional Description Possible Phylogenetic origin Number Genes in the family Representative gene ID Subcellular Localizationa HGT time Estimateb Disintegrin Disintegrin proteobacteria 16 SPRG_14051 secreted recent Laminin like Associated with cell surface - 1 SPRG_08424 secreted recent CHAP CHAP domain - 7 SPRG_15528 secreted Endonuclease DNA/RNA non-specific endonuclease bacteria 6 SPRG_08128 secreted HylE Haemolysin E enterobacteria 9 SPRG_04818 membrane recent a Subcellular localization is predicted by the N-terminal signal peptide, mitochondrial targeting motif and transmembrane domains. b The time of horizontal gene transfer is estimated by the presence in other oomycetes and coding potential of a given gene. ‘a recently acquired gene’ refers to a gene occurring only in Saprolegnia and having an uncharacteristic coding potential. Several groups of extracellular enzymes were potentially acquired from bacteria; for example, the CHAP (cysteine, histidine-dependent amidohydrolases/peptidases) family and a family of secreted nucleases (Table 3). The presence of the CHAP family (pfam hit E value 4-fold expression differences (p = 4-fold difference; p 4 fold differences, p 4 fold differences, p 4 fold differences, p 4-fold differences, p 4-fold differences, p 4-fold differences, p 0.75 && DP>40 ∥ DP>500 ∥ MQ0>40 ∥ SB>−0.10. SNP calls for CBS 223.65 and VI-02736 can be retrieved from the Broad Institute Saprolegnia parasitica genome database website (http://www.broadinstitute.org/annotation/genome/Saprolegnia_parasitica/MultiDownloads.html). Loss of heterozygosity and haplotype region identification SNP calls and depth of read coverage information were parsed from the VCF file described above and analyzed in non-overlapping 5 kb windows (Figure S1). Using this information genome segments were partitioned into three groups: separated haplotypes (coverage depth ranging 20–55-fold and SNP rate = 1%). Coverage of separated haplotype regions peaks at ∼40-fold and the regions are mostly devoid of SNPs. The region corresponding to the diploid consensus exhibits ∼60-fold coverage and nearly a 3% SNP rate. The peaks in Figure S1B corresponding to the separated haplotype and consensus diploid regions are connected by a small ridge, which correspond to windows spanning boundaries between the different kinds of regions. The coverage for the diploid consensus regions is not exactly double as compared to the predicted separately assembled haplotype regions, and is less than the diploid homozygous regions; most likely this results from the relative difficulty of aligning short Illumina reads to diploid consensus sequences in the context of the high polymorphism rate observed. Individual genes located in haplotype contigs were assigned as likely allelic pairs based on SNP rate, depth of coverage, and taking into consideration best reciprocal blast matches and synteny between separately assembled haplotype contigs. Genome annotation Gene finding used both evidence-based (including EST, RNA-Seq and homology data) and ab initio methods. Gene-finding algorithms FGenesH, GeneID and GeneMark were trained for S. parasitica using existing gene and EST datasets. Then a statistical sampling of gene calls as well as genes of interest were manually curated, and the results were used to validate gene calls and fine-tune the gene caller. RNA-Seq data was incorporated into gene structure annotations using PASA [77] as described in Rhind et al. [78]. Subsequently, the annotated total gene set was subjected to Pfam domain analysis, OrthoMCL clustering analysis and KEGG metabolic pathway analysis. Expression analysis with RNA-Seq Illumina RNA-Seq data was processed as follows. Sequencing adaptors were identified and removed from reads by exact match to adaptor sequences. Reads were aligned to S. parasitica gene transcripts using Bowtie (allowing up to 2 mismatches per read, and up to 20 alignments per read). Transcript levels were calculated as FPKM (fragments per kilobase cDNA per million fragments mapped). The program EdgeR [79] was used to identify differentially expressed transcripts. Transcripts with significantly different levels (p 0.999 (P 50%) are shown in solid, while those with less similarity (between 50% and 30%) are shown with a line. Sequences with less than 30% are not shown. (PDF) Click here for additional data file. Figure S3 Mobile element comparison between S. parasitica and Phytophthora. (A) Mobile elements in S. parasitica and Phytophthora. The average copy number in P. infestans, P. sojae and P. ramorum is used as the copy number for Phytophthora. The elements are sorted based on the estimated copy number. (B) The S. parasitica element LTR-Sp1 is similar to the Copia-like family (Q572G9_PHYIN) in Phytophthora species. (C) The S. parasitica line element Line-Sp1 shows most homology with LINE elements found in fish and amphibian species (no other similar elements were found in other animal species). SwissProt protein species codes were used to name the sequences. The phylogenetic tree was constructed by using the neighbor joining method with 5000 replicates for bootstrap analysis. (PDF) Click here for additional data file. Figure S4 The expanded kinome of S. parasitica. (A) The distribution of S. parasitica kinases compared to other organisms. The kinases are named after the Standard Kinase Classification Scheme at kinase.com. TK = tyrosine kinase; TLK = TK-like; STE = STE7,11,20 family of MAP kinases; CMGC = (CDK, MAPK, GSK3 and CLK) family; CK1 = cell (casein) kinase 1 family; CAMK = Calmodulin/Calcium modulated kinase family; AGC = Protein Kinase A, G, and C families. The unclassified kinases are indicated in black. (B) S. parasitica contains a large number of protein kinases that contain trans-membrane helices. Pies are scaled to the total number of kinases in each species. (C) Kinase genes that are induced in the germinating cyst stage in S. parasitica compared to mycelia. Transcripts elevated more than four-fold relative to vegetative stages are considered to be induced. (PDF) Click here for additional data file. Figure S5 (A) Phylogram of PLCYc domains of S. parasitica PLC1 and PLCs from various organisms. For phylogenetic analysis, the PLCYc domains were determined by Smart (http://smart.embl-heidelberg.de), alignments were made and regions containing gaps were eliminated resulting in a total of 88 positions in the final dataset. The optimal tree was inferred using the Neighbor-Joining method with 5000 replicates and constructed using MEGA version 4. PLC sequences were derived from NCBI (*), JGI databases (http://genome.jgi-psf.org/,#), the Sanger Institute (http://www.genedb.org,), http://bioinformatics.psb.ugent.be,). Arabidopsis thaliana AtPLC1 (Q39032*), AtPLC2 (Q39033*), AtPLC3 (Q56W08*), AtPLC4 (Q944C1*), AtPLC5 (Q944C2*), AtPLC6 (UPI000034EE4D*), AtPLC7 (Q9LY51*), AtPLC8 (Q9STZ3*), AtPLC9 (Q6NMA7*); Aureococcus anophagefferens (Auran; 18506#); Ciona intestinalis (Cioin; XP_002129990*); Cryptosporidium parvum (Crypa; Q5CR08*); Danio rerio (Danre; XP_689964*); Ectocarpus siliculosus (Ectsi; Esi0000_0131$); Emiliania huxleyi (209393#); Fragilariopsis cylindrus (Fracy;186252#); Homo sapiens (as described by [9]*); Naegleria gruberi (Naegr; 1225#); Paralichthys olivaceus (Parol; ACA05829*); Paramecium tetraurelia, (PLC1, see [10]*); Phaeodactylum tricornutum (Phatr; 42683#), Plasmodium falciparum (Plasmo; PF10_0132@); Salmo salar (Salsa; NP_001167177*); S. parasitica: Sap-PLC (SPRG_04373#), Thalassiosira pseudonana (Thaps; 263246#), Toxoplasma gondii (Toxgo; XP_002367229*). (B) Gene structure of PLC genes. PLC is missing from other sequenced oomycete genomes, but present in S. parasitica. Multiple introns have been identified in the S. parasitica PLC gene. (PDF) Click here for additional data file. Figure S6 Sterol biosynthetic pathway inferred in S. parasitica. (A) The pathway from acetyl-CoA to lanosterol. (B) The pathway from lanosterol to zymosterol. The red box shows CYP51 sterol demethylase, a target of azole anti-fungal chemicals. (C) Pathways from zymosterol to cholesterol and fucosterol. (PDF) Click here for additional data file. Figure S7 Phylogenetic distributions of infection-related molecules. (A) Classes of infection-related molecules. Two groups of PAMPs, elicitin-like and cys-rich-family-3 proteins are present in both animal- and plant-pathogenic oomycetes (colored red). The gray dots indicate infrequent occurrences. (B) Elicitin-like proteins in S. parasitica and Phytophthora. The canonical Phytophthora and Pythium elicitins are colored green. S. parasitica elicitin-like proteins are divergent and form species-specific clades. (PDF) Click here for additional data file. Figure S8 Distribution of rates of polymorphisms. (A) Summary of SNP content across 5 kb regions for Saprolegnia CBS and N12 strains. (B) Density of SNPs according to 5 kb regions of the CBS genome. The mode for the SNP rate is 2.6%. The bulge on the left side of the distribution likely corresponds to 5 kb regions of the assembly that are mosaic between haplotype and consensus diploid, as can be seen having overlap in the distribution shown in the contour plot (Figure S1C). (C) Distribution of rates of polymorphisms between strains CBS and N12. Both heterozygous and homozogous polymorphic sites were considered across 5 kb regions of the CBS genome with Illumina reads aligned from strain N12. The mode for the %SNP was computed as 3.1%. (D) Distribution of rates of polymorphisms within strain N12. Only heterozygous sites were examined in the alignments of Illumina N12 reads to the CBS strain's genome. The mode for the %SNP was computed to be 1.7%. (PDF) Click here for additional data file. Figure S9 (A) Nucleotide substitution rate between S. parasitica strain CBS223.65 and N12. Asterisks indicate significant differences between the gene family and the core orthologs (* p<0.001; ** p<10−5) based on a non-parametric Z-test. Phytophthora data is based on the published results of Raffaele et al. (2010). GSR (Gene Sparse Region), GDR (Gene Dense Region). (B) Nucleotide substitution rate between the separated haplotypes of the strain of S. parasitica strain CBS223.65. (PDF) Click here for additional data file. Figure S10 Predicted disintegrin SPRG_14052 does not enter fish cells in vitro. (A) Amino acid sequence of fusion protein SPRG_14052_mRFP-His6. The CRxxxxxCDxxExC disintegrin motif is shaded in red. The mRFP sequence is indicated in blue, the His-tag is in green. (B) RTG-2 cells were exposed to 3 µM of mRFP, SpHtp1 or SPRG_14052_mRFP-His6 and incubated for 30 min, before photography. (PDF) Click here for additional data file. Figure S11 Stage specific gene expression detected by RNA-Seq in the total gene set. (A) Genes differentially expressed during fish cell interaction. (B) Differentially expressed genes in different life stages. (C) The correlation coefficients of pairwise comparisons between RNA-Seq data sets (p<0.001). (D) Transcript levels of a subset of disintegrin-encoding genes in various life stages of S. parasitica determined by RNAseq and qPCR. For RNAseq, the log2 value of RKPM of a gene is plotted. For qPCR, transcript levels are relative to the transcript levels of SpHtp1 in cysts and normalized against the reference gene SpTub-b encoding for tubulin. Error bars correspond to four biological replicaties. (PDF) Click here for additional data file. Text S1 Supplementary information on: phospholipid modifying enzymes and signaling enzymes, sterol metabolism, disintegrin-like proteins and supplementary methods. (DOCX) Click here for additional data file. Table S1 Saprolegnia parasitica genome assembly statistics. (DOCX) Click here for additional data file. Table S2 Assembled S. parasitica genome partitioned into classes based on coverage and polymorphisms. (DOCX) Click here for additional data file. Table S3 S. parasitica genes and SNPs percentages. (TXT) Click here for additional data file. Table S4 Chitin biosynthesis, modification and degradation in oomycetes. (DOCX) Click here for additional data file. Table S5 Gene encoding nitrogen and sulphur assimilation enzymes in oomycetes. (DOCX) Click here for additional data file. Table S6 Phospholipid modifying and signaling enzymes in Saprolegnia parasitica and other oomycetes. (DOCX) Click here for additional data file. Table S7 Unique and expanded domains in Saprolegnia parasitica proteome. (DOCX) Click here for additional data file. Table S8 Saprolegnia parasitica secretome. (XLSX) Click here for additional data file. Table S9 Protein domain combinations in Saprolegnia parasitica. (XLSX) Click here for additional data file. Table S10 Candidate effectors in Saprolegnia parasitica. (XLSX) Click here for additional data file. Table S11 Polymorphism statistics for Saprolegenia strains. (XLSX) Click here for additional data file. Table S12 Genes that are significantly induced in cyst and germinating cyst stages as compared to the mycelial tissue by RNA-seq experiment. (XLSX) Click here for additional data file. Table S13 Genes encoding lineage-expanded domains that are Iinduced in cysts or germinating cysts stage. (XLSX) Click here for additional data file. Table S14 Expression of peptidase genes in S. parasitica by RNA-seq experiments. (XLSX) Click here for additional data file.