Aphanomyces spp. genome sequencing and phylogeny relationship

Genome sequences and assemblies of the crayfish pathogen A. astaci and the opportunistic pathogen A. stellatus were generated, giving an estimated size of 45 Mb and 62 Mb (Table 1) with 16,479 and 25,573 predicted genes, respectively. As expected, genome assemblies of A. astaci and A. stellatus compiled only from short reads generated by Illumina are more fragmented and must be considered as drafts.

These data were used to perform a phylogenetic analysis with a gene set recently used to estimate timescale evolution of oomycetes [38] (see Methods for further details). Aphanomyces species form a group which diverged from other Saprolegniales (Saprolegnia parasitica, Achlya hypogena, Thraustotheca clavata) more than 100 Mya (Fig. 1). This analysis also shows that divergence times are very deep within the Aphanomyces genus, with the three Aphanomyces species diverging from each other more than 50 Mya.

Transposable Elements (TE)

TEs are known to play a prevalent role in the evolution of eukaryotic pathogens [2, 3]. A de novo characterization of TEs was thus performed on all Aphanomyces genomes. The aim of the TE analysis was to uncover and annotate all TE copies in all three Aphanomyces genomes to facilitate future studies that will investigate the origin, evolution, and genomic impact of these elements in more detail. The fraction of Aphanomyces genomes occupied by TEs (5–13%; Additional file 1: ST1d) is lower than that for all other sequenced oomycetes (17–74%) in which TEs have been mined, such as Phytophthora species [32, 39], the downy mildew Plasmopara viticola [35], and the white rust Albugo laibachii [34], with the exception of Pythium ultimum (5%) [40]. It is noteworthy that in oomycetes genome size is strongly correlated with TE content (R² = 0.85), supporting the hypothesis that TEs are a major determinant of genome size in these taxa [32]. In all three Aphanomyces species, the most abundant TEs are DNA transposons of the Tc1-mariner superfamily (up to 2.5% of the genome in A. astaci). Interestingly, none of the TEs that we identified are shared between all three Aphanomyces species, and very few TEs are shared between two species: two are shared between A. astaci and A. euteiches (74 and 83% identity over 741 and 170 bp respectively), two are shared between A. stellatus and A. euteiches (84 and 86% identity over 334 and 1262 bp respectively), and four are shared between A. astaci and A. stellatus (between 78% identity over 1820 bp and 84% identity over 195 bp) (Fig. 2). Furthermore, BLASTN searches using all TEs identified in Aphanomyces spp. as queries against all Repbase TEs did not reveal any significant hit. Together, these results suggest that Aphanomyces species TE families are characterized by a high rate of turnover.

Comparative analysis of A. euteiches and A. astaci proteomes

To identify conserved and specific features of Aphanomyces proteomes, an OrthoMCL analysis [41] was performed using A. euteiches and A. astaci proteomes and nine deep-sequenced oomycete proteomes (S. parasitica, P. infestans, P. sojae, P. parasitica, P. ramorum, Hyaloperonospora arabidopsidis, A. laibachii, Py. ultimum, and Py. irregulare) (Additional file 1: ST1e). A total of 2296 orthologous groups (34,404 genes) were detected in the 11 genomes defining a “core proteome” set (Additional file 1: ST1f). In this set 104 groups corresponding to 1528 genes defined the “core secretome” (Additional file 1: ST1f). A focus on the A. euteiches secretome (1240 genes, ~ 6% of the proteome, Additional file 1: ST1g) shows that 70% of the sequences harboring a Gene Ontology (GO) term are related to the “hydrolase activity” category (GO:0016787) (Additional file 1: ST1h and Fig. 3a). This category is enriched in enzymes with glycosyl hydrolase or peptidase functions and other typical oomycete secreted proteins, such as protease inhibitors. The secretome of A. astaci is predicted to contain 744 genes (< 5% of the proteome). While this set is probably not complete, 323 sequences harbor a GO annotation and more than 65% (217 sequences) fall into the “hydrolase activity” category with a predominance of “peptidase GO:0008233” function (160 sequences, 73%). A closer examination of both secretomes showed that numerous genes detected in A. euteiches are not reported in A. astaci (Fig. 3a). In addition, about 72% of the A. euteiches-specific secretome (i.e., sequences not present in other oomycetes including A. astaci, 506 genes) did not display a putative functional domain (368 genes), and 80% encoded proteins below 300 amino residues in size (296/368, Additional file 1: ST1i, Fig. 3b). The same observation was made for the A. astaci-specific secretome, where 160 specific genes did not harbor a functional annotation and 86% (138/160, Additional file 1: ST1i) coded for small proteins (Fig. 3b). Interestingly, 59% of these small proteins were predicted to be putative effectors from A. euteiches or A. astaci using EffectorP (Additional file 1: ST1i), a software based on properties shared by fungal effectors [42]. This suggests that Aphanomyces species presented a set of SSPs similar to fungal effectors.

Carbohydrate-Active enZymes (CAZy) and carbohydrate-binding modules (CBMs)

