Genomes and virulence difference between two physiological races of Phytophthora nicotianae

Background

Phytophthora nicotianae, also known as Phytophthora parasitica var. nicotianae, is a soil-borne bi-flagellated oomycete plant pathogen, which causes black shank in cultivated tobacco (Nicotiana tabacum), and root rot, leaf necrosis, and stem lesions in a variety of plants [1]. P. nicotianae is able to infect a wide range of hosts, spanning 255 genera in 90 different plant families. It devastates the production of a number of economically important plants, and causes millions of dollars worth of economic losses each year in the tobacco industry alone [2, 3]. So far, management strategies for P. nicotianae are limited to non-host crop rotation, cultivation of pathogen-resistant breeds, and the use of chemical control (e.g. mefenoxam) [4]. The primary reason for the difficulty in controlling P. nicotianae is the production and survival of chlamydospores in unfavorable growth conditions, as well as the production of motile zoospores. The ability of P. nicotianae to infect specific tobacco cultivars with different resistance genes defines four physiological races (0, 1, 2 and 3). The predominant physiological races, 0 and 1, are widely distributed throughout China, the United States and other major tobacco cultivation areas [5, 6]. Previous studies using tobacco cultivars with moderate or high levels of resistance have found that race 0 has better pathogenic and ecologic fitness levels than race 1, suggesting that the difference in virulence between the two races is affected by additional genetic factors [7]. To discover better and more efficient ways to control the pathogen, we undertook a global examination of the genes involved in the infection process from different races. Although five strains of P. parasitica are already public available [8], none of these includes any of the four physiological races of P. nicotianae. Here we report the genomes of P. nicotianae physiological races 0 and 1, sequenced using a combination of PacBio single-molecule real-time (SMRT) sequencing technology, and Illumina HiSeq and MiSeq sequencing technologies, and identify candidate genes that may cause the difference in virulence between them.

Isolation of P. nicotianae races and genomic DNA extraction

Tobacco plants infected by either P. nicotianae race 0 or race 1 were obtained from Yunnan Tobacco Research Institute. Any surface dirt on the infected plant was washed off under tap water. After drying, stem tissue from the lesion margin were cut into 5 × 5 mm squares, sterilized using 70 % ethanol for 1 minute, and then rinsed three times using sterile water. Sterilized tissue squares were then placed in lima bean agar (LBA) plates amended with 50 μg/ml ampicillin, 100 μg/ml rifampicin, and 50 μg/ml of pentachloronitrobenzene to suppress possible contaminant. LBA plates were incubated for 2–3 days in darkness at 25 °C. Color and texture of the colony and mycelium were used to confirm the identity of P. nicotianae. Mycelium was transferred to LBA slants and cultured for 7 days in darkness at 25 °C. Genomic DNA was extracted using the modified cetyltrimethyl ammonium bromide method [9].

Sequencing and quality control

Assembly

Because of its relatively high heterozygosity, we used a hybrid assembly approach to assemble the genome of race 0 (Fig. 1). We also compared the performance of different assemblers including Velvet (v1.2.09) [12], ABySS (v3.81) [13], JR-Assembler (v1.0.3) [14], EULER-SR (v1.1.2) [15], SPAdes (v3.0.0) [16], SOAPdenovo2 (r240) [17], Celera Assembler (v8.0) [18] and Minimus2 (v 3.1) [19] on P. nicotianae race 0. Comparison showed that assemblies from PacBio reads were generally of better quality than those from HiSeq reads (Fig. 2). The final assembled genome sizes for race 0 and race 1 were 80 Mb and 69 Mb, respectively, which is slightly different from the previous estimation of P. nicotianae (90 Mb) [20]. The corresponding contig N50 sizes were 23 kB and 30 kB, respectively (Table 2). Over 95 % of core eukaryotic genes could be mapped to the two genomes using CEGMA [21], and over 90 % of Illumina HiSeq reads could be mapped back to the genome assemblies.

Races	Categories	Number	N50 (kb)	Longest size (kb)	Size (Mb)	Percentage of the assembly
P. nicotianae race 0	Total contigs	6142	23	203	81	-
	Total scaffolds	6139	23	203	81	-
	Genes annotated	17,797	1.8	26	27	33 %
	Transposable elements	-	-	-	32	39 %
P. nicotianae race 1	Total contigs	6116	22	196	71	-
	Total scaffolds	5317	30	318	71	-
	Genes annotated	14,542	1.7	15	20	28 %
	Transposable elements	-	-	-	27	38 %

Annotation

