There are several members of Venturiaceae family known as scab-causing pathogens including Venturia inaequalis, V. nashicola, V. pyrina, V. carpophila, and V. effusa (Gonzalez-Dominguez et al., 2017). Interestingly, these species have distinct and narrow host range (Sivanesan, 1977). Venturia inaequalis, infecting apple but not pear, is the most studied species of Venturiaceae family (Bowen et al., 2011; MacHardy, 1996). Several genome analyses have been reported in recent years, providing a platform for development of scab resistant apple cultivars (Deng et al., 2017; Passey et al., 2018; Shiller et al., 2015). In contrast, significantly less effort was given on other Venturiaceae family members. Recently, draft genome sequence of Venturia pyrina which infects European pear but not Asian pear has been published (Deng et al., 2017). In addition, draft genomes of the peach scab causing fungus V. carpophila and pecan infecting V. effusa have been reported (Bock et al., 2016; Chen et al., 2017). More recently, genome sequences of several V. nashicola isolates have been presented (Johnson et al., 2019; Prokchorchik et al., 2019).

Plant pathogens deliver specific effector proteins into the host plant cells in order to promote virulence. On the other hand, plants have evolved specific immune receptors to recognize corresponding effector proteins directly or indirectly, and activate disease resistance (Jones and Dangl, 2006). To date, many effectors were discovered using genetic map-based cloning approach and proteome analysis in fungal and oomycete plant pathogens (Selin et al., 2016; Uhse and Djamei, 2018). In addition, recent progress on plant pathogen genomics enabled researchers to accelerate prediction of effectors. For example, comparative genomics approach helped to identify the Ave1 effector from a fungal pathogen Verticillium dahliae that is recognized in tomato plants carrying Ve1 disease resistance gene (de Jonge et al., 2012).

In V. inaequalis-Malus pathosystem, 17 genetically characterized disease resistance loci were reported (Bus et al., 2011). Of these, two resistance loci, Rvi6 and Rvi15, were studied in detail but the avirulence determinants from V. inaequalis have not been identified to date (Joshi et al., 2011; Schouten et al., 2014; Vinatzer et al., 2001). The genetic basis of scab resistance in Asian and European pears has not been fully established, although some molecular markers linked to scab resistance have been developed (Bouvier et al., 2012; Won et al., 2014).

In this study, we used the whole genome sequences of several major Venturiaceae species to reconstitute the phylogenetic relationships between them. Furthermore, we predicted the small secreted proteins that may function as virulence or avirulence effectors during host plant infection. The data presented in this study would be valuable to better understand the host-specific virulence mechanism of V. nashicola.

