Comparative genomics of downy mildews reveals potential adaptations to biotrophy

Read processing and assembly

Approximately 230 million 100 base-par (bp) paired-end reads were generated for both isolates that were reduced by 32% (R13) and 18% (R14) after quality and k-mer trimming. Read filtering identified approximately 58 and 97 million read pairs associated with other oomycete assemblies and 65 and 77 million un-associated read pairs for R13 and R14, respectively. Optimal k-mer lengths for highest output statistics were empirically found to be between 91 and 99 nt for both isolates. The top five BLASTn [26] filtered assemblies of R13 were assembled using k-mer sizes 92 to 96 and had scaffold N₅₀s ranging from 20.9 k-base-pair (kb) to 23.8 kb with between 12 and 16 scaffolds over 100 kb with BUSCO [27] completeness (protist library) ranging from 94.9 to 97.9%. The top five R14 assemblies (k-mer sizes 92–95 + 97) had scaffold N₅₀ ranges of 18,157 to 18,478 bp, with 18 to 23 scaffolds over 100 kb and BUSCO completeness ranging from 95.3 to 97.4%.

Post mitochondrial fragment filtering, pairwise merging of assemblies, from lowest scaffold count size to highest, and removal of fragments under 1 kb, the scaffold N₅₀ of isolates R13 and R14 rose to 52.8 kb and 42.8 kb respectively and the highest BUSCO [27] score previously observed was retained in both isolates, although the duplication rate had increased. A single round of redundancy removal with Redundans [28] and Haplomerger2 [29] resulted in the final assemblies, with little redundancy detected by BUSCO. Complete statistics for the intermediate assemblies generated are provided in Additional file 2.

P. effusa R14 was assembled in to 1275 contigs and 880 scaffolds. The contig N₅₀ was 51.71 kb and the scaffold N₅₀ was 61.4 kb. The assembly totaled 30.8 Mb, contained 0.56% gaps and had 47 scaffolds over 100 kb, one of which was over 250 kb. JELLYFISH analysis produced an estimated haploid genome size of 41.2 Mb with 27.9 Mb inferred as single copy (Additional file 3). BUSCO reported a completeness of 97.0%, with no duplicated predictions and 0.4% fragmented predictions. The remaining 2.6% BUSCOs were reported as missing (Table 1).

Assembly quality

SyMap [31] plots showed that the two P. effusa assemblies had a high degree of collinearity although they are highly fragmented (Fig. 1). There were 237 and 214 scaffolds over 50 kb for R13 and R14, respectively; inter-isolate alignments were identified for 141 R13 and 134 R14 of these scaffolds. When compared to P. sojae, 177 R13 and 167 R14 scaffolds over 50 kb could be aligned with 23 P. sojae scaffolds over 50 kb.

Repeat masking and LTR analysis

Repeat libraries for both isolates were generated using RepeatModeler [32], identifying 98 elements in R13 and 88 in R14. In both isolates most elements identified were Long Terminal Repeat (LTR) Retrotransposons (Additional file 4). The LTR profiles of the two isolates of P. effusa were similar to one another, although the density profile indicates that isolate R14 may have more similar Copia elements than R13 (Additional file 5).

Annotation

Initial SNAP [33] hidden Markov model (HMM) training with MAKER [34] resulted in 7853 gene predictions for R13. Subsequent runs with intermediate SNAP bootstrapping produced predictions ranging from 6493 to 9601 models, the mean residue length ranged from 407 to 482 amino acids. Equivalent bootstrapping of the same HMM run on P. effusa R14 produced predictions ranging from 6691 to 9034 gene models, the mean residue length ranged from 407 to 476 amino acids (Additional file 6). Further assessment of each individual unfiltered run showed mean BLAST [26] scores back to the training data ranging between 569 to 671 for R13 and 603 to 678 for R14 (Additional file 6). When orthology analysis was performed for each run against the training data, the number of orthogroups detected varied from 4803 to 7117 for R13 and 5363 to 7073 for R14 (Additional file 6). The mean e-value for Pfam [35] domains, scoring under 1e^− 5, of individual MAKER runs ranged from 1.99 e^− 7 to 2.20 e^− 7 for R13 and 1.99 e^− 7 to 2.23 e^− 7 for R14 (Additional file 6). For R13 the MAKER run with the highest mean BLAST score, second highest orthogroup count, and lowest mean Pfam e-value was considered the best. For R14, the best annotation set scored the highest mean BLAST score, highest orthogroup count, and lowest mean Pfam e-value (Additional file 6). Investigating gene models at unique loci from alternative MAKER runs did not produce a high scoring set of gene models; therefore, integrating unique models from multiple runs was not performed.

