Gene Gain and Loss during Evolution of Obligate Parasitism in the White Rust Pathogen of Arabidopsis thaliana

1For more than 150 years, attempts to culture downy mildews, powdery mildews, and rusts on artificial nutrient media have been unsuccessful. The terms obligate parasitism and obligate biotrophy are used to denote organisms that live in such an obligatory association with living hosts,. Recent research on the obligate biotroph powdery mildew fungus Blumeria graminis or downy mildew oomycete Hyaloperonospora arabidopsidis reveals a close correlation between the biotrophic life style and massive gene losses in primary and secondary metabolism,. Obligate biotrophs form an intimate haustorial interface with parasitized cells. Haustoria are differentiated intercellular hyphae, but little is known about their functionality and evolution beyond their involvement in nutrient uptake,.

2The obligate biotroph oomycete Albugo laibachii is a member of the Chromalveolata, which incorporates Dinophyta, Ciliophora, Heterokonts (containing the oomycetes), Haptophyta, Cryptophyta, and Apicomplexa (which includes human parasites like Plasmodium falciparum and Toxoplasma gondii,).

3Within the oomycetes, A. laibachii belongs to a lineage known as peronosporalean, which includes the hemibiotrophic pathogen of potato Phytophthora infestans and the necrotroph pathogen Pythium ultimum. Within this lineage, obligate biotrophy evolved twice independently in white blister rusts (Albuginales) and downy mildews (part of the Peronosporaceae). The downy mildew pathogen H. arabidopsidis and A. laibachii are both pathogens of the model plant Arabidopsis thaliana. While both show similar infection structures within the host,, A. laibachii releases motile zoospores from asexual spores and sexual oospores, while H. arabidopsidis lacks all motile stages,. Both pathogens are regularly found to co-infect plants and sporulate on the same leaf.

4A remarkable consequence of infection by Albugo sp. is enhanced host plant susceptibility to other parasites to which the host is resistant in the absence of Albugo infection, and also impairment of cell death mechanisms. Albugo sp. infect 63 genera and 241 species, including economically important Brassica rapa (canola), B. juncea (oilseed mustard), and B. oleracea (cabbage family vegetables),. Recent analysis of oomycete evolutionary history suggest that Albugo is more closely related to necrotrophs such as Pythium than to downy mildews, and thus provides a unique system to study the evolution and consequences of biotrophy, and to identify new defence-suppressing effectors and their host targets.

5Since prolonged culture of pathogen strains can result in genetic changes, we sequenced a fresh highly virulent field isolate of A. laibachii. The strain was selected from a heavily infected Ar. thaliana field plot (Norwich, United Kingdom), and strains were single zoospore purified. Isolate Norwich 14 (Nc14) was determined as A. laibachii and used for further analyses. In contrast to Nc14, A. laibachii isolate Em1 (formerly Acem1, A. candida East Malling 1 ) is an established Albugo strain that was collected 15 y ago,,, and we resequenced this strain. Both strains show identical ITS (internal transcribed spacer of ribosomal RNAs) and COX2 (cytochrome C oxidase subunit II) sequences. To ensure that sequence differences observed between these strains are of biological relevance not just the result of background mutations, we tested the host range for both isolates on 126 Ar. thaliana accessions and identified 12 that show resistance to only one of the A. laibachii isolates (Table S1). Nc14 is virulent on more accessions than the Em1 isolate is (Table 1).

6The A. laibachii Nc14 genome was sequenced using Illumina 76-bp paired reads with ∼240-fold coverage (Figure 1). In order to assemble the diploid heterozygous genome, an assembly pipeline was developed using Velvet as primary assembler and Minimus as meta-assembler (Figure S1). Short read assembly programs are sensitive to heterozygous positions depending on read depth and kmer-length. Reads not aligning to bacterial or plant sequence in public databases were used to estimate the genome size as ∼37 Mbp. Using the estimated genome size, 50% of the resulting assembly is contained in 164 contigs with an N50 of 56.5 kbp. A comparative analysis of contig size classes versus frequency indicates that 90% of the assembled genome shows a high degree of continuity in only 585 contigs, while 10% of the genome is fragmented in 3,231 contigs (Figure 2A). Read depth indicates that this 10% of the genome shows elevated levels of nucleotide coverage that are likely to comprise unresolved repeats (Figure 2B). Aligning Illumina cDNA reads from different stages of infection to reveal transcriptionally active regions in the assembly shows that few transcripts arise from the unresolved repetitive regions of the genome (Figure 2D), suggesting that the gene space of a genome can be reliably defined using Illumina-only approaches. A CEGMA analysis revealed a high degree of completeness of assembly of core eukaryotic genes, as well as a continuity within the core genes comparable to high-quality Sanger read assemblies (Figure S2; Table S2). We designed 32 primer pairs for regions between 0.6 and 5 kb based on our assembly (Table S3). Thirty-one genomic regions could be amplified and were Sanger sequenced from both ends. All PCR products had the predicted size, and sequences showed 100% identity to the genome assembly.

7The mitochondrial draft genome was assembled in a separate attempt because of its high repeat content and therefore higher coverage compared to the core genome. The assembled genome comprises 26.7 kb in 11 contigs and shows a high degree of synteny to the P. infestans mitochondrion Ia and the Py. ultimum mitochondrion (Figure S3). Considering the node coverage of the Velvet primary assembly (∼150×), 15.6 kb of the mitochondrial genome have >300× node coverage and seem to be duplicated. This might indicate, comparable to the Py. ultimum mitochondrion genome, that ∼50% of the genome is duplicated, leading to an estimated genome size of ∼43 kb. While the highly repetitive tRNAs are not resolved within the A. laibachii mitochondrial genome, regions of high synteny between the Py. ultimum and the P. infestans mitochondrial genome are found in ribosomal proteins and subunits of the NADH dehydrogenase as well as cytochrome C oxidase.