Known transposable elements (TEs) were identified with RepeatMasker (version 3.2.6) [22] using the Repbase TE library (v16.10) [23] and default parameters. Tandem repeats were predicted using TRF [24]. gypsy and copia types of long terminal repeat (LTR) were the main contributors to the repeat, making up 12.5 % and 3.5 % of the genome for race 0, and 11.5 % and 3.6 of the genome for race 1. For gene structure prediction, gene sets from 9 species including Phytophthora infestans [25], Phytophthora sojae [26], Phytophthora ramorum [26], Hyaloperonospora arabidopsis [27], Pythium aphanidermatum [28], Pythium arrhenomanes [28], Pythium irregulare [28], Pythium vexans [28], Pythium iwayamai [28] and Pythium ultimum [29] were used for homology-based prediction. GENSCAN [30], AUGUSTUS [31] and GlimmerHMM [32] were used for de novo gene prediction. Evidence derived from homology-based and de novo predictions were then integrated in GLEAN to generate a consensus gene set. A total of 17,797 and 14,542 protein-coding genes were annotated in P. nicotianae race 0 and race 1, respectively. Over 97 % of these genes could be aligned against KEGG [33], Swiss-Prot and TrEMBL databases [34]. Mean exon numbers per gene in P. nicotianae and related species varied between 2.2 and 2.8, suggesting that homology and de novo-based prediction were appropriate for annotation (Additional file 1). We also used publicly available expressed sequence tags (ESTs) from the appressorium [35] and mycelium [36, 37] of P. nicotianae to validate the annotation. We retrieved a total of 10,524 ESTs from the dbEST database. Using the threshold of match length >200 bp and E-value <1e-5, we aligned 8,043 ESTs to the race 0 genome and 7,618 ESTs to the race 1 genome. Additionally, 4,454 genes in race 0 and 3,604 genes in race 1 were supported by at least one EST (Additional file 2). Whole genome comparison using NUCmer [38] found that average identity was 99 % for 1-to-1 alignment, and 98.84 % for m-to-m alignment between P. nicotianae races 0 and 1. Using KaKs_Calculator, mean synonymous mutation ratio (Ks) was estimated to be 0.075 between race 0 and race 1 [39], and four genes were identified to be positively selected (Additional file 3).

Gene family clustering and evolution

Gene family clustering using OrthoMCL [40] revealed that over 72 % of gene families were shared between species pairs among P. nicotianae race 0, race 1, and related species (Additional file 4). The average number of genes per gene family was 1.19 to 1.50 in Phytophthora and 1.14 to 1.26 in Pythium, suggesting more copies of genes exist in the Phytophthora genus. A total of 1,604 single-copy genes were identified between P. nicotianae and the other 9 related species (Additional file 5). Gene family expansion and contraction estimated using CAFÉ [41] found that 1,237 gene families expanded and 294 gene families contracted in race 0, while 217 gene families expanded and 508 gene families contracted in race 1 (Fig. 3).

ABC transporter expanded in P. nicotianae

The ATP-binding cassette transporter (ABC transporter) superfamily facilitates the transport of ions, proteins, lipids and toxins across plant membranes [42]. Interestingly, a domain-centric study found this gene family to be enriched in the oomycete plant pathogen genomes [43]. It was proposed that an important function of ABC transporters in pathogens involves exporting toxic phytoalexins [44, 45]. Based on the result of CAFÉ analysis, we found the ABC transporter gene family to be significantly expanded in the branch of P. nicotianae (likelihood ratio test, p-value < 0.05), but not in the branch of P. infestans (likelihood ratio test, p = 0.9). To verify this result, we used Pfam to annotate ABC transporter domains (PF00005.22, PF00664.18, PF01061.19) between P. infestans, and P. nicotianae races 0 and 1 (Additional file 6). The portions of ABC transporters in P. nicotianae were significantly larger than that in P. infestans (chi-square test, p < 0.05). This result suggests that the ABC transporter family plays important roles in P. nicotianae in its adaptive evolution to the host.

Distribution of effectors and their differences in races 0 and 1

Plant pathogens have evolved to secrete effectors, which can manipulate the host immune system and suppress host defense. Based on their target sites in the host plant, effectors can be classified into two classes: (1) apoplastic effectors, which are secreted into plant extracellular spaces; and (2) cytoplasmic effectors, which are translocated into the plant cell. Some effector genes, e.g. ATR5 in H. arabidopsidis, are found to be avirulence genes [46]. These genes are under selective pressure to evade host recognition while maintaining their original functions.

RxLR effectors are important cytoplasmic effectors that contain a conserved N-terminal motif (Arg-X-Leu-Arg). The RxLR motif is involved in translocation into host cells [47]. During infection, the RxLR family functions to suppress host immunity. This process usually involves manipulating plant immunity-associated signaling pathways. For example, PexRD2 can perturb MAPKKKε signaling pathways to suppress NB-LRR-mediated immunity in P. infestans [48]. In addition, a set of RxLR effectors from P. infestans can suppress the signaling pathway induced by flg22, a kind of microbe-associated molecular pattern (MAMP) [49], while some RxLR effectors such as Avr1b-1 and Avr1k can be recognized by NB-LRR immune receptors to confer resistance [50]. Research in 2015 also found that the evolution of RxLR effectors varies between the genus Phytophthora and downy mildews: more conserved RxLR effectors were observed in the genus Phytophthora [51]. We performed a whole-genome scan for RxLR effector genes in race 0 and race 1. The analysis showed that most RxLR genes were distributed in repeat-rich, gene-sparse regions (Fig. 4), suggesting rapid evolution of RxLR effectors. Specifically, a total of 308 RxLR effector genes were predicted in race 0, and 199 genes in race 1 (Additional files 7 and 8). The difference in RxLR effector gene numbers between race 0 and race 1 may contribute to the variance in the aggressiveness of these pathogens in cultivated tobacco.

Crinkler (CRN) effectors are another important class of effectors that cause leaf crinkling in plants [52]. To investigate CRN effectors in P. nicotianae races 0 and 1, we first used EMBOSS getorf (−minsize 300) to extract open reading frames (ORFs) from the whole genome, and then used HMMer (−E 1e-5) with existing profiles [25]. Predicted CRN effectors were filtered by the presence of the LxLFLAK motif. A total of 32 and 26 CRN effectors were annotated in P. nicotianae races 0 and 1, respectively. However, the number of CRN effectors may be underestimated, given the model we used [53].