Comparative analysis

Orthology analysis between the two isolates of P. effusa identified 7314 overlapping orthogroups, of which 6833 were single copy. Nine orthogroups, containing 54 genes and a further 653 singletons unassigned to orthogroups were found exclusive to R13, while 10 orthogroups, containing 70 genes and a further 638 singletons were exclusive to R14. These isolate-specific genes represented 8% of the total predicted genes (Fig. 2a). Reads were mapped between isolates to determine the inter isolate coverage of coding regions and to determine if SNPs/short-indels could be detected in genes to determine why the proteins were inferred as absent. Mapping reads of R14 to R13 revealed 40 R13 genes that had low coverage, one of which encoded an RxLR effector. Seven genes contained indels within their coding regions. Additionally, 28 genes, including one encoding a necrosis inducing protein (NPP1), had SNPs which introduced a premature stop codon. SNPs in 14 genes resulted in the loss of their start codon and the loss of the stop codon in 22 genes, relative to the genes predicted in R13. When R13 reads were mapped to R14, ten genes had low coverage, including two RxLR effectors, which had a single residue difference between one another. Fifteen genes, including one RxLR contained indels within their coding regions. Additionally, 31 genes, including one CRN, had SNPs which introduced premature stop codons, 22 genes, including one RxLR and a gene encoding a protein kinase domain, had SNPs which caused the loss of the stop codon and nine lost their start codon, relative to the genes predicted for R14. These proteins, especially those with putative virulence/necrotic function may be useful as diagnostic markers for each race (Additional file 11). Repeating the analysis on the two previously sequenced isolates of P. tabacina [23] revealed 6818 overlapping orthogroups, 3869 of which were single copy and 19% of the gene models reported as unique to either of the isolates (Fig. 2b).

There were 709 orthogroups containing at least one protein encoding a Pfam domain inferred as under-represented in Peronospora species/H. arabidopsidis and grouped under one of the previously defined categories (Table 4). These were visualized to investigate the absence of the orthologs in downy mildew species (Fig. 3). While a core component of each orthogroup, typically the high gene number groups in all categories except phytopathogenicity, are retained across all oomycetes tested (Table 6), the majority of orthogroups for each category have missing orthologs for Peronospora sp. and H. arabidopsidis. When these are grouped with P. halstedii over 50% of orthogroups have detectable orthologs for downy mildew species in every category except phytopathogenicity, with a much larger fraction of motor and calcium associated domain encoding orthologs being identified than carbohydrate and transporter orthologs in P. halstedii (Table 6). The counts and IDs of the proteins contained within each orthogroup are supplied as Additional file 13.

	All		>  1 Phytophthora species		>  1 downy mildew		>  1 Peronospora lineage species
	# Orthogroups	# Single Copy Orthogroups	# Orthogroups	# Single Copy Orthogroups	# Orthogroups	# Single Copy Orthogroups	# Orthogroups	# Single Copy Orthogroups
Calcium	174	127	174 (100%)	127 (100%)	133 (76.4%)	92 (72.4%)	62 (35.6%)	29 (22.8%)
Carbohydrate	41	22	41 (100%)	22 (100%)	23 (56.1%)	7 (31.8%)	20 (48.9%)	4 (18.2%)
Flagella/Motor	77	50	77 (100%)	55 (100%)	57 (74.0%)	30 (60.0%)	27 (35.1%)	8 (16.0%)
Phytopathogenicity	165	41	163 (98.8%)	40 (98.6%)	75 (45.5%)	21 (51.2%)	61 (37.0%)	15 (36.6%)
Transporter	252	127	251 (99.6%)	126 (99.2)	144 (57.1%)	59 (46.5%)	117 (46.4%)	43 (33.9%)

K-mer analysis and heterozygosity