8Approximately 22% of the A. laibachii Nc14 genome assembly consists of repetitive regions (Figure 3; Tables S4 and S5). The majority of repeats are represented by transposable elements (96%), while 4% of all repeats are A. laibachii-specific (Table S5). Compared to other obligate biotrophs, the number of repeats is low. H. arabidopsidis, for example, with an estimated genome size of 100 Mb, contains ∼43.3% repeats, while transposable elements account for 64% of the ∼120-Mb Bl. graminis (powdery mildew) genome. We identified 45 contigs carrying telomeric repeats; amongst these, 25 contigs have telomeric repeats located at one end of a contig. We therefore postulate that the A. laibachii Nc14 genome is distributed over 12 or 13 chromosomes (Table S6). tRNA genes are difficult to resolve because of their high copy number. Within our Illumina assembly, 153 tRNA genes were detected with 48 distinct anticodons (Figure S4; Table S7). Our ability to resolve all these repeats within the Illumina short read assembly illustrates its quality.

9Based on read depth, both Nc14 and Em1 isolates possess ∼6 Mbp of hemizygous or highly heterozygous regions (6.2 and 5.6 Mbp for Nc14 and Em1, respectively) (Figure 1B and 1D) as well as ∼13,000 heterozygous loci (13,116 and 13,523 for Nc14 and Em1, respectively) (Figure 2C). Remarkably, most of the hemizygous/highly heterozygous regions are shared between Nc14 and Em1.

10Compared to other sequenced oomycetes like P. infestans (240 Mbp), H. arabidopsidis (100 Mbp), or even Py. ultimum (42.8 Mbp), A. laibachii has a highly compact genome structure (Figure 4A). Approximately 50% of the A. laibachii genome assembly matched cDNA reads, and transcriptionally active regions are further clustered, resulting in transcriptional hot spots and silent genomic regions (Figure 4B).

11A reference set of 13,032 gene models was generated incorporating cDNA reads from different stages of infection (Figure S5A). From extensive cDNA sequencing of infected Arabidopsis leaves, approximately 20 M (∼1.5 Gbp) unique Illumina reads match the Nc14 genome assembly but not Ar. thaliana TAIR 9.0, and these were used to generate training sets for ab initio gene predictions and as evidence sets for consensus gene prediction. In all, 88.3% of all gene models are supported by at least three cDNA hits.

12For validation of these gene models, a set of 860 annotated core eukaryotic orthologous groups (KOGs) was compiled and tested. In all, 75% of these groups are present in the current annotation. For comparison, 78% of KOGs were present in P. infestans, 73% in H. arabidopsidis, 42% in Pl. falciparum, and 85% in Ar. thaliana (Figure S5B). In addition, 49.9% of all gene models show Pfam support, resulting in 2,505 Pfam domains, and 803 genes were functionally assigned to pathways using ASGARD and manual annotation. Transcriptional units show an even more compact, clustered occurrence than P. sojae or P. ramorum and an occurrence pattern clearly different from that of P. infestans (Figure 4C).

13From our annotations using ASGARD we identified major enzymes of the lipopolysaccharide biosynthesis pathway, as have been described for P. infestans. These analyses revealed, in addition, the possibility that A. laibachii is able to synthesize brassinosteroids. We identified potential homologues to the Ar. thaliana brassinosteroid biosynthesis genes Dwf4 and DET2 (Table S8). Although ASGARD identified homologues of Br6ox, D2, and CPD, manual annotation revealed that assigning function to members of the superfamily of cytochrome P450 enzymes in A. laibachii is difficult based on homology alone (Table S8). It has been hypothesized that the frequency of functionally redundant genes is reduced in obligate biotrophs, as reported for Bl. graminis. Combining ASGARD and manual annotation we identified the absence of the whole steroid biosynthesis pathway, and, like other oomycetes, A. laibachii probably relies on the host as a source of sterols. We hypothesize that A. laibachii would need to take up campesterol from the plant as a precursor for brassinosteroid synthesis.

14During evolution, plastids of both red algae and green algae were transferred to other lineages by secondary endosymbiosis. How often and when secondary endosymbiosis occurred is difficult to address but of importance to clarify the origin of chromalveolates and their gain and loss of endosymbionts. There are two distinct hypotheses for what took place. The monophyletic hypothesis posits that a red alga was taken up only once, followed by repeated losses of this algal genome, giving rise to the highly divergent group of chromalveolates. An alternative and more common view hypothesizes polyphyletic origins of the Chromalveolata, with in some cases multiple events of secondary endosymbiosis –.