Since comparative analyses suggested that the repertoires of glycosyl hydrolases are divergent between A. euteiches and A. astaci, a global prediction of CAZy domains (glycosyl hydrolases (GH), carbohydrate esterases (CE), polysaccharide lyases (PL), and carbohydrate-binding modules (CBM)) was performed. This analysis revealed that the A. euteiches genome is enriched in gene coding for proteins with CAZy domains compared to the A. astaci genome (Fig. 4, Additional file 1: ST1j). For example, more than 300 genes coding for proteins with at least one predicted GH domain were detected in the A. euteiches genome, a repertoire size similar to those predicted in oomycete genomes such as Phytophthora species [43], whereas only 109 genes were found in the A. astaci genome.

A closer examination of this repertoire revealed that several enzyme families targeting plant-specific polysaccharides (pectins, hemicellulases) are expanded in A. euteiches and absent in A. astaci. Clearly, the A. astaci genome lacks genes coding for hemicellulases (e.g., GH, 10, 11, CE4 families) and pectinases (e.g., GH28, PL1, PL3, and PL4), which are largely represented in the A. euteiches genome (Fig. 4, Additional file 1: ST1j). Based on the CAZy database, eight GH families correspond to eukaryotic cellulases. The GH5, GH6, and GH7 families, which encode endocellulases and cellobiohydrolases, are highly represented in A. euteiches as compared to other oomycetes. Since oomycete cell walls contain cellulose, these enzymes can play a role in oomycete cell wall remodeling. Interestingly, while GH5 and GH6 cellulases are present in A. astaci and A. euteiches, a large family of GH7 was found only in A. euteiches (39 GH7-encoding gene), suggesting that this class of enzymes could play a role in plant pathogenesis (Fig. 4).

In phytopathogenic fungi it is suggested that α-l-arabinofuranosidases such as GH62 and GH54 are involved in plant penetration and pathogenesis [44, 45]. While GH62 are absent in oomycetes belonging to the Peronosporalean lineage analyzed so far, 12 putative GH62-coding sequences are detected in A. euteiches, while GH54 are missing (Fig. 4, Additional file 1: ST1j). Among GH62, seven are predicted to be secreted. It has also been suggested that α-l-arabinofuranosidases, by degrading arabinofuranose side chains, improve the accessibility of the xylan backbone of hemicellulose to xylanases as GH10 and GH11 [46], which are present in A. euteiches (22 and 3 sequences respectively). Thus, the presence of genes coding for these enzymes suggests that A. euteiches can express a complete repertoire of hemicellulose-degrading enzymes, distinct from those produced by other oomycetes.

The availability of CAZyme repertoires from distantly related oomycetes offers the opportunity to investigate their origin. Since most of these genes are closely related to fungal genes, it has been suggested that some of them were acquired by oomycetes by horizontal gene transfer (HGT) events from true fungi [47]. In certain cases, this hypothesis was sustained by phylogenetic analyses. This is the case for a pectate lyase gene from P. infestans (EEY64154) [48] for which seven orthologous genes, corresponding to the PL1 family, were detected in the A. euteiches genome but not in A. astaci (Additional file 1: ST1j). Phylogenetic analysis using these sequences (Additional file 1: ST1k, Additional file 2: Figure S1) supports the hypothesis that these genes were acquired early in the evolution of oomycetes, before the Saprolegniale-Peronosporale divergence, and probably lost during the specialization of Aphanomyces species to animal hosts.

Numerous CAZymes are associated with CBMs to facilitate enzyme activity [49]. In A. euteiches, many CAZymes are associated with CBMs, and 130 CBM1 (cellulose-binding modules) predicted proteins are detected in this species, but only 7 genes in A. astaci (Additional file 1: ST1j, Fig. 4). In contrast, A. astaci harbors a set of genes encoding proteins with a putative chitin-binding domain (PF03067, IPR004302, CBM33, now reclassified as the AA10 domain in CAZy) that are not detected in A. euteiches (Additional file 1: ST1j, Fig. 4). Interestingly, > 80% (19/23) of these domains are associated with a predicted signal peptide. These data suggested that CBM33 from A. astaci might play a role in the interaction with the crustacean shell, as reported for the PlCBP49 CBM33-containing virulence factor of the honey bee pathogen Paenibacillus larvae that degrades chitin [50].

Aphanomyces spp. cytoplasmic effectors (RxLR and CRN oomycete effectors)

An important discovery regarding oomycete genome sequencing projects resides in the identification of predicted secreted proteins that could be delivered into host cells to aid pathogenicity. Two types of cytoplasmic effectors dominate the oomycete secretome: the RxLR effectors and the Crinklers (CRNs). These effector proteins are characterized by an amino acid signature located at the N-terminus of the sequence. We thus investigated the presence of these host-targeting signals in A. euteiches and A. astaci species. The genomes of both Aphanomyces were searched for RxLRs and CRNs as previously reported [17, 40]. This approach did not allow the identification of RxLR effectors on both Aphanomyces genomes, as previously suggested upon analysis of EST libraries from the same strain of A. euteiches [17]. This analysis sustains the view that RxLR effectors seem to be absent in the Saprolegniale lineage [20]. For predicting CRNs, a combination of automatic searches (hidden Markov model (HMM), regular expression) and manual curation methods was used, and 160 and 31 CRNs and CRN-like genes were detected in the A. euteiches and A. astaci genomes respectively (Additional file 1: ST1k). For A. euteiches, this number is similar to the one reported in phytopathogenic Phytophthora species (> 80 P. capsici; > 200 P. sojae; > 400 P. infestans [32, 51]).