KAT density plots [36] were made to investigate heterozygosity in P. effusa. Two clusters of 21-mers are visible in both plots in addition to the many low-frequency 21-mers due to sequencing errors and contaminants (Fig. 4a). R14 has strong homozygous k-mer signal and a weak heterozygous signal at half coverage. R13 has strong homozygous k-mer signal and a significant signal at higher than half coverage. The same analysis of two P. tabacina isolates [23] detected strong homozygous k-mer signal and two heterozygous signals at half and quarter coverage consistent with the presence of multiple distinct haplotypes due to polyploidy, heterokaryosis, or mixtures for both isolates (Fig. 4a).

Spectra-cn plots that are also generated by KAT [36] were used to investigate the frequencies of k-mers in common between the read sets and assemblies of P. effusa. R14 contains one significant peak with the majority of 21-mers being represented once in the assembly, consistent with a high-quality assembly of a predominantly homozygous organism (Fig. 4b). The majority of 21-mers of the R13 read set were represented once, regardless of the k-mer frequency in the read set; this indicates that neither cluster of k-mers in R13 was heterozygous because most of the k-mers in the lower coverage peak were incorporated in the assembly rather than the anticipated proportion (50%) of the k-mers being absent. Therefore, the assembly of R13 seems to be of high quality but these k-mer plots are not consistent with a simple diploid genome. In contrast, spectra-cn plots of P. tabacina confirmed that both isolates contained three clusters of k-mers with a high proportion of k-mers from the first peaks absent in the assembly, approximately half the k-mers absent from the second peaks, and a small fraction absent from the homozygous peaks. Both isolates of P. tabacina had high levels of k-mer duplication; a significant proportion in the homozygous peak was represented twice as much as expected in the assembly indicating that both P. tabacina assemblies are under-assembled (Fig. 4b).

Reads were mapped back to the respective assemblies to identify single nucleotide polymorphisms (SNPs) and investigate the frequency of reads supporting the alternative allele in P. effusa. In R13, 106,714 heterozygous SNP sites were identified and 74,690 in R14, indicating 0.33 and 0.24% heterozygosity for R13 and R14, respectively (Fig. 4c). Plots of the frequency of reads supporting the alternative allele at each SNP revealed a clear peak at 0.5 in both isolates as expected in a diploid; however, a smaller second peak was present for both isolates at ~ 0.33. When only SNPs in genes were considered, the SNP count was reduced to 31,041 and 17,943 inferring a 0.23 and 0.13% heterozygosity in the predicted gene space of R13 and R14, respectively. Plots of the frequency of reads supporting the alternative allele of genic SNPs retained the peak at 0.5, although the peak at 0.33 was greatly reduced, containing ~ 5.7% (R13) to 8.6% (R14) of the genic SNPs (Fig. 4c). Collectively, these data indicate that the allele frequencies of numerous SNPs in P. effusa are not consistent with the 1:1 ratio expected in a diploid organism.

The normalized coverage of each predicted gene in R13 was more variable than in R14 (Fig. 4d). While the majority of the ~ 8600 predicted genes had around 1 to 1.2x normalized coverage in both isolates, 851 gene models in R13 had 0.6 to 0.8x coverage; in contrast, only 67 genes in R14 had 0.6 to 0.8x coverage. There was no significant deviation in representation of the 513 Pfam domains encoded by these genes.

Phylogenetics

A maximum likelihood tree with 1000 bootstraps was produced from 49 concatenated, single copy genes predicted by BUSCO [27] from both P. effusa isolates and all published downy mildew assemblies available from NCBI, plus three diverse Phytophthora spp. and rooted with Pythium ultimum as the out-group (Fig. 5). Two downy mildew clades were evident. P. effusa was in the larger clade with P. tabacina and clusters with Pseudoperonospora cubensis, H. arabidopsidis and Sclerospora graminicola. The second downy mildew clade is made up of the three isolates of two Plasmopara spp., which were more closely related to P. infestans than to the other downy mildews.

Mitochondria