15Molecular divergence of A. laibachii from other species within the Chromalveolata was assessed by examining the percentage of amino acid identity between orthologous gene pairs (Figure 5). These analyses demonstrate that the green alga Chlamydomonas reinhardtii, the brown alga Ectocarpus siliculosus, and the diatom Phaeodactylum tricornutum show the same distribution of percentage amino acid identity to A. laibachii Nc14 regarding the cumulative frequency of orthologous pairs. In contrast, previous systematic analyses suggested that brown algae and diatoms are the closest relatives of oomycetes and that secondary endosymbiosis occurred with a red alga, although there are suggestions that oomycetes diverged before this event. Using a set of >1,700 genes that are of “green” origin (from green algae) or “red” origin (from red algae) and that have been integrated into the diatom nuclear genome, we found more oomycete genes that show significant BLAST hits to green algae than to red algae (34 “green” compared to five “red”) (Figure S6; Table S9). These findings are consistent with the results published by Moustafa et al. for diatoms. In a separate approach we identified genes showing high similarity between oomycetes, green algae, and red algae that are absent from diatoms (32 “green”; 11 “red”) (Tables S10 and S11). This result might indicate the presence of all these genes in a common ancestor, followed by loss or expansion of the gene family depending on adopted live style. To address this question, we further analysed genes absent from A. laibachii Nc14 and studied their presence/absence in three other oomycetes, Pl. falciparum, and the brown alga E. siliculosus (Table S12). The majority of genes absent from A. laibachii Nc14 are absent from other oomycetes and from Pl. falciparum but are present in the brown alga. These genes are involved in the photoautotrophic, aquatic life style of diatoms and algae, such as a sodium/bile acid cotransporter, a haloacid dehalogenase-like hydrolase, fatty acid biosynthesis genes, a zeaxanthin epoxidase and a fucoxanthin chlorophyll a/c binding protein. In contrast to the genes lost, we found that certain gene families like aspartic proteases or proteases containing MORN (membrane occupation and recognition nexus) repeats show expansion in A. laibachii Nc14 compared to in diatoms. Although our results fit the hypothesis of a common ancestor, we cannot exclude horizontal gene transfer and uptake of an endosymbiont after the divergence between a brown algal ancestor and an oomycete ancestor, given the low number of diagnosed genes that we could analyse.

16Potentially green-algae-derived proteins carrying MORN repeat domains (Figure S7) are involved in the complex process of internal budding in apicomplexans, which may be similar to the zoospore formation of oomycetes within oospores or zoosporangia or gamete formation in diatoms. While oomycetes with a motile zoospore stage like A. laibachii and P. infestans carry the MORN repeat proteins, these proteins are absent in the non-motile H. arabidopsidis and absent in the non-motile red alga Cyanidioschyzon merolae. We therefore hypothesize that loss of this gene of hypothetical green algal origin could have led to the evolutionary loss of the whole flagellum apparatus in H. arabidopsidis. However, we cannot rule out that depletion of any major flagellar protein could have caused evolutionary loss of the whole flagellum apparatus. Inspection of the flagellar inner arm dynein 1 heavy chain alpha, which is absolutely necessary for flagellum function, reveals that genomic regions carrying flagellar inner arm dynein 1 heavy chain alpha genes show a high degree of synteny between oomycetes like Py. ultimum and A. laibachii. In contrast, a syntenic region in H. arabidopsidis shows replacement of the flagellar dynein by Mariner- or Gypsy-like transposable elements (Figure S8).

17Since within the peronosporalean lineage, biotrophy evolved twice independently, we compared A. laibachii with the other obligate biotroph H. arabidopsidis, hemibiotroph P. infestans, and necrotroph Py. ultimum (Figure 5; Tables S13 and S14). We found that H. arabidopsidis is the most diverged from A. laibachii. H. arabidopsidis shares the fewest (4,826) orthologous genes with A. laibachii, versus the average of 5,722 in A. laibachii/P. infestans and A. laibachii/Py. ultimum comparisons. Meanwhile, H. arabidopsidis genes show the highest amino acid identity with the genes of P. infestans, on average 73% of amino acid identity between all single copy orthologous pairs.

18Py. ultimum shares the highest number of orthologous genes with A. laibachii (5,910 pairs). P. ultimum proteins also have a slightly higher percentage of amino acid identity with A. laibachii proteins than with other oomycetes (Figure 5). Yet, Py. ultimum itself is closer to H. arabidopsidis and P. infestans than to A. laibachii, sharing with them more orthologous genes with higher mean amino acid identity.

19These analyses support the hypothesis that A. laibachii and H. arabidopsidis evolved biotrophy independently; genes missing in one or the other genome compared to the necrotroph Py. ultimum or hemibiotroph P. infestans may be correlated with biotrophy (Table S15). One of these genes is that for molybdenum-cofactor-dependent nitrate reductase. Nitrate reductase catalyzes pyridine-nucleotide-dependent nitrate reduction for nitrogen acquisition. Both biotroph pathogens have a set of transporters showing homology to amino acid transporters, but other uptake mechanisms or sources could also enable nitrogen acquisition from their hosts. While H. arabidopsidis lost only the nitrate reductase, A. laibachii also lost the sulphite oxidase and the whole molybdopterin (a cofactor required for nitrate reductase and sulphite oxidase function) biosynthesis pathway. In Pl. falciparum, which shows a high degree of adaptation to parasitism, nitrate reductase, sulphite oxidase, and the whole molybdopterin biosynthesis pathway are also missing. Most likely the loss of the two Mo-containing enzymes and the Mo-cofactor biosynthesis is the outcome of biotrophy and not the reason for biotrophy, though conceivably there may have been selection against this pathway if other nitrogen or sulphate sources are less energy-consuming and therefore enhance fitness during parasitism. Molybdenum has been reported to interfere with function of chaperones like Hsp90,. Avoiding the uptake of molybdenum might prevent this Hsp90 inhibition and increase fitness on Ar. thaliana accessions with high molybdenum levels like Col-0. H. arabidopsidis therefore could be in a less advanced stage of host adaptation compared to A. laibachii and Pl. falciparum.