Consistent with previous data on oomycete genomes [7, 52, 53], fewer than 25% of A. euteiches CRNs and none in A. astaci were predicted to be secreted by means of SignalP analysis. By contrast, around 60% of AeCRNs and AeCRN-like genes harbor a predicted nuclear-localization signal (NLS). Among them 16 have both a predicted signal peptide and an NLS. A majority of AeCRNs (> 60%) harbor a LYLAK motif at the N-terminal rather than the canonical Phytophthora LxLFLAK motif, as previously reported upon complementary DNA (cDNA) annotation of A. euteiches [28, 29, 54]. We noticed that the N-terminal trafficking signal is less conserved in A. astaci as compared to CRNs from A. euteiches, and only five sequences harbor a putative header signal in combination with a C-terminal domain (Additional file 1: ST1k). The headers of CRNs are followed by a more diverse C-terminal domain that confers effector activity [32]. Based on sequence similarity, 36 domains were initially defined for the CRN repertoire of P. infestans, and new C-termini domains have been characterized upon the CRN repertoire annotation of various oomycete species [53, 55]. HMM searches and manual assignment on AeCRNs and AeCRN-like sets showed that the necrotic DXZ, D5, and D2 kinase domains are among the most widespread C-termini domains in A. euteiches (Additional file 1: ST1k). The DC domain, reported as a putative DNA-binding helix-hairpin-helix (HhH) motif in PsCRN108 from P. sojae [56], is also largely represented in A. euteiches but not in A. astaci. The DXZ and DN17 domains frequently detected in phytopathogenic oomycetes and the chytrid Batrachochytrium dendobraditis [32, 55, 57] are well represented in A. astaci.

Dynamic changes of A. euteiches transcriptome during infection of M. truncatula roots

To get an insight into the expression and regulation of A. euteiches genes during different life stages, RNA-seq analyses were performed on RNAs isolated from A. euteiches zoospores and mycelium grown in liquid culture as well as in infected M. truncatula roots harvested 1, 3, and 9 days post inoculation (dpi). For the zoospores and the in vitro mycelium libraries, around 50 M reads were obtained, of which 70–77% were mapped in pairs to gene-coding regions (Additional file 3: ST2a). For M. truncatula-infected roots, 46—71 M reads were obtained, of which most could not be mapped to the A. euteiches genome, since a high amount of M. truncatula material was present in the samples. In addition, since less than 1% of the reads of 1 dpi could be mapped to the A. euteiches genome, this time point was excluded from the analysis (data not shown). For time points 3 and 9 dpi, 6–13% of the reads could be mapped in pairs to gene-coding regions (Additional file 3: ST2a). After filtering for lowly expressed genes and data normalization, a multi-dimensional scaling plot revealed a clear separation between expression patterns of the different life stages of A. euteiches, while the biological repeats of each life stage (mycelium, zoospore, 3 dpi, and 9 dpi) clustered together (Additional file 4: Figure S2). This plot suggested that different gene subsets are expressed during M. truncatula infection.

For the 16,786 (85.8%) genes that were expressed, differential expression between the different life stages was determined. Around 36–37% of the genes were differentially expressed between zoospores and the other life stages, while only 10–11% of the genes showed differential expression between the in vitro grown A. euteiches mycelium and the infected root material (Additional file 3: ST2b). This could reflect a transcriptome remodeling during host infection. To further investigate in what processes the differentially expressed genes are involved, a Gene Ontology (GO) enrichment analysis was performed. During infection, genes that are upregulated between 3 and 9 dpi are mainly involved in “carbohydrate metabolic (GO:0005975)” and “proteolysis (GO:0006508)” processes, while downregulated genes belonged mainly to transport processes (GO:0055085, GO:0030001). The opposite situation is detected when comparing zoospores with in vitro grown mycelium (Additional file 3: ST2c–e), suggesting important physiological differences between zoospores and mycelium life stages of A. euteiches as reported for Phytophthora infestans [58].

We then investigated the expression pattern of the A. euteiches secretome with a focus on CAZymes, proteases and protease inhibitors, CBMs, and CRNs/CRN-like effectors (Fig. 5a). When zoospores and mycelium stages are compared, most groups of genes encoding secreted proteins show a similar percentage of differential expression with the exception of the “polysaccharide lyase (PL)” category, which displays the highest percentage of differentially expressed genes. In contrast, during infection (3 and 9 dpi) around three times as much of the genes coding for secreted proteins are upregulated, with proteases, glycoside hydrolases (GH), and PL (86% and 76%) showing the highest percentages of upregulated genes (Fig. 5a). Less than 2% of CRN genes are upregulated at 3 and 9 dpi, while 13% are upregulated in zoospores as compared to in vitro grown mycelium, suggesting that a subset of AeCRN is present at the early stage of Medicago infection (Fig. 5a).