The mitochondrial genomes of race 13 and 14 isolates were both 41,318 bp in size with a GC content of 22.8% (GenBank accessions MH142315 and MH325167, respectively). The coverage of the mitochondrial assemblies was 2081 and 1003x for R13 and R14, respectively. Sequences of both genomes were identical and had the same organization as P. tabacina (KT893455) with the exception of an inverted repeat (IR) that was present in P. effusa. Coding regions constituted 93.9% of the genome with 13.3% of this total representing hypothetical coding regions. A total of 35 known genes (encoding 18 respiratory chain proteins, 16 ribosomal proteins, and an import protein, ymf16 of the secY-independent pathway), the rnl and rns, and 25 tRNA genes encoding for 19 amino acids were present. In addition, there were five hypothetical proteins (ymf98, ymf99, ymf100, ymf101 and orf32) in common with other oomycete mitochondrial genomes [23, 37–42] and five putative ORFs that were unique to P. effusa. Four of the unique ORFs (orf 181, orf131, orf201 and orf277) were located between the rns and cox2 genes (in the same location as the unique putative ORFs of P. tabacina but encoded in the opposite orientation) and another between the atp1 and nad5 genes (orf209). BLAST queries to GenBank identified no significant sequence similarities for orf181, orf201, and orf209; however, there was moderate sequence similarity of orf131 and orf 277 with putative ORFs in mitochondrial genomes of P. tabacina and Phytophthora sojae (DQ832717). In orf131 of P. effusa, bases 287 to 393 were 75% identical to the 5′ end of orf269a of P. tabacina, while bases 7 to 262 were 69% identical to the 3′ end of this putative ORF. The 3′ end of P. tabacina orf269a (bases 323 to 743) shared similarity with putative ORFs from P. sojae with 85% identity to the spacer and 5′ end of orf116, while bases 59 to 320 are 89% identical to part of orf101. Bases 8–677 of P. effusa orf277 were 73% identical to all but the last 100 bp of P. tabacina orf269b and bases 8–785 sharing 70% identity with the 5′ end of P. tabacina orf290.

A feature of the P. effusa mitochondrial genome that was not present in P. tabacina was the presence of an IR. The first arm of the IR was located between orf201 and orf277 with 382 bp of the 5′ end representing the 5′ end of orf201. The second arm of the IR was located between atp1 and orf209 with 382 bp of the 5′ end of the IR the 5′ end of orf209. The sequence for both arms of the IR is not a perfect match because there is a 2 bp deletion after base 667 in the first arm relative to the second repeat, hence the sizes are 872 bp and 874 bp, respectively. In addition to the IR there was also a 47 bp repeat present in two head to tail copies between orf277 and cox2.

Phenotyping of isolates and DNA extraction

Genomic DNA samples for sequencing were obtained from two isolates of P. effusa collected from grower fields in Monterey County, California, in 2012 and 2013. The pathotypes of these isolates were determined by inoculation onto a differential set of spinach cultivars as previously described [8]. For both isolates, leaves of plants of a single cultivar showing heavy sporulation were collected and the spores were scraped off the leaf surface into water in a 50 ml tube as well as vortexed to remove additional sporangia. The suspension of the sporangia was transferred to a microfuge tube and spun at 21000x g for 1.5 min. The resulting pellet was washed in 1 ml of 95% ethanol, spun at 21000x g for 3 min and the pellet frozen at − 80 °C. Four hundred microliters of the Macherey-Nagel NucleoSpin Plant II kit (Düren, Germany) buffer PL1 and a single microfuge tube cap full of glass beads (Sigma G8772) was added and vortexed. The suspension was heat shocked at 65 °C with 10 μl RNase A solution, followed by another high-speed vortex. A 100 μl volume of chloroform was added, followed by a brief vortex and centrifugation for 5 min at 21000 g. The supernatant was added to a NucleoSpin® column (Machery-Nagel NucleoSpin Plant II kit), and the manufacturer’s plant DNA extraction protocol was followed.

Flow cytometry

Flow cytometry was performed on sporulating and pre-sporulating spinach leaves mixed with 1 cm² of young leaf tissue from Oryza sativa cv. Kitaake (2C = 867 Mb), which was sufficiently different from the genome size of P. effusa for use as the internal reference. The O. sativa 2C DNA content was determined by calibrating against nuclei from flower buds of Arabidopsis thaliana Col-0 which has a known absolute DNA content of 2C = 314 Mb [62]. Nuclei extraction and staining with propidium iodide was done using the Cystain PI absolute P kit (Sysmex, Lincolnshire, IL). Flow cytometry was done on a BD FACScan (Becton Dickinson, East Rutherford, NJ). For each measurement, 10,000 nuclei were assessed. Data was analyzed using FlowJo (Ashland, OR).