20Besides biotrophy, the formation of haustoria and haustorium-like structures evolved several times in peronosporalean biotroph and hemibiotroph pathogens. Haustoria in fungi are sites of enhanced nutrient uptake and metabolism, such as thiamine biosynthesis. In the oomycetes, all haustorium-forming species have lost the thiamine biosynthetic pathway. We infer that haustorial oomycetes obtain thiamine from the host.

21We therefore hypothesize that evolution to biotrophy is initiated not by gene loss, but rather from the ability to build a haustorium and therefore differentiate a sophisticated interface with a host. The critical step to adopting biotrophy is likely to be efficient defence suppression to enable persistence of functioning haustoria; subsequent loss of biosynthetic pathways is likely to be secondary.

22Well-adapted human pathogens like Pl. falciparum and plant pathogenic fungi like Ustilago maydis have small secretomes (320 and 426 proteins, respectively) compared to necrotrophic fungi like Aspergillus fumigatus (up to 881 proteins ). We found that the same is true for oomycetes. Using SignalP to predict potential secretion signal peptides and MEMSAT to predict transmembrane (TM) domains, we identified 2,473 (2,136 without TM domains) potentially secreted proteins in the hemibiotroph P. infestans and 1,636 (1,222 without TM domains) in the necrotroph Py. ultimum. For H. arabidopsidis only 1,350 (1,054 without TM domains) and for A. laibachii 949 (672 without TM domains) were identified. Analysing the secretome for pathogenicity-related proteins like proteases, glucosyl hydrolases, and potential elicitins or lectins reveals a significant reduction in the H. arabidopsidis and A. laibachii secretome (Tables 2 and S16). We postulate that biotrophs reduce their activation of host defence by reducing their inventory of secreted proteins, particularly cell wall hydrolyzing enzymes.

23The ability to establish a sophisticated zone of interaction like the parasitophorous vacuole in Pl. falciparum or the haustorium in oomycetes and fungi requires sophisticated host defence suppression, which is predominantly achieved via secreted proteins delivered into the host cell,. The A. laibachii secretome comprises 672 secreted proteins without TM domains. Genetically identified oomycete avirulence (Avr) proteins are secreted proteins that have signal peptide and RXLR motifs,. In many oomycete genomes the RXLR motif is over-represented and positionally constrained within the secreted protein. We identified 25 RXLR and 24 RXLQ effector candidates in the A. laibachii secretome. To determine the likelihood that RXLR or RXLQ motifs occur merely by chance in the A. laibachii secretome based on amino acid content, we performed in silico permutation of the motifs (Figure 6A and 6B). We concluded that the RXLR and RXLQ motifs were not likely to occur merely by chance, and that the likelihood of occurrence by chance is higher in the proteome as a whole than among secreted proteins. It was shown for P. infestans that effectors are often located in gene-depleted repetitive regions of the genome. We therefore investigated RXLR candidate proteins in highly repetitive regions of the genome. We identified two RXLRs, one in a highly conserved repeat region with ∼10 repeats in Nc14 and one in a more diverged repeat region with >80 repeats within the genome. The first region also exists in A. laibachii isolate Em1; the diverged repeat of the second identified region exists but without the RXLR gene-containing region (Figure S9). There are 563 RXLR effector candidates identified in P. infestans, so RXLR effectors are less likely to be relevant for A. laibachii virulence.

24Similar conclusions can be drawn for the CRN protein family, which shows expansion in P. infestans, but not A. laibachii, where only three members of the CRN family could be identified with signal peptides. Eight additional CRN-like proteins were identified where no signal peptide has been predicted.

25To identify new classes of effectors in the Albuginales clade, the secretome of A. laibachii was computationally screened for genes either showing heterozygosity or showing nucleotide polymorphisms between Nc14 and Em1. We identified a new class carrying a “CHXC” motif by inspection of the first 80 amino acids after the signal peptide cleavage site. CHXC candidates are significantly enriched within the secretome (Figure 6C). Comparisons of the N-terminal part of the CHXC proteins revealed additional conserved amino acids, particularly a glycin at +6 to the CHXC motif (Figure 6D).

26In host–pathogen interactions, intraspecies comparisons enable the search for virulence alleles that undergo positive selection and fixation within the population,. Secreted proteins with close contact to the host cell, such as effector proteins, often show enhanced levels of positive selection,. By comparing the two A. laibachii isolates Nc14 and Em1, we identified a significantly higher frequency of non-synonymous to synonymous mutations within the predicted secretome compared to the rest of the proteome. Our analyses showed that this was particularly true for heterozygous positions and less convincing for homozygous SNPs (Table S17). Genes that are highly conserved between species, like KOGs, showed comparable non-synonymous and synonymous substitution rates, with a slight excess of synonymous mutations. There are significantly more genes within the KOGs showing a non-synonymous/synonymous ratio less than 1 than genes with values greater than 1. Comparing this to candidate effector classes like RXLRs, RXLQs, and CHXCs reveals that in particular the CHXCs show significantly higher frequencies of non-synonymous to synonymous mutations. This supports the idea that the CHXC sub-class of secreted proteins is under positive selection, similar to other described oomycete effectors like ATR1 or ATR13 from H. arabidopsidis,.