Heatmaps to visualize expression level of genes rather than differential expression show that, for GH, several groups can be identified (Fig. 5b, Additional file 3: ST2f) with groups II, IV, and VI containing genes that are upregulated during infection compared to mycelium. These groups contain a large number of genes coding for proteins with predicted cellulase (GH5, GH6, and GH7) and polygalacturonase (GH28) activity. They also include numerous hemicellulase genes (GH10, GH11) and one secreted arabinofuranosidase (GH62). Most of the polysaccharide lyases (PL) are also strongly induced during infection (Fig. 5d, Additional file 3: ST2g). Secreted PL1 are detected in group II or III, where gene expression is less abundant in mycelium as compared to zoospore and infection stages. This pattern is similar to the one observed in Phytophthora capsici, where the PL1 family in combination with PL16 and PL20 (absent in the Ae genome) account for nearly all of the contribution of the 22 PL genes to Phytophthora virulence [59]. Most of the carbohydrate esterases (CE) showed a consistent low to medium expression in the different life stages (Additional file 5: Figure S3A, Additional file 3: ST2i). For the CBM-containing proteins without a predicted catalytic site mostly composed of CBM1 domains (cellulose-binding), several groups of genes can be identified with group II being strongly induced and expressed during infection while groups III and V are specifically upregulated or downregulated respectively in the zoospore stage (Additional file 5: Figure S3B, Additional file 3: ST2j). This expression pattern suggested that A. euteiches-secreted CBM1 domains may contribute either to A. euteiches virulence or microorganism development, as reported for the CBM1-containing protein CBEL from P. parastica that mediates adhesion to cellulosic substrates and contributes to Phytophthora cell wall architecture [60, 61].

For the secreted peptidases, several groups can also be identified (Fig. 5b, Additional file 3: ST2h). Many peptidases show constitutive expression during all life stages (groups I, VII, and VIII), while some are induced in zoospores (groups V and VI) or during M. truncatula root infection (groups II, III, IV, and IX). These last groups contain a high number of serine proteases (S1 family) that could be considered as candidates for pathogenicity factors. Indeed, serine proteases have been characterized as the major extracellular proteolytic enzymes secreted by Phytophthora spp. and indispensable for necrosis [62]. Notably, none of the secreted protease inhibitors is differentially expressed during infection. However, most of the secreted protease inhibitors are highly expressed during all life stages of Ae (Additional file 5: Fig. 3C, Additional file 3: ST2k).

Together, these analyses indicate that expression of several families of genes coding secreted proteins are induced in zoospores and during infection of M. truncatula roots compared to saprophytically grown mycelium. This supports an important role in the pathogenesis of cell wall-degrading enzymes specifically found in A. euteiches and not in A. astaci targeting plant cell wall polysaccharides (pectins, hemicellulose).

Identification and functional characterization of small secreted proteins in A. euteiches (AeSSPs)

As reported above, 40% of A. euteiches-specific genes that encoded putative secreted proteins did not display a functional annotation and are below 300 amino acid residues in size (< 300 genes, Fig. 3b, Additional file 1: ST1i). Consequently, putative AeSSPs represented 24% of the secretome. Around 5% are organized in clusters (≥ 3 adjacent SSP genes, Additional file 6: ST3a). The clusters are distributed all over the genome and comprise 3 to 9 genes each. The clusters mainly contain AeSSP genes from the same OrthoMCL group, indicating that they might have arisen by duplication.