Sequencing

Illumina TruSeq DNA libraries were prepared and sequenced at the Center for Genome Research & Biocomputing, Oregon State University (Corvallis, OR; http://cgrb.oregonstate.edu/core). DNA was quantified using a Qubit HS dsDNA assay (Invitrogen, Carlsbad, CA) and sheared by sonication followed by end repair, adenylation of 3′ ends, and adapter ligation. Fragments were purified by excision from an agarose gel, enriched by PCR, and the library was quantified with a Qubit HS dsDNA assay (Invitrogen). Sizing of the library was done using Agilent Bioanalyzer HS-DNA chip (Agilent Technologies, Waldbronn, Germany), with final quantification by qPCR using a KAPA Library quantification kit. The median library fragment sizes were 516 bp and 365 bp, for R13 and R14, respectively. Paired end libraries were prepared from genomic DNA of two P. effusa isolates (R13 and R14) and sequenced, 100 bp paired-end on an Illumina HiSeq 2000.

Assembly

Reads were adapter and quality trimmed using BBMAP [63] and mapped to a reference containing bacterial and oomycete genomes available from NCBI, using BWA MEM, v0.7.12 [64], with flags -aC. Paired reads which mapped to oomycete genomes or failed to map to any organism were then advanced to assembly with MaSuRCA v2.3.2 [65]. Assemblies were done, specifying a JELLYFISH [30] size of 1 × 10¹⁰, in iterative k-mer steps of 10, ranging from 31 to 91, and 99 with additional assemblies performed at single step k-mer sizes flanking the highest scoring assemblies (measured on N₅₀, assembly size, number of scaffolds > 100 kb, total number of scaffolds, BUSCO [27] score). The top five assemblies of each isolate were then positively filtered for oomycete scaffolds against NCBI nt, and negatively filtered for mitochondrial associated scaffolds with the mitochondrial assembly (produced as described below) using BLASTn [26]. Assemblies made up of scaffolds with a top BLASTn hit against oomycete scaffolds and with a minimum scaffold size of 1 kb, were then merged in a step-wise manner from lowest scaffold number to highest using Quickmerge [66]. Repeat libraries were generated with RepeatModeler [32], assemblies soft-masked with RepeatMasker [67] and secondary haplotypes collapsed, first with Redundans [28], then Haplomerger2 [29]. Final assembly statistics were obtained using BBMAP [63] and compared to previously assembled genomes. Completeness statistics were obtained with BUSCO.

Assemblies were aligned against one another with NUCmer [68] (−l 100), and both were independently aligned against P. sojae v3.0 with PROmer [68] (−l 30). These were ran and visualized as part of SyMap [31], with a minimum dot requirement of 5 and allowing merged blocks. REAPR [69] summary scores of the final and intermediate assemblies were calculated as previous [70]. REAPR summary scores of both isolates increased as post-assembly processing was performed on each isolate with the final draft assemblies presented here having the highest score (Additional files 2 and 14).

Annotation

The annotation workflow is depicted in Additional file 8. Ab initio annotation was performed with SNAP [33] as part of the MAKER pipeline [34]. Initial predictions were generated with MAKER from all oomycete EST and protein data available on NCBI (est2genome = 1, protein2genome = 1) masking with the above generated repeat library. These predictions were used to produce a HMM using SNAP default settings (fathom, forge, hmm-assembler.pl) and then bootstrapped through subsequent rounds of MAKER, the first using the est2genome and protein2genome evidence for prediction, with subsequent runs turning this off to use only the SNAP HMM for gene prediction (est2genome = 0, protein2genome = 0). All maker runs surveyed single exon proteins with a minimum nucleotide length of 240 being considered (single_exon = 1, single_length = 240). The optimal run was identified through comparative analysis of alternative predictions, namely runs were scored for BLAST [26] hits to the Oomycete training protein sequences, % orthology detected with the Oomycete training database as detected by OrthoFinder [71] and average e-value of Pfam [35] domains detected with InterProScan [72] with a value under 1e^− 5.

Additional putative effector gene models were identified from all ORFs predicted from the genome over 80 residues in length with no missing data. All ORFs were surveyed for secretion signals using SignalP v4.1 (−u 0.34) and independently with SignalP v3.0 (default settings) [73, 74]. ORFs were considered secreted if a positive result was obtained through either or both approaches and no trans-membrane domains were detected by SignalP v4.1. Further filtering, removing ORFs targeted to the mitochondria was performed using TargetP [75]. All ORFs predicted (i.e. regardless of signal peptide prediction) were surveyed for Crinklers (CRN) and WY repeats using hidden Markov models (HMMs) with HMMER v3.1b1 [76]. The CRN HMM was produced with hmmsearch, inputting an alignment of all labelled Phytophthora spp. CRNs from NCBI. The WY HMM was described previously [77, 78]. Predicted secreted ORFs were surveyed for RxLR motifs and previously described/queried degenerate [GHQ]xLR or RxL[GKQ] motifs [53, 79–82] linked with either degenerate [DE][DE][ER] motif or a WY repeat [78]. Predicted effectors were integrated into the Maker produced GFF file by identifying over-lapping gene models on the same strand with BEDTools [83]. Gene models were characterized as to whether they shared either or both start and stop codon positions.

Candidate effectors at unique genomic loci were further curated to ensure that they did not overlap or neighbor one another within 1 kb on the same strand. In instances when not the case, these models were manually refined, with the aid of BLAST [26], to produce multi-exonic RxLR-(EER)-WY models [84, 85].

The final predicted protein sequence was then generated from the GFF using GAG and putative effectors were re-identified on the entire annotation set to obtain final putative effector counts, which include those independently predicted by MAKER [34]. Conserved functional domains, including pathogenicity associated domains, of the final gene models were identified with InterProScan [72] and putative gene model names were assigned through stringent BLAST to the UniProt reference database [86]. These features were added to the NCBI table file with Annie and GAG [87, 88].

Comparative genomics

De novo identification of genes encoding pathogenicity associated domains for P. effusa, other downy mildew and Phytophthora species was performed by running InterProScan on all available gene models. Gene models were inferred to contain a pathogenicity associated domain if such a domain was identified by CDD, Gene3d, PANTHER, Pfam, PRINTS, ProDom, ProSitePatterns, ProSiteProfiles, SMART or SUPERFAMILY [35, 72, 89–94]. The total frequency of gene models encoding pathogenicity associated domains was calculated as the number of unique gene models identified with any of the specified domains divided by the total number of gene models surveyed. T-tests were performed on the frequency of models containing each domain and on the total frequency of pathogenicity associated gene models, segregating the results as downy mildew vs. Phytophthora and Peronospora spp., and H. arabidopsidis vs. Phytophthora spp. and P. halstedii. Chi-squared testing was performed on Pfam [35] domains detected across all seven species. In this test, a protein was counted once for each unique domain it encoded (i.e. it was not weighted if it encoded multiple of a single domain type). Expected scores were weighted for the number of proteins detected as having Pfam domains in each species. Tests were performed on all domains as P. effusa vs. all, Peronospora spp. vs. all, Peronospora spp. and H. arabidopsidis vs. all, and downy mildews vs. Phytophthora spp. A Bonferroni-adjusted p-value was calculated based on the number of domains tested, with Pfam domain scoring below this p-value manually investigated. Instances where an over-representation of domains in P. tabacina appeared to score a significant result were investigated, owing to the high duplication in the P. tabacina assemblies detected in this study.

Orthogroups were inferred using OrthoFinder [71]. First gene models of P. effusa isolates were compared in a two-way analysis. Proteins detected as absent between isolates were investigated by mapping reads of one isolate to the other isolates assembly using BWA MEM v0.7.12 [64] and calling SNPs/Indels with SAMtools mpileup v0.1.18 [95]. For comparison, the same was replicated on the two P. tabacina isolates. These isolates were then combined into a set for their respective species and a seven-way comparison was performed, including H. arabidopsidis, P. infestans, P. ramorum, P. sojae and P. halstedii. Intersects were visualized with upsetR [96] and orthology coefficients (OC) were calculated as OC = (C/T₁) x (C/T₂), where C is the number of overlapping orthologous groups, T₁ is the total number of orthologous groups identified in sample 1 and T₂ is the total number of orthologous groups identified in sample 2. Orthogroups which contained gene models encoding Pfam domains, inferred as under-represented by the Chi-squared analysis, in downy mildew assemblies were identified and analyzed for presence/absence of each species. An orthogroup was considered single copy if only 0–1 models were present for each Phytophthora spp., P. halstedii and H. arabidopsidis and 0–2 models detected for P. tabacina and P. effusa, owing to a) two isolates being used in the analysis for both Peronospora spp. and b) the high level of duplication detected in P. tabacina assemblies. Finally, orthogroups containing previously identified P. infestans nitrogen and sulphur assimilation enzymes [45] were identified to obtain orthologs of these enzymes in downy mildew species.