27Further to this we identified Nc14 genes absent or highly diverged from the Em1 complement. We defined a gene as absent or highly diverged if >10 bp showed 0 coverage in the Em1 alignment. Out of the 672 secreted proteins without TM domains, we identified seven as absent from Em1 (1.04%). We also detected two with a predicted TM domain (0.73%) that are absent from Em1. Regarding all gene models, 96 were absent (0.74%). This finding is a further indication for a greater selection pressure on secreted than on non-secreted proteins, as has been found in species or interspecies comparisons in Phytophthora sp. and Ustilago/Sporisorium.

28We tested A. laibachii effector candidates (one CHXC, one RXLR, and one CRN effector candidate) for their host delivery efficiency using a P. capsici–Nicotiana benthamiana translocation assay. Briefly, N-terminal domains of candidate effectors were fused to the P. infestans Avr3a effector domain, transformed into P. capsici, and tested for whether they confer translocation of Avr3a into N. benthamiana carrying R3a, resulting in avirulence. Statistical analyses of the delivery efficiency (Figure 7) clearly indicate that the A. laibachii CRN3 N-terminus and CHXC9 N-terminus are as efficient as the Avr3a N-terminus in Avr3a translocation, while the RXLR1 N-terminal domain is less efficient. An alanine replacement construct of the CHXC motif supports the importance of this motif for delivery efficiency. The Avr3a C-terminus alone confers a low basal delivery level without the need for the N-terminal enhancer. These findings reveal the potential of the CHXC proteins to be delivered into the host cell, similar to RXLRs and CRNs, though the delivery mechanism for all these effector classes requires further investigation.

29To assay the effectors for virulence function, we used Pseudomonas syringae pv. tomato (Pst) DC3000 luciferase carrying “effector detector vector” (EDV) constructs to deliver effectors into the plant cytoplasm via type III secretion (Figure 8). Tests on Ar. thaliana Nd-0 plants revealed that several selected A. laibachii RXLRs, CRNs, and CHXCs enhance virulence compared to a non-functional AvrRps4 (AvrRps4[AAAA]). On Ar. thaliana Col-0, in contrast, the CRN and one RXLR (RXLR1) do not enhance virulence while RXLR2 and CHXCs still do. These tests indicate that CHXCs carry the capacity to enhance virulence in phytopathogenic bacteria, perhaps by suppression of host resistance mechanisms,. These virulence assays together suggest that A. laibachii uses at least three different major effector classes.

30To try to identify the evolutionary source of CHXCs, we investigated enrichment of CHXC-motif-containing proteins in the secretomes of P. infestans, Py. ultimum, H. arabidopsidis, Saprolegnia parasitica, Thalassiosira pseudonana (diatom), Pl. falciparum (Apicomplexa), E. siliculosus (brown alga), C. merolae (red alga), Ch. reinhardtii (green alga), Volvox carteri (green alga), and Ar. thaliana. Only A. laibachii contained a significant enrichment of CHXCs in its secretome. Although not significantly enriched, both the fish pathogen S. parasitica and the land plant Ar. thaliana contained more than ten CHXC proteins carrying potential secretion signals (14 and 11, respectively) (Figure S10). In contrast to CHXC-containing proteins, almost all inspected organisms show a high number of CXHC-containing potentially secreted proteins; a common CXHC protein is protein disulphide isomerase (Table S18).

31Given that A. laibachii CHXCs show the closest clustering with S. parasitica, V. carteri, Ch. reinhardtii, and Ar. thaliana CHXCs (Figure 9), conceivably this candidate effector class evolved from an ancestral green-alga-derived gene.

32Whatever their origin, we conclude that CHXC proteins are present in all organisms analysed but evolved effector function only in Albuginales and possibly Saprolegniales. In Albuginales, one N-terminal sub-class of CHXCs (CHxCLx(4)Gx(5–6)L) shows significant expansion, with 23 members, while other CHXCs are distinct from this clade. S. parasitica CHXCs are distinct from this major A. laibachii clade and therefore remain to be tested in future experiments.

33The A. laibachii genome assembly sheds light on the evolution of biotrophy since it allows the first comparison, to our knowledge, of two oomycete obligate biotroph pathogens (A. laibachii and H. arabidopsidis) that evolved biotrophy independently. In addition, A. laibachii shows the highest overall amino acid identity to the necrotroph pathogen Py. ultimum and the hemibiotroph P. infestans. One of the striking results of this comparison is that all organisms able to build haustoria have lost their thiamine biosynthesis pathway, presumably because thiamine is easily obtained from hosts via the haustorial interface. A closer interface requires effective host defence suppression. We therefore hypothesize that the evolution of biotrophy involves a series of steps: step 1, involving progressively more effective effectors to suppress defence, step 2, attenuated activation of defence by reduction in the inventory of cell wall hydrolyzing enzymes, resulting in, step 3, weak selection to maintain certain biosynthetic pathways if the products of the pathways can be directly obtained from the host. This results in progressively more comprehensive auxotrophy and culminates in irreversible biotrophy (Figure 10).

34An infected leaf was harvested from an Ar. thaliana plant grown in a heavy infected field plot in Norwich (UK; 52.6236,1.2182) in December 2007. Zoosporangia were washed off the leaf surface and used to infect Ar. thaliana Ws-0-eds1 plants. After 1 wk one pustule was punched out, and spores were placed on ice for 30 min to release zoospores. Unhatched zoosporangia were removed by filtration, and zoospores were diluted to ∼10 zoospores/ml and sprayed on Ar. thaliana Ws-0 plants (∼100 µl/plant). This procedure was repeated 4× until spores were bulked up on Ar. thaliana Ws-0 plants. Zoosporangia were harvested using a home-made cyclone spore collector.