To ascertain the expression of AeSSPs during host infection, we checked the expression data presented above (e.g., mycelium, zoospore, interaction 3 dpi (T3), interaction 9 dpi (T9)). Since we could not exploit RNA-seq data generated at 1 day post infection, we included data from our previous set of cDNA generated from mycelium (MYC library) and A. euteiches grown in contact with roots during 1 or 2 days (INT T1 + T2, [28] (Additional file 6: ST3a). As shown in Fig. 6a, a large set of AeSSPs are upregulated in zoospores, in mycelium in contact with plant roots, or during infection. In total, 120 AeSSPs have been found to be induced in at least one condition.

We selected in AphanoDB v2.0, the AeSSP cluster from scaffold2_1449549_1653277 (Fig. 6b) and scaffold2_134062_1449527 to evaluate the effector function of AeSSP genes. These clusters are unique, since in addition to their expression at the beginning of the interaction, they contain AeSSPs with a predicted NLS (AeSSP1251 and AeSSP1254, AeSSP1256), suggesting that these proteins could be translocated to the host cell to target nuclear components. To investigate if the AeSSPs can play a role in virulence, Nicotiana benthamiana leaves transiently expressing the AeSSPs were infected with P. capsici (Fig. 7a). In the assay, one side of the leaf was infiltrated with Agrobacterium tumefaciens containing an AeSSP:green fluorescent protein (GFP) fusion construct (AeSSP1250, AeSSP1254, AeSSP1255, and AeSSP1256), while on the other side of the leaf GFP was delivered alone. Quantification of the lesion surfaces at 3 dpi showed that AeSSP1256 significantly enhanced the susceptibility of Nicotiana compared to GFP alone, while no differences were observed for AeSSP1250, AeSSP1254, and AeSSP1255 (Fig. 7b). This indicates that AeSSP1256 may efficiently contribute to oomycete pathogenicity.

To investigate the subcellular localization of AeSSPs, chimeric genes encoding AeSSPs independently with (full-length version) or without (matured form) their own signal peptide under the control of the CaMV35S promoter were generated in fusion with a GFP marker. As shown in Fig. 7c, both versions of AeSSP1250 are detected either in the nucleus or the cytoplasm of Nicotiana cells, while AeSSP1255 versions are mainly cytoplasmic. The NLS-containing AeSSP1254 and AeSSP1256 are nuclear localized without any labeling of the nucleolus. The pattern of fluorescence is similar with the full version or the matured form of AeSSP1254 and AeSSP1256. We noticed that the expression of both versions of AeSSP1256 led to the detection of fluorescence as a ring surrounding the nucleolus in addition to filament-like structures. A closer view of AeSSP1256 in fusion with its own signal peptide confirmed the accumulation around the nucleolus and also as filament-like structures (Fig. 7d). This unexpected nuclear localization is reminiscent of the one observed with the cytoplasmic effector CRN79_188 from Phytophthora capsici upon its transient expression in Nicotiana leaf [63].

Since in all cases the presence of the native signal peptide did not alter the localization of the corresponding AeSSP, we tested whether the AeSSP entered the plant secretory pathway. A functional endoplasmic reticulum (ER) retention signal (KDEL motif) was added at the C-terminal of the full-length form of AeSSP1256 to trap the protein in the ER. In Nicotiana leaves expressing a +SP:AeSSP1256:GFP:KDEL construct, an extensive network throughout the cytoplasm was labeled as observed with an ER-marker [64] (Fig. 8a). In contrast, expression of the matured form –SP:AeSSP1256:GFP:KDEL produced fluorescence in the nucleus. This indicates that +SP:AeSSP1256 is delivered to the ER secretion pathway and that the predicted signal peptide of AeSSP1256 is active in N. benthamiana. To confirm AeSSP1256 trafficking in Nicotiana cells, we used the Brefeldin A (BFA) drug. BFA inhibits transport from the ER to the Golgi and causes the formation of membranous islands throughout the cell [65]. As expected, the localization of the full-length version of AeSSP1256 is disrupted by BFA, while the treatment did not affect the localization of the matured form. Thus, only +SP:AeSSP1256 transits through the Golgi in Nicotiana (Fig. 8b). Finally, we checked the effect and subcellular localization of AeSSP1256 in Medicago composite plants expressing the +SP:AeSSP1256:GFP construct. The expression of AeSSP1256 did not provoke any necrosis or alteration of root system development (data not shown). As illustrated in Fig. 8c, AeSSP1256 in fusion with its own signal peptide is nuclear localized 15 days after root transformation by A. rhizogenes, as observed previously in Nicotiana leaf.

Aphanomyces spp. isolates and DNA preparation

The mycelia of A. euteiches isolate ATCC201684, A. astaci (Genotype E, strain Li07, provided by A. Petrusek, Czech Republic), and A. stellatus isolate CBS 578.67 were grown for 4 days in liquid YG medium (2.5% yeast extract, 5% glucose) at 23 °C in the dark. Biological samples were frozen in liquid nitrogen and the DNA extracted as reported in [31].

Aphanomyces euteiches genome sequencing and assembly

CNS-Genoscope, Evry, France performed the sequencing and assemblies of the A. euteiches ATCC201684 genome, using a combination of Sanger and Illumina technologies. For 454 libraries of A. euteiches, DNA was extracted and fragmented to a range of 5–10 kb or around 20 kb using a HydroShear instrument. Fragments were end-repaired, and the extremities were ligated to 454 circularization adapters. Fragments were size selected respectively to 8 kb or 20 kb through regular gel electrophoresis and circularized using Cre-Lox recombination. Circular DNA was fragmented again by nebulization or using the Covaris E210 instrument (Covaris, Inc., Woburn, MA, USA). Fragments were end-repaired and ligated with library adapters. Mate-pair libraries were amplified and purified. Single-stranded libraries were isolated, then bound to capture beads and amplified in an oil emulsion (emPCR). The libraries were then loaded on a picotiter plate and sequenced using a GS FLX sequencer according to the manufacturer’s protocol. We also prepared 454 single-end read libraries according to the Roche standard procedure using RL (GS FLX Titanium Rapid Library Preparation Kit, Roche Diagnostics, Indianapolis, IN, USA). The libraries were sequenced using a 1/4 Pico Titer Plate on a 454 GS FLX instrument with Titanium chemistry (Roche Diagnostics). For the Illumina library, 2 μg of genomic DNA (gDNA) of A. euteiches was sheared to a 150–700-bp range using the Covaris E210 instrument (Covaris, Inc.). Sheared DNA was used for Illumina library preparation according to a semiautomatized protocol. Briefly, end repair, A tailing, and Illumina-compatible adapter (BiooScientific) ligation were performed using the SPRIworks Library Preparation System and SPRI TE instrument (Beckman Coulter, Simsbury, CT, USA) according to the manufacturer’s protocol. A 300–600-bp size selection was applied in order to recover most of the fragments. The DNA fragments were amplified by 10 cycles of PCR using Kapa HiFi HotStart DNA Polymerase (Life Technologies) and Illumina adapter-specific primers. The libraries were purified with 0.8× AMPure XP beads (Beckman Coulter). After library profile analysis with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and qPCR quantification, the library was sequenced using 100 base-length read chemistry in paired-end mode on the Illumina HiSeq2000 (Illumina, Canton, MA, USA). All reads were assembled with Newbler version vMapAsmResearch-04/19/2010-patch-08/17/2010. The 6319 contigs were linked into 349 scaffolds. The sequence quality of scaffolds from the Newbler assembly was improved as described previously [77] by automatic error corrections with Solexa/Illumina reads (146-fold genome coverage), which have a different bias in error type compared with 454 reads. Finally, the assembly was gap closed using Illumina data and GapCloser [78]. Putative misassemblies were identified and corrected using default parameters of REAPR (version 1.0.17) [79]. Assembly completeness was estimated using BUSCO v3 [35] based on a set of common fungal genes (F) or Alveolata/Stramenopiles genes (AS), aka benchmarking universal single-copy orthologs (BUSCOs). We found 83.1% (F)/78.6% (AS) complete BUSCOs, 3.8% (F)/0% (AS) duplicated BUSCOs, and 3.1% (F)/1.7% (AS) fragmented BUSCOs, leading to 13.8% (F)/19.7% (AS) missing BUSCOs in A. euteiches. Statistics for different sequencing technologies performed in this study for the A. euteiches ATCC201684 genome are presented below.

Aphanomyces astaci and Aphanomyces stellatus genome sequencing and assemblies

The GeT-PlaGe core facility, Toulouse, France realized Illumina sequencing of the A. astaci and A. stellatus genomes, and the assemblies were performed by the Genotoul Bioinformatic Platform, Toulouse, France. DNA-seq libraries were prepared using an Illumina TruSeq DNA v2 Library Prep Kit following the manufacturer’s instructions. DNA was fragmented by sonication, size selection was performed using E-Gel 0.8% (Thermo Fisher Scientific, Waltham, MA, USA), and the adapters were ligated. Ten PCR cycles were applied to amplify the library before final purification with Agencourt AMPure XP beads (Beckman Coulter). Library quality was assessed using an Agilent Bioanalyzer, and the libraries were quantified by qPCR using the Kapa Library Quantification Kit. DNA-seq experiments were performed on an Illumina HiSeq2000 Sequencer using a paired-end read length of 2 × 100 pb with the HiSeq v3 Reagent Kit. The raw reads have been quality checked and stored in NG6 [80]. FastQC (version 0.10.0) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to produce quality metrics and bwa aln (version 0.6.1-r104) to search Escherichia coli-, yeast-, and phage-contaminated reads. The reads of A. astaci and A. stellatus were assembled using MaSuRCa, version 2.0 [81], and the assembly metrics were calculated using the assemblathon_stats.pl script [82].

Structural and functional annotation of Aphanomyces genomes

The assembled data of Aphanomyces were annotated with Augustus v2.757 [83] trained with the assembled RNA-seq transcript generated in this study and the publicly available ESTs previously obtained from A. euteiches [28] using autoAug.pl [84] and PASA [85]. RNA-seq reads were first assembled with trinityrnaseq_r2012–10-05 [86]. The accuracy of the prediction was evaluated by mapping the RNA-seq reads to the genomes using bwa [87]. Genes were annotated using BLASTP against the RefSeq database [88]. For protein family classification, InterProScan [89] and the Pfam protein domain database [90] were used. Gene ontologies were classified based on InterProScan annotation IDs. Kyoto Encyclopedia of Genes and Genomes (KEGG) and Enzyme Commission (EC number) data were obtained with KEGG Automatic Annotation Server (KAAS) predicted protein sequences [91] and were also mapped to the Clusters of Orthologous Groups (COG) classification system [92]. Predictions of signal peptides were performed using SignalP 4.1 [93]. Carbohydrate-Active enZymes (CAZymes) in the protein models were predicted using the dbCAN pipeline [94]. Orthology comparisons between the predicted Aphanomyces genomes and protein datasets from nine deeply sequenced oomycete genomes were performed via OrthoMCLcompanion [95] using standard parameters. Interspersed repeat sequences were de novo identified in the genomes of A. euteiches, A. stellatus, and A. astaci using the RepeatModeler 1.0.8 pipeline [96]. Short (< 1000 bp) Class II transposons that were devoid of recognizable open reading frames were classified as miniature inverted-repeat transposable elements (MITEs). When possible, we assigned them to a superfamily based on the nature of their target site duplication (TSD): piggybac (TTAA motif), Tc1-mariner (TA motif), and PIF-Harbinger (TTA or TAA motif) [97]. The remaining RepeatModeler unclassified consensus sequences (producing no BLASTX hits and having no distinct boundaries) were classified as “unclear”. A library of repeated sequences was then constructed for each species and used to mask each genome with RepeatMasker 4.0.5 [96]. To characterize the global evolutionary dynamics of Aphanomyces TEs, we plotted the frequency distribution of the percent divergence between each consensus and all their cognate copies. To assess the extent to which the TEs are related to each other and to other known TEs, we used each library as a query to perform BLASTN searches on the other libraries and on the Repbase [98] library of TEs (e-value cut-off = 10–20). Gene annotations were visualized in Apollo [99]. All the data generated in this work were incorporated in an updated version (version 2.0) of AphanoDB [37, 100].

Preparation of RNA material

Seeds of M. truncatula Gaertn. F83005.5 were scarified, sterilized, and in vitro cultured as previously described [101, 102]. Roots of 10-day-old plants were inoculated with A. euteiches ATCC201684 zoospores as described in [103]. Infected roots were harvested 1 dpi, 3 dpi, and 9 dpi. In parallel, flasks with 50 mL PG medium (20.0 g/L animal peptone; 5.0 g/L glucose) were inoculated with 1000 zoospores of A. euteiches. Mycelium samples were collected from the PG medium 1 dpi, 3 dpi, and 9 dpi. Zoospores were collected by centrifugation at 14.000× g at 4 °C for 45 min. For all samples three biological replicates were collected. RNA was isolated using the RNA Plant Mini Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s protocol, except for the DNAse treatment (RNase-free DNase, QIAGEN), which was done on the column for 20 min. The RNA concentration was determined, and the quality of the RNA was verified using a Fragment Analyzer (Advanced Analytical, Ankeny, IA, USA).

RNA sequencing

Illumina sequencing of RNA samples generated from A. euteiches mycelium (MY) and mycelium grown in contact with M. truncatula roots (INT) [28], was performed by CNS-Genoscope, Evry, France. Starting with 2 μg of total RNA, double-stranded cDNA was first generated using the TruSeq RNA sample prep kit (Illumina, Canton, MA, USA), and then paired-end libraries were prepared using NEBNext Sample Reagent Set (New England Biolabs, Ipswich, MA, USA). Briefly, the messenger RNAs (mRNAs) were polyA-selected, chemically fragmented, and converted into single-stranded cDNA using random hexamer priming. The second strand was generated to create double-stranded cDNA. The cDNA was then end-repaired and 3’-adenylated, and Illumina adapters were added. The ligation products were purified and the DNA fragments (> 200 pb) were PCR-amplified using Platinum Pfx DNA Polymerase (Life Technologies) and Illumina adapter-specific primers. After library profile analysis by the Agilent 2100 Bioanalyzer (Agilent Technologies) and qPCR quantification (MxPro, Agilent Technologies), the libraries were sequenced on an Illumina HiSeq2000 instrument using 101 base-length read chemistry in a paired-end mode. Statistics for INT and MY libraries are reported below.

In this study new RNA samples were generated from A. euteiches mycelium, zoospores, or infected M. truncatula roots (1 dpi, 3 dpi, and 9 dpi) according to [103, 104]. RNA-seq libraries were prepared according to Illumina’s protocols using the Illumina TruSeq Stranded mRNA Sample Prep Kit to analyze mRNA, at the GeT-PlaGe core facility, Toulouse. Briefly, mRNA was selected using poly-T beads, and cDNA was generated using random hexamer priming, and adapters were then added. Ten cycles of PCR were applied to amplify the libraries. Library quality was assessed using an Agilent Bioanalyzer, and the libraries were quantified by qPCR using the Kapa Library Quantification Kit. RNA-seq experiments have been performed on an Illumina HiSeq2500 using a paired-end read length of 2 × 100 pb with the Illumina TruSeq SBS Sequencing Kits v3. Statistics are presented in the manuscript.

RNA-sequencing data analysis

Reads obtained from the different life stages of A. euteiches (zoospore, mycelium) or infected roots (1 dpi, 3 dpi, and 9 dpi), were mapped to the A. euteiches genome with CLC Genomics Workbench 10.1.1 (QIAGEN). Predefined settings were used except for the similarity fraction, which was set at 0.98, and the read counts per gene (paired reads count as one) were exported. Using the R package edgeR [105], the data were filtered for genes with low counts across all libraries (counts per million (cpm) > 2 in at least three libraries) and then trimmed mean of M (TMM) normalized [106], the tagwise dispersion was estimated, and the differential expression was calculated. Using the generalized linear model likelihood ratio test, p values were obtained, and multiplicity correction was performed by applying the Benjamini-Hochberg method [107, 108]. The multi-dimensional scaling plot was created using the 500 genes with the highest dispersion over all libraries. The enrichment of GO terms in up- or downregulated genes was tested using a Fisher’s exact test with the classical algorithm in the topGO R package [109]. Heatmaps were created using the gplots R package [110], and clustering was performed according to the complete linkage method with Euclidean distance measure.

Phylogenetic and divergence time analysis of Aphanomyces species

To assess the phylogenetic position and divergence times of A. astaci, A. euteiches, and A. stellatus within oomycetes, we followed the methodology of Matari and Blair [38], who produced a robust timetree of oomycetes. We used the amino acid alignments of 40 genes that were assembled for 17 oomycetes and one outgroup species (Tetrahymena thermophila), to which we added the sequences of the three Aphanomyces species. These 40 datasets were selected among 70 initial alignments of genes involved in regulation of gene expression, because there was strong evidence supporting their orthology relationships and they contained minimal missing data [38]. Amino acid sequences were retrieved from the proteomes of the three Aphanomyces species using sequences of the three Saprolegniales (Thraustotheca clavata, Saprolegnia parasitica, and Achlya hypogyna) as queries in BLASTP searches. Most genes could be identified for the three Aphanomyces species, except six in A. astaci (HMG-CBF-BFY, MnmA, p15, Ssl1, TAF6, TAF12) and one in A. stellatus (TFIIB). The 40 alignments were then concatenated and subjected to a maximum likelihood phylogenetic analysis in PhyML v3.0 [111] using the WAG model of amino acid substitutions as in [38]. The robustness of the tree was evaluated by performing 1000 bootstrap replicates. Divergence times were conducted in BEAST v1.8.2 [112] using the same parameters as in [38]. Briefly, we performed 50 million generations run using the PhyML-generated tree as a guide tree and treating each of the 40 datasets as a separate partition with the WAG substitution model, a Yule speciation process (uniform distribution; 0–100; initial value 0.01), and the random local clock model, which was shown to best fit the data by [38]. The calibration strategy follows that of [38] and consists of three calibration points modeled with a gamma prior distribution: (1) diatom node: 5–95% quantiles = 74–100 Myrs, (2) diatom + Ectocarpus node: 5–95% quantiles = 176–202 Myrs, and (3) oomycetes + ochrophytes node: (5–95% quantiles = 418–550 Myrs). The root age was modeled using a uniform prior distribution (408–1750 Myrs; initial value 635). We used Tracer v1.6 [113] to visualize convergence and determine the burn-in and FigTree v1.4 [113] to generate the dated tree of oomycetes.

Construction of plasmid vectors and Agrobacterium-mediated transformations

The AeSSPs sequences were amplified by PCR from A. euteiches gDNA with specific primers (Additional file 6: ST3b). The CACC cloning site was added to each forward primer. PCR products were purified using the PCR Clean-Up Kit (Promega, Madison, WI, USA) and introduced in the pENTRY-D-TOPO vector (pENTR/D-TOPO Cloning Kit, Invitrogen). Positive clones were introduced in a pK7FWG2 vector (Invitrogen) or pAMpAT/YFP vector. After sequencing, positive clones were introduced in A. tumefaciens GV3101 and A. rhizogenes ARqua1. For the KDEL fusion, the GFP construct was amplified from the obtained pK7FWG2 vector by PCR using primers that introduced a C-terminal SEKDEL sequence (Additional file 6: ST3b). Cloning was performed as previously reported using a pK2GW7 vector (Invitrogen). For leaf infiltration, A. tumefaciens GV3101-transformed strains were syringe-infiltrated as described in [61]. For M. truncatula roots transformation, A. rhizogenes ARqua1 strains were used and confocal imaging was performed at 28 dpi, as described in [29].

Brefeldin A treatment and P. capsici infection assay

Brefeldin A (BFA, Sigma-Aldrich, St. Louis, MO, USA) treatment was performed 24 h after agroinfiltration of N. benthamiana. Leaf discs were vacuum infiltrated with TBS1X (50 mM Tris-HCl; 150 mM NaCl pH 7.6) and 15 μg/mL BFA for 30 min, and incubated for 24 h at room temperature. Leaf discs were washed in TBS1X before confocal imaging. For the infection assay, Phytophthora capsici LT3112 was grown on V8 agar plates for 7 days at 22 °C. Zoospore preparation and inoculation on agroinfiltrated N. benthamiana leaf and symptom measurement were performed as reported in [29].

Confocal microscopy

Imaging was performed on Leica DM6B-Z-CS or Leica AOBS TCS SP2 SE laser scanning confocal microscopes. Excitation wavelengths and emission filters were 488 (GFP) or 514 nm (YFP) long-pass. Images were acquired with a 40× water immersion lens or a 25× Fluotar Visir water objective, and corresponded to Z projections of scanned tissues. Image processing was performed using ImageJ software including three-dimensional reconstruction to compute projections of serial confocal sections.