K-mer and read analysis

JELLYFISH [30] 21-mer hashes for individual read files were generated and histograms were obtained to estimate the genome size of both isolates. These were plotted with R [97] to obtain the 21-mer boundaries of single copy regions of the genome. K-mer based genome size estimates were calculated by summing the results of k-mer density multiplied by its frequency. Estimates of the size of the single copy portion of the genome were produced by limiting this calculations to k-mers with densities between the limits of the k-mer peak profiles (R13; 82–340, R14; 142–560; Additional file 13). Hashes were visually inspected through KAT density plots [36], for both P. effusa isolates and comparatively for two previously sequenced P. tabacina isolates [23]. Hashes generated from pairs of read files were also compared to assembly hashes generated from resulting isolate assemblies and visually inspected through KAT spectra-cn plots.

Heterozygosity of each isolate was estimated by calling SNPs from reads mapped back to either assembly, using SAMtools v0.1.18 mpileup [95] and calculating the number of heterozygous sites, defined as those with an allele frequency between 0.2 and 0.8. Plotting the frequency of the alternative allele of heterozygous SNPs was performed by extracting the number of reads supporting the reference and alternative allele at each bi-allelic SNP site and obtaining the frequency of reads supporting the alternative allele and binned to the nearest hundredth. For SNP sites to be counted they had to be covered by a minimum of 50 reads. The number of SNPs per bin were summed and plotted in R [97]. Variant call files containing the SNP coordinates were intersected with genic lines of the annotated gff using BEDTools [83] to obtain SNP counts within genic regions and the allele frequency plots were produced as above. Read coverage per gene was calculated by obtaining the number of reads per genic locus, using BEDTools multicov [83] requiring the assembly file and an indexed binary alignment map file, multiplied by the read length, divided by the length of the gene. The coverage was normalized for each isolate (calculated from a BEDTools genomecov plot) and plotted in R using ggplot2 [98] applying a multiplicate bandwidth adjustment of ten.