35Zoospores were suspended in water (105 spores/ml) and incubated on ice for 30 min. The spore suspension was then sprayed on plants using a spray gun (∼700 µl/plant), and plants were incubated in a cold room in the dark over night. Infected plants were kept under 10-h light and 14-h dark cycles with a 20°C day and 16°C night temperature.

36High molecular weight DNA was extracted from zoosporangia using a phenol/chloroform-based purification method after grinding in liquid nitrogen, adapted from. Library preparation for Illumina sequencing was performed as described. All data were generated using paired-end reads. 800 bp and 400 bp paired-end sequencing libraries were constructed, and 8.8 Gbp of usable data were generated (for read and insert length, see Figure 1A).

37Figure 1A lists all reads after purification from plant and bacterial contamination as well as all reads aligned to the assembly. In summary, 91.6% of all reads can be aligned to the contigs, suggesting 2.8 Mbp missing from the assembly. Since 32.7 Mbp are in the assembly, the genome can be estimated to 35.5 Mbp. In another approach considering all reads and their read length, 8.8 Gbp (∼7% correction for lower quality of second read pair) were generated, which would lead to an expected coverage of the 32.7 Mbp genome of ∼270×. The mean coverage using single copy genes (glycolysis and TCA) is 240×. Considering the 2.5 Mbp of repeats (Figure 1B, right side, coverage underestimated) with an average coverage of 1,086×, which is ∼4.4 times more than the mean coverage of the contigs, this repeat region corresponds to 10.9 Mbp. In contrast to this, the genome contains ∼6.2 Mbp of hemizygous regions (Figure 1B, left side, coverage overestimated). These calculations suggest a genome size of ∼43 Mbp, given all repeats resolved, or an effective genome size of ∼37 Mbp.

38A. laibachii–infected Ar. thaliana Ws-0 plants were harvested 0 (after cold room, see plant inoculation), 2, 4, 6, 8 and 10 d after infection. Total RNA was extracted using TRI Reagent RNA Isolation Reagent (Sigma), and Dynabeads (Invitrogen) were used to enrich for mRNA. First and second strand cDNA synthesis was performed according to manufacturer's instructions using the SMART cDNA Library Construction Kit (Clontech), and cDNA was normalized using the Trimmer kit from Evrogen. cDNA samples were mixed in equal amounts and fragmented using a Covaris sonicator (Covaris). Illumina libraries were prepared as described for fragmented genomic DNA.

39Data for comparative genomics were downloaded from the sources listed in Table 3.

40First Velvet was used, running different kmer-lengths and different sequencing library subsets (kmer-length: 23, 31, 41, 45, 49, 55, 61, 67, and 73; subsets: 400-bp insert only, 800-bp insert only). N50 number and length were determined for each of the assemblies, and the best assembly was selected as the matrix to be used with the Minimus2 genome merge pipeline. For the current assembly the 400-bp only subset with kmer-length 61 was used as matrix, and for kmer-lengths 49, 55, 61, 67 and 73, all 400- and 800-bp assemblies were added (Minimus parameters: consensus error <0.001; minimum identity >99%; 20-bp maximum trimming). A set of genes showing high heterozygosity was used to ensure that contigs were properly joined. Parameters were changed through several rounds, and minimum overlap, in particular, was lowered from 100 bp to 15 bp. An overlap of 15 bp was found to be the optimum for difficult heterozygous regions. After each Minimus assembly, all reads were back aligned to the contigs using MAQ aligner. Regions showing less than 3× average coverage were removed, and redundant fragments were removed using BLASTN with an e-value cut-off of 1e−20 and 99.9% identity. After this step a next round of Minimus was started, with changing minimum overlap in steps of 20 bp down from 100 bp. Below 20 bp steps were changed by 5 bp (See Figure S1 for work flow).

41Since it is impossible to cultivate obligate biotrophs under sterile conditions, plant and bacterial contaminations were removed by using BLAST against genome sequences of the host plant Ar. thaliana (TAIR 9.0), fungal genomes (Neurospora crassa), oomycetes (H. arabidopsidis), and diverse bacterial genomes (Xanthomonas sp. and Pseudomonas sp.).

42To identify heterozygous loci, Illumina reads were aligned using MAQ, and the SNP detection pipeline was used according to the manual, with default parameters and minimum coverage greater than 180× for the Nc14 alignment and greater than 20× for the Em1 alignment. From the MAQ SNP file, positions were selected where two bases are possible and maximum coverage was less than 350×.