Previously we have sequenced and annotated the genomes of two Venturia nashicola isolates, MAFF615029 and PRI2 collected from Japan and Korea, respectively (Prokchorchik et al., 2019). To define their phylogenetic relationships to the other members of Venturia genus, we aimed to construct a phylogenetic tree consisting of the strains for which whole genome sequence is available. As our study focuses on very closely related Venturia isolates we aimed to use the whole genome level data in order to be able to reliably resolve the closest relatives, which is in contrast to limited gene set alignments used previously for phylogenetic analysis (Yun et al., 2013; Zhao et al., 2012). We used whole genome predicted protein sets of 10 Venturia isolates for constructing phylogenetic tree. The protein sets were either retrieved from Mycocosm database (https://genome.jgi.doe.gov/mycocosm/home) for V. pyrina ICMP11032 and V. inaequalis ICMP13258, or generated in-house using Fungap pipeline applied to genome sequences of V. carpophila JP3-5 (NCBI BioProject: PRJNA321389), V. effusa 3Des10b (NCBI BioProject: PRJNA285422), V. inaequalis 1389 (NCBI BioProject: PRJNA261633), V. inaequalis 1639 (NCBI BioProject: PRJNA261633), V. inaequalis EU-B04 (NCBI BioProject: PRJNA407103), and V. inaequalis 05/172 (NCBI BioProject: PRJNA354841). Protein sequences of all the genes annotated in the Venturia isolates were clustered using OrthoMCl v2.0.9 software (Li et al., 2003) with BLAST e-value 10⁻⁷. Only clusters conserved in all the isolates were considered for phylogenetic analysis (total 6,650 clusters). Hidden Markov Model (HMM) profiles of each gene cluster were built using HMMER3 v3.2.1 (Mistry et al., 2013). These were then compared between strains and neighbor-joining phylogenetic trees were built based on these comparisons using the pHMM-Tree (Huo et al., 2017). Each tree was bootstrapped 100 times and the best tree was used for further analysis. All cluster trees were used to build a majority rule consensus super tree with 70% branch support threshold using Geneious Prime 2019.04 (Fig. 1). Similar to recent studies (Yun et al., 2013; Zhao et al., 2012), the five V. inaequalis isolates were grouped together and V. nashicola isolates analyzed in this study showed close relationship to V. pyrina (Fig. 1). These data indicate that V. pyrina and V. nashicola form a separate clade from V. inaequalis isolates (Fig. 1).

In order to further confirm the genetic relatedness of V. naschicola MAFF615029 and PRI2 isolates, we performed a whole genome comparison. We first aimed to discover single nucleotide polymorphisms (SNPs) using Snippy (v. 4.3.8) (Seemann, 2015). To begin with, comparison between V. nashicola MAFF615029 and PRI2 only was run, then Snippy was run in pairwise manner comparing different Venturia isolates to V. nashicola MAFF615029. Finally, snippy core command was used to call the core SNPs between all the strains compared. We could identify total 105,564 SNPs when comparing the MAFF615029 and PRI2 isolates only. To compare the amount of SNPs between more distant relatives of V. nashicola we generated the core SNPs shared in all the genomes compared to V. nashicola MAFF615029. We could identify 542 core SNPs when comparing MAFF615029 and PRI2, 14,304 SNPs when comparing V. nashicola MAFF615029 and V. pyrina, 35,549 SNPs when comparing V. nashicola MAFF615029 and V. inaequalis ICMP13258, 46,735 SNPs when comparing V. nashicola MAFF615029 and V. effusa and 48,263 SNPs when comparing V. nashicola MAFF615029 to V. carpophila. These results confirmed the close relation of V. nashicola strains between each other suggested by our phylogenetic analysis.

We further performed structural genome variations discovery in PRI2 genome in comparison to MAFF615029 isolate. We used NucMer with following parameters ‘-maxmatch -l 100 -c 500’ (Delcher et al., 2002) and then visualized using Assemblytics (Nattestad and Schatz, 2016). We could identify 777 significant variants. Majority of them were small, mostly 50-500 nucleotides (nt) deletions and, interestingly, insertions, deletions, repeat expansions and contractions were represented relatively equally among more significant (500-10,000 bp) variants (Fig. 2A and B). However, this could be the consequence of short-read assembly of PRI2 isolate genome.

Next, we analyzed V. nashicola MAFF615029 and PRI2 genomes for larger genomic re-arrangements. For that we used Sibelia with following parameters ‘-s fine -q -v --gff --correctboundaries’ (Minkin et al., 2013) and visualized the results using Circos (Krzywinski et al., 2009). There were no large genomic rearrangements detected inside the large scaffolds (Fig. 2C). Overall the genomes surveyed show very high similarity to each other and most of the differences are located in the unaligned short scaffolds. For example, there are 21 predicted genes in V. nashicola PRI2 and there are 586 predicted genes in V. nashicola MAFF615029 located in such scaffolds.

Bioinformatic prediction of effectors may help to identify the virulence and avirulence determinants in fungal species. We used a pipeline for effector annotation based on multiple criteria outlined in Fig. 3. Predicted genes were translated and checked for the presence of signal peptide using SignalP v4.1 (Nielsen, 2017). Proteins with predicted signal peptides were further surveyed for the absence of transmembrane helices using TMHMM (Krogh et al., 2001). All the proteins longer than 400 AA were removed and the short-listed proteins were further assessed for potential targeting in the plant cell using TargetP (Emanuelsson et al., 2007). Resulting set of proteins was subjected to several analyses to score and prioritize effector candidates. In short, the short-listed candidates were subjected to EffectorP2 analysis (Sperschneider et al., 2018) and each effector candidate was assigned a score between 0 and 10 based on the effector prediction confidence. Some fungal and oomycete effectors were shown to be localized to the plant cell nucleus (Ahmed et al., 2018; Vargas et al., 2016) and therefore the short-listed candidates were evaluated for the presence of nuclear localization signal (NLS) using NLStradamus (Nguyen Ba et al., 2009), based on the NLS prediction confidence a score from 0 to 5 was assigned. In addition, short-listed proteins were blasted against Fungal Stress Related Database (FSRD) (Karányi et al., 2013) and PHI-base (Pathogen-Host Interaction database) (Urban et al., 2017) databases. Depending on BLAST hit coverage and similarity level a score from 0 to 10 was assigned for comparison with each of two databases. Finally, short-listed proteins were assessed for the known fungal/oomycete effector motifs, such as RxLR, DEER, RxLX[EDQ], [KRHQSA][DENQ]EL, [Y/W]xC, RSIVEQD (Sonah et al., 2016) using a custom python script and a score from 0 to 6 was assigned. A sum of individual tool scores was used as an integrated effector prediction confidence score and was used to rank the effector candidates in a final short-list. Effector candidates with overall score 0 were discarded from the final short-list. As the results we identified the following numbers of potential effectors in Venturia strains analyzed: V. nashicola MAFF615029, 451; V. nashicola PRI2, 424; V. pyrina ICMP11032, 477; and V. inaequalis ICMP13258, 673 (see Supplementary Tables 1-4, Supplementary FASTA-files 1-4). We focused on these 4 representative isolates in order to identify the putative effectors that are specific to the Venturia species infecting European pear, Asian pear or apple. The putative effectors were clustered using OrthoMCL v2.0.9 with BLAST e-value of 10⁻⁷ and resulting effector families and single effector genes were assessed for presence or absence in the Venturia strains. We identified the core effector repertoire of two V. nashicola strains as 353 conserved effectors (Fig. 4A). When we expanded effector comparison to include two strains of V. nashicola and one V. pyrina strain, we could identify 235 core effectors. In addition to core effectors we identified 118 potential effectors specific to V. nashicola strains and 183 potential effectors specific to V. pyrina (Fig. 4B). These putative effectors may act as host-range determinants between V. nashicola and V. pyrina enabling an avirulent recognition in non-host plant species. Finally, we expanded the effector comparison to include V. inaequalis ICMP13258 putative effectors. In this case we could identify 150 prospective core effectors for Venturia genus and 82 effector families specific for V. nashicola strains infecting Asian pear (see Supplementary Table 5), 134 effector families specific for V. pyrina infecting European pear (see Supplementary Table 6) and 200 effector families specific for V. inaequalis infecting apple (Fig. 4C, Supplementary Table 7). Overall presence and absence patterns of effector families in repertoires of V. nashicola, V. pyrina and V. inaequlis seem to correlate with the host range of the pathogen (Fig. 4D).

Based on our sequencing data we were able to analyze the phylogeny of representative isolates of Venturia genus using the whole genome supertree approach. We also employed the information from multiple gene clusters using the HMM profile comparison instead of protein alignment of single copy genes. Overall, our phylogenetic tree supports the phylogenetic analysis reported previously (Yun et al., 2013; Zhao et al., 2012). We also confirm a close relation between V. nashicola and V. pyrina species.

Fungal effector prediction is one of the most challenging tasks in contemporary plant pathology research. Availability of whole genome assemblies and utilizing a range of computational tools should increase the effector prediction sensitivity. Here we used integrated approach for identifying effectors where results of different analytic tools are evaluated to rank and shortlist the effector candidates for further analysis. We analyzed the presence of effector families in various Venturia species infecting different plant hosts. This provides a valuable resource for future research aiming to confirm effectors that play important roles in host specialization of Venturia species.