Phylogenomics

Phylogenetics of BUSCO [27] orthologs across a panel of published, publicly available downy mildew assemblies was carried out with select, high quality Phytophthora assemblies and rooted with Pythium ultimum. BUSCO was run independently on every assembly to be surveyed and single copy orthologs from each were identified. Amino acid sequences of orthologs present in all isolates were then aligned independently using MAFFT v7.123b [99] (Additional file 15). Alignments were concatenated and RAxML v8.0.26 [100] was run with 1000 bootstraps and the PROTGAMMAAUTO substitution model. The resulting tree was rooted with P. ultimum and visualized in Geneious and labels were manually placed to improve legibility.

Mitochondrial assembly and annotation

Contigs from a de novo genomic assembly in CLC Genomics Workbench (v8; Qiagen, Redwood City, CA) were identified as mitochondrial due to sequence similarity with P. tabacina mitochondrial sequences (KT893455) by BLAST analysis. These were used as templates for further assembly with SeqMan NGen (v4.1.2, DNASTAR, Madison, WI, USA). The resulting assemblies were evaluated for uniformity and depth of coverage. Contigs were broken when gaps/low coverage or inconsistencies were observed and the set of smaller contigs reassembled using the small templated assembly option of SeqMan NGen to extend the ends of the contigs and the close gaps. ORFs were predicted and annotated with DS Gene v1.5 (Accelrys, San Diego, CA) using the universal genetic code with confirmation of gene identities using BLAST [26] analysis against mitochondrial genome sequences published for P. tabacina, Pythium, and Phytophthora spp. [23, 38, 40] tRNA coding regions were placed using tRNAscan-SE v1.3.1 [101].