43Assembled repetitive elements were identified using the RepeatScout program (http://bix.ucsd.edu/repeatscout/) with a seed size of 14. The frequency of elements and their location in the assembly were estimated with RepeatMasker using a library of repetitive elements built up by RepeatScout. A sequence was considered to be repetitive if it occurred in the genome assembly on at least three different contigs. The resulting library was searched for the sequences homologous to the known transposon elements using TBLASTX (e-value cut-off of 1e−5) and a database of transposons, RepBase. Consensus repeats that matched predicted Nc14 protein coding genes were filtered out. The remaining consensus repeats that do not match any sequences deposited in the NCBI database or any known transposon element and that do not overlap with Nc14 protein coding genes represent either Albugo-specific repeats or simple repeats.

44tRNA genes were predicted with the program ARAGORN using first default parameters and second options allowing introns in the gene sequences.

45CEGMA was used according to the manual with a local installation.

46For the combined ABySS and Oases assembly, adaptor sequences from the SMART kit cDNA synthesis were removed for the ABySS assembly, and the ABySS program was used according to the manual. Different kmer-lengths were tested, and a length of 61 used for the final assembly.

47Untrimmed cDNA sequences were assembled using Velvet and a kmer-length of 51, 57, 61, and 71. Oases was used for the final assembly of the contigs according to the manual, using default parameters.

48MUMmer in maxmatch mode was used to combine all ABySS and Velvet assemblies. Redundant contigs were removed using BLAST.

49Since the assembled cDNA is not strand specific but orientation is needed for gene prediction, cDNA 5′ tags were generated by Illumina sequencing (E. Kemen, A. Balmuth, J. D. Jones, unpublished data). Using Bowtie aligner, cDNA 5′ tags were aligned onto the assembled cDNA and, based on tag counts, orientated in the 5′ to 3′ direction.

50To map assembled cDNA against the genome, either BLAT in trimT and fine mode or PASA with default settings was used.

51Illumina reads were directly mapped to the genome using the Bowtie aligner, in “best” mode and with strand correction (strandfix mode). Pileup files were generated using bowtie-maqconvert and maq pileup allowing four mismatches per 76-bp read. To incorporate this data as hints files for gene prediction, regions with greater than 3× coverage were extracted.

52To generate a reliable gene set to train further programs, GeneMark was used for ab initio gene prediction. ORFs plus 50 bp on the 3′ end and 50 bp on the 5′ end were extracted, and Illumina-sequenced cDNA was aligned to the ORFs using Bowtie. Gene models were selected if the coverage within the ORF didn't drop below three. This dataset with more than 2,000 genes was used as “traingenes” for the automated training program provided with the Augustus package (autoAug.pl). The trained Augustus program was then used for gene prediction including the combined Oases/ABySS-assembled cDNA (mapped using BLAT) as evidence. Default parameters (extrinsic.ME.cfg) were used for all predictions.

53For consensus gene predictions with P. infestans, SGP2 was used according to the manual.

54ASGARD alignments were converted into GFF files to be used for consensus predictions.

55Consensus gene models were generated using Evigan. cDNA from assemblies and alignments was converted into GFF files and combined with Augustus, GeneMark, SGP2, and ASGARD predictions. The genome was than screened for gene-free regions, and Augustus gene predictions were added if available. In a third round, regions that did not contain consensus gene models or Augustus gene models were extracted, and GeneMark annotations were added if available.

56A set of genes was further tested by 5′ and 3′ RACE to validate start and stop sites.

57Molecular divergence of A. laibachii from other species was assessed by examining the percentage of amino acid identity between orthologous gene pairs.

58Orthologous pairs were identified using the OrthoMCL program with an e-value cut-off of 1e−5. Alignments of protein pairs were performed with MUSCLE.

59Amino acid identity was calculated only for the single copy genes by either excluding alignment gaps from calculations or taking gaps into account. The results show similar trends, so we present only results for the calculations when alignment gaps were excluded.

60The total number of orthologous groups identified between species and the number of one-to-one orthologous pairs, as well as a mean amino acid identity, are shown in Table S7. In the comparison of T. gondii and A. laibachii, we found few orthologous pairs represented by the single copy genes (23 pairs); therefore, we excluded this pair of species from the analyses of sequence divergence.

61We also estimated the levels of amino acid identity for the core eukaryotic genes (orthologous genes shared by all examined species); these data are presented in Table S8.

62To identify A. laibachii genes with sequence similarity to green- or red-algal-derived diatom genes, a set published by Moustafa et al. was used. All A. laibachii proteins showing homology to genes identified by Moustafa et al. were further blasted (BLASTP) against the Ch. reinhardtii gene set, the E. siliculosus gene set, the U. maydis gene set, and the Fusarium oxysporum gene set with an e-value cut-off of 1e−20. Genes were considered to be green-alga-derived only if the protein was absent from U. maydis and F. oxysporum but present in Ch. reinhardtii, and was considered red-alga-derived if not in U. maydis or F. oxysporum but in E. siliculosus. The same analyses were performed on the Saccharomyces cerevisiae, Pl. falciparum, H. arabidopsidis, P. infestans, Py. ultimum, V. carteri, Ch. reinhardtii, C. merolae, C. merolae, Th. pseudonana, and Ph. tricornutum gene sets.

63A. laibachii candidate genes with significant sequence similarity to green or red algae and other oomycetes (e-value cut-off of 1e−20) but not to fungi, brown algae, or diatoms were identified using the criteria in Table 4. Representative organisms for each group are as follows: green algae: V. carteri, Ch. reinhardtii; red algae: C. merolae, Galdieria sulphuraria; fungi: F. oxysporum; brown algae: E. siliculosus; diatoms: Ph. tricornutum, Th. pseudonana; oomycetes: P. sojae, Py. ultimum, H. arabidopsidis.

64Homologues between oomycetes, fungi, brown algae, and diatoms were identified using OrthoMCL (e-value cut-off of 1e−20 or 1e−5).

65Synteny between multiple species was analysed using the Artemis Comparison Tool. Alignments between genomic sequences were performed using TBLASTX with a score cut-off of 210. Annotations of P. infestans, Py. ultimum, and H. arabidopsidis were transferred using TBLASTN with an e-value cut-off of 1e−30. LTR_FINDER was used to annotate long terminal repeats (LTRs) within the genomic sequences, and coordinates were manually added. Regions between LTRs were blasted against RepBase to identify the presence and/or type of transposon.

66Secreted proteins were predicted using a local installation of SignalP 3.0. Proteins were considered to be secreted if both the neural networks and hidden Markov model methods predicted the protein to have a signal peptide. Predictions of TM domains were performed after removing the predicted secretion signal. TM domains were identified using MEMSAT3. Proteins were considered to be without a TM domain with p non-TM>0.0004 or, for high stringency, p non-TM>0.01.

67To identify new motifs, subsets of secreted proteins were selected and analysed using MEME with default parameters. Identified motifs were tested against the whole gene set and the Swiss-Prot database using MOTIF Search. In a second step, motifs were selected only if they were positioned within 50 amino acids after the secretion signal.

68Tests for over-representation of an identified motif were done using motif and sequence shuffling. Secreted proteins were predicted as described in the previous section, and the signal peptide was removed prior to further analyses.

69Each of the sequences without secretion signal was randomly shuffled 30 times. After each shuffling the sequences were screened for the motif in question. If the motif was identified after shuffling, the sequence was excluded from the next round. If the motif was never identified within the 30 times shuffling, the motif in the original protein was counted as “unique empirical”. All possible combinations of the amino acid sequence within the motif were calculated. For each of these permutations, the “unique empirical” proteins were calculated.

70The 30 times shuffling was repeated 1,000 times to calculate background levels. Background levels were defined as how often a sequence was found again having the motif or the permutated motif. This was called “background (mean)”. Motifs that were above this background were considered for further analyses.

71The second criterion was if a motif was significantly enriched in the secretome compared to all non-secreted proteins. For statistical validations we calculated the cumulative hypergeometric probability.

72Candidates for further experiments were evaluated according to a ranking list. Maximum possible score was nine points, and the following scores were given: one point for being on a shorter, repetitive contig (≤3,000 bp) or end of contig, since we assumed that effector candidates might be in repetitive regions as shown for P. infestans effectors; one point for having cDNA support; two points for being a short protein (≤400 amino acids); two points for carrying one of the identified motifs (RXLR, RXLQ, CHXC, CRN); one point for being expressed before day 10 after infection; one point for being expressed before day 4 after infection; and one point for showing SNPs in the Em1 comparison.

73N. benthamiana plant genotypes and P. capsici strain LT1534 were grown and cultured as described by Schornack et al.. P. capsici transformation was performed as described by Schornack et al..

74Phytophthora transformation constructs SP_AVR3aC, RXLR1_AVR3aC, CRN3_AVR3aC, CHXC9_AVR3aC, CHXC9AAAA_AVR3aC, CHXC7_AVR3aC, and CHXC7AAAA_AVR3aC were synthesized and cloned into pTOR by Genscript. Fusion genes were flanked by ClaI (5′) and SacII (3′), and internal AscI sites were inserted between the N-terminal effector domain and AVR3aC coding domain. N-terminal domains used are listed in Table 5.

75Phytophthora infection assays were performed according to Schornack et al. with slight modifications.

76Vector pTOR::Avr3a and pTOR::Avr3a (AAAA-AAA) were obtained from Dr. Steve Whisson. The control construct SP_Avr3aP was synthesized using the signal peptide of P. infestans Avr3a, fusing the signal peptide directly to the Avr3a C-terminus (GenBank accession number ACX46530.1). Translocation fusion constructs were synthesized using N-terminal coding sequences of Albugo RXLR1 (Gene name: AlNc14C278G10072, GI: 325190660), CRN3 (Gene name: AlNc14C196G8578, GI: 325188975), CHXC9 (Gene name: AlNc14C832G12555, GI: 325193652), and CHXC7 (Gene name: AlNc14C191G8449, GI: 325188831), and the AVR3a C-terminus (GenBank accession number ACX46530.1, GI: 260594559). CHXC9_AAAA_AVR3aC and CHXC7_AAAA_AVR3aC fusion constructs were synthesized by replacing the CHXC motif with a quadA motif. All plasmid suspensions used for P. capsici transformation were prepared using the Qiagen Midi Prep kit (Qiagen). For a summary of constructs see Table 5.

77Candidate RXLR effectors were cloned from RXLR to stop; all other candidate effectors were cloned from SP cleavage site to stop into pENTR D-TOPO (Invitrogen) and mobilized into pEDV6. The resulting effector∶pEDV6 constructs were conjugated into Pst DC3000 luxCDABE and Pst DC3000 ΔAvrPto/ΔAvrPtoB. The contribution of an individual effector was assessed by spray inoculating 4- to 5-wk-old short day grown plants as previously described.

78Growth of Pst DC3000 luxCDABE effector∶pEDV6 was calculated by measuring whole plant luminescence using a Photek camera system and normalizing this to plant fresh weight.

79To assess the virulence of Pst DC3000 ΔAvrPto/ΔAvrPtoB effector∶pEDV6, bacterial colony counts were performed as previously described.

80All Illumina sequence reads generated during this study have been submitted to the Sequence Read Archive at EBI and are accessible under the accession number ERA015557. Individual studies are available with accession numbers ERP000440 (Alias: albugo_laibachii_nc14_dna_sequencing, http://www.ebi.ac.uk/ena/data/view/ERP000440), ERP000441 (Alias: albugo_laibachii_nc14_cdna_sequencing, http://www.ebi.ac.uk/ena/data/view/ERP000441), and ERP000442 (Alias: albugo_laibachii_em1_dna_resequencing, http://www.ebi.ac.uk/ena/data/view/ERP000442).

81All contigs and annotations are available through EBI or NCBI. The accession range is from FR824046 to FR827861 (3,816 contigs including annotations) and can be accessed through the ENA browser (http://www.ebi.ac.uk/ena/).