Novel Aspects on The Interaction Between Grapevine and Plasmopara viticola: Dual-RNA-Seq Analysis Highlights Gene Expression Dynamics in The Pathogen and The Plant During The Battle For Infection

Mgaloblishvili, a Vitis vinifera cultivar, exhibits unique resistance traits against Plasmopara viticola, the downy mildew agent. This offers the unique opportunity of exploring the molecular responses in compatible and incompatible plant-pathogen interaction. In this study, whole transcriptomes of Mgaloblishvili, Pinot noir (a V. vinifera susceptible cultivar), and Bianca (a resistant hybrid) leaves, inoculated and non-inoculated with the pathogen, were used to identify P. viticola effector-encoding genes and plant susceptibility/resistance genes. Multiple effector-encoding genes were identified in P. viticola transcriptome, with remarkable expression differences in relation to the inoculated grapevine cultivar. Intriguingly, five apoplastic effectors specifically associated with resistance in V. vinifera. Gene coexpression network analysis identified specific modules and metabolic changes occurring during infection in the three grapevine cultivars. Analysis of these data allowed, for the first time, the detection in V. vinifera of a putative P. viticola susceptibility gene, encoding a LOB domain-containing protein. Finally, the de novo assembly of Mgaloblishvili, Pinot noir, and Bianca transcriptomes and their comparison highlighted novel candidate genes that might be at the basis of the resistant phenotype. These results open the way to functional analysis studies and to new perspectives in molecular breeding of grapevine for resistance to P. viticola.

interfere with the host's immune system by delivering cytoplasmic effectors inside the plant cell and secreting apoplastic effectors into the plant extracellular space [19], leading to susceptibility (Effector Triggered Susceptibility-ETS). Effectors are pathogen-encoded proteins that lack of clear sequence similarity to known function, do not always have enzymatic activity and possess, in most cases, a high sequence diversity [20]. Among the known cytoplasmic effectors, encoded by genes expressed by P. viticola, are the RXLR, CRN, and YxSLK classes [21][22][23][24][25][26][27]. The RXLR class of effectors is defined by a conserved N-terminal RXLR motif and a diverse C-terminal domain, responsible for the effector activity inside the host cell [20]. The CRNs ("Crinkler") effectors are widespread across the oomycete lineage, and can cause cell death in the host plant in some cases [28,29]. They consist of modular factors that are secreted and translocated inside host cells by means of a conserved N-terminal domain, with a characteristic LXLFLAK motif, and execute their presumed effector function in the host nucleus through the C-terminal domain [29]. The YxSL[RK] proteins, first discovered in Pythium ultimum Trow [30] and recently found in P. viticola [21], show a modular organization with a conserved amino-terminal region with four conserved motifs and a highly variable C-terminus that make them candidates to be considered as an effector family. Oomycete apoplastic effectors include a large number of hydrolytic enzymes, which are involved in the degradation of host cell components, enabling penetration of host cells [31]. Apoplastic effectors act at the extracellular side by inhibiting the activity of antifungal enzymes, such as proteases and 1,3-glucanases, produced by the plant [32] or can be translocated inside the host cell, through a poorly understood mechanism [33], that can be mediated by haustoria as in the case of EPIC1 in Phytophthora infestans (Mont.) de Bary [34]. Apoplastic effectors known in the P. viticola group include trypsin, elicitin, and NPP1 effectors [21,27]. Serine proteases and specifically trypsin-like enzymes are involved in the response mechanism of oomycetes against plant defenses [35]. Elicitins are structurally conserved extracellular proteins in Phytophthora and Pythium species known to sequester sterols from the host plant, to overcome their inability to synthesize these lipids [31]. NPP1 effectors belong to the necrosis-inducing proteins, known as Nep1-like proteins or necrosisand ethylene-inducing peptide 1-like proteins (NLPs) [20,36]. In addition to the above-cited effectors, proteins possessing a signal peptide for secretion could putatively act as effectors themselves [20].
Effector recognition by the host may occur through direct effector-receptor binding or by sensing the perturbing activity of an effector on host components [37]. Effectors delivered to the host cytoplasm can be recognized by intracellular resistance (R) proteins of the cytoplasmic nucleotide binding site leucine-rich repeat (NBS-LRR) class [38] leading to immunity (Effector Triggered Immunity-ETI). The pathogen growth can be blocked through programmed cell death (hypersensitive response) at the infection site, as occurs in American grapevine species, or by the synthesis of antifungal compounds and structural barriers mediated by a complex response pathway, that involves hormones such as salicylic acid, jasmonic acid, and ethylene [39], as occurs in Mgaloblishvili [17].
To get insight into the resistance mechanism of Mgaloblishvili, a comparative transcriptomic study, where the differentially expressed genes from Mgaloblishvili leaves inoculated or non-inoculated with P. viticola were compared to those of Pinot noir (susceptible control) and Bianca (resistant control harboring Rpv3 QTL of American origin), was undertaken. Since the genomes of Mgaloblishvili and Bianca are still not available, the reads were mapped against predicted mRNAs of a Pinot noir selfing reference transcriptome (PN40024 12X v2). Transcriptomic analysis led to the discovery of unique resistant traits of Mgaloblishvili associated with the overexpression of genes related to pathogen recognition, the ethylene signaling pathway, synthesis of antimicrobial compounds and enzymes, and development of structural barriers [17].
In this study, we used the same transcriptomic data to further dissect the molecular mechanisms underlying the interaction between P. viticola and grapevine. Based on the analysis of P. viticola transcriptome on Mgaloblishvili, Pinot noir, and Bianca, we first aimed at identifying the pathogen genes and pathways involved in the interaction with resistant (Mgaloblishvili and Bianca) and susceptible (Pinot noir) grapevines. Then, we explored the pathways involved in susceptibility or resistance to the pathogen based on network analysis of Mgaloblishvili, Pinot noir, and Bianca transcriptomes. Finally, we performed a de novo assembly of the grapevines transcriptomes to discover novel transcripts, not present in the reference genome, that might be at the basis of the resistant phenotype.

RNA-seq Data
Whole transcriptomes of leaves collected from Mgaloblishvili (M), Pinot noir (P), and Bianca (B) plants inoculated and non-inoculated with P. viticola at 0, 1, 2, and 3 days after inoculation (dai) published in [17] (Study: PRJEB24540; European Nucleotide Archive) were used to i) analyze the P. viticola transcriptome; ii) analyze grapevine coexpression gene patterns; and iii) perform de novo transcriptome assembly. Alignments against the P. viticola genome were performed with the following Tophat command: tophat-G PVIT_GENES_v2.gff3 -segment-length 17 -transcriptome-index \\ PVITGENESv2 -i 20 -I 100,000 -g 10 -p 3 -library-type fr-firststrand \\ PVITv2 R1_trimm.fastq R2_trimm.fastq where PVIT_GENES_v2.gff3 contains all the predicted genes in the P. viticola genome and PVITv2 the genome sequence for the strain sequenced in [21] and are both available at the Plasmopara viticola multi omics dissection project [40]. Here, R1_trimm.fastq R2_trimm.fastq indicates the mates after trimming. Alignments were then filtered to retain only unique mapping reads, by exploiting the flag NH:i:1 that indicates alignments involving reads with only one alignment.
Transcript abundance was quantified by using htseq-count and by setting -s reverse to take into account the correct strandedness of the libraries.
Analysis of differentially expressed genes (DEGs) was performed only for the 3 dai time point, where there is the largest number of reads mapping on the P. viticola genome, for each variety; for this task we considered only genes with at least 20 reads in total over the 9 libraries at 3 dai. We did not apply a more stringent threshold as DESeq2 itself will remove low counts genes; in DESeq2 we used the following thresholds: alpha = 0.1, minimum absolute value of the log2 fold change (LCF) = 1, therefore we will report only differentially expressed genes that have an absolute change of their expression level across two samples of at least 2-fold.

GO Enrichments for P. viticola DEGs
It was evaluated whether the set of DEGs from specific functional categories are enriched in the different samples by using the GO (Gene Ontology) annotation of P. viticola [21] and the webserver AgriGO [41]. Enriched biological processes from this analysis were then reduced by exploiting ReviGO [42] and taking into account only GO processes with a p-value for enrichment smaller than 1E-03.
The annotation of cytoplasmic and apoplastic effectors and genes having a signal peptide for secretion made in [21] was intersected with the differentially expressed genes from this work to identify the transcripts differentially expressed by P. viticola when infecting Mgaloblishvili and Pinot noir (PvM_vs_PvP), Mgaloblishvili and Bianca (PvM_vs_PvB), and Pinot noir and Bianca (PvP_vs_PvB). Three-way Venn diagrams showing the overlap among different contrasts were built using the jVenn web-server [43], using DEG lists as input data. The function heatmap.2 of gplots R package [44] was used for the graphical representation of gene expression and hierarchical clustering.

Network Analysis of Grapevine Transcriptome: Gene Coexpression Network Analysis
Matrix of raw count data of gene expression including all replicates was imported to DESeq2 [45] and a model.matrix was defined according groups. After normalization by library size and RLE (Relative Log Expression) methods a matrix of average by group was prepared in DESeq2. After filtering out the low expressed genes, resulted combined gene expression matrix was used for gene coexpression network analysis using the log2-transformed expression values of the most variable genes (3000) and the WGCNA package v1.51 [46]. Briefly, a softpower β was chosen using the function pickSoftThreshold to fit the signed network to a scale-free topology. Next, an adjacency matrix was generated. Topological Overlap Matrix (TOM) was used as an input in the function hclust ("average" method) to construct a hierarchical clustering tree (dendrogram). TOM is a measure of topological similarity between the genes within a network, i.e., it evaluates whether two or more nodes share links within the network and groups them into the same module [47,48]. A threshold of 0.15 (correlation > 85%) was chosen to merge similar modules, and only modules having at least 20 genes were kept.

Network Analysis of Grapevine Transcriptome: Functional Enrichment Analyses
GO and pathway enrichment analyses of gene sets (within networks or modules) were performed by in-house scripts using GOseq and KEGG as the data sources. p-values were determined through a Fisher Exact test and adjusted via the Benjamini-Hochberg's method. A threshold of adjusted p-value < 0.05 was used to determine the statistical significance of enrichment results.

Network Analysis of Grapevine Transcriptome: Transcription Factor Enrichment
To identify potential common transcription factors (TFs) that may control transcription of module genes, the promoter region of the genes of each module were investigated for transcription factor binding sites (TFBS) using Plant TFDB database and TF enrichment tools and a hypergeometric test was applied to enrich the statistically significant elements and related TFs (p-value ≤ 0.05). p-values were adjusted for multiple testing corrections using Benjamini and Hochberg's false discovery rate (FDR) method (p ≤ 0.05).

De novo Transcriptome Assembly of Grapevine Cultivars
The raw RNA-seq datasets for the 63 paired ends (PE50) libraries available for Pinot noir (P), Bianca (B), and Mgaloblishvili (M) were retrieved and prepared for a de novo assembly approach. For each genotype, three biological replicates were available for the non-inoculated plants (ni) at zero dai (ni0dai; untreated control plants), day one (ni1dai), day two (ni2dai), and day three (ni3dai), and the inoculated plants (i) at day one (i1dai), day two (i2dai), and day three (i3dai). The triplicate libraries of each sample were merged in one single dataset, respectively, but keeping right and left reads into two separate files. The merged datasets were used for the de novo assembly procedure performed with Trinity [49]. The parameters used were -seqType fq -max_memory 60G -left -right -include_supertranscripts -SS_lib_type RF -CPU 4 -trimmomatic. The assemblies were run on a HPC cluster [50].
The de novo assembled transcriptomes were then merged in a single dataset and a clustering was performed using cd-hit-est (c = 0.95) [51]. The clusters formed only by transcripts belonging to either Bi, Bni, Mi, Mni, or combinations of these but not clustering with Pi or Pni were selected. Similarly, clusters formed only by Pi or Pni or combinations of these were retrieved separately. The cluster consensus sequences with a length ≥ 300 bp were retrieved obtaining a dataset of transcripts for each genotype. These transcripts were then subjected to blastx (evalue 1E-7) against a dataset formed by all the plant proteins having a UniProt database entry.
The blastx outputs were filtered, considering only the hits having a ratio (alignment length/subject length) ≥ 0.8 and pident ≥ 70.00. The UniProt IDs of the retrieved blastx hits for B and M were compared with the IDs retrieved for Pi/Pni for both the inoculated and non-inoculated state, to find unique IDs characterizing Bi, Bni, Mi, Mni, Pi, and Pni. The sequences corresponding to the transcripts associated with the IDs obtained were subjected to Gene Ontology and enrichment analyses using the Blast2GO suite (version 5.2.5) [52].

Statistical Analysis of Differentially Expressed De novo Transcripts
The de novo transcriptomes generated for each genotype and condition were clustered separately using cd-hit-est with the parameter c = 0.95. The RNA-seq reads passing the quality filter set in the Trinity step and used to assembly the transcriptome de novo were re-aligned on their respective clustered transcriptome. Bowtie2 (version 2.3.4.3) [53] was used to align reads with the default parameters and -L 31. The read count number aligning on each transcript was calculated using an in-house script counting the alignments having the SAM flag 99/147, 83/163, 81/161, and 97/145. A read count matrix was generated and used for the subsequent analysis. DESeq2 R package [45] was used to perform DEG analysis, by a multifactor design method, among the different treatments (non-inoculated vs. inoculated) of each variety at each time point. Transcripts with more than 5 reads were retained for the analysis. For each transcript, the log2 fold change (FC), p-value and adjusted p-value were evaluated, and only transcript IDs with false discovery rate (FDR)-adjusted p-value < 0.1 were retained. Three-way Venn diagrams showing the overlaps among different time points per genotype were built using the jVenn web-server [43], using DEG lists as input data.

Availability of Data and Materials
The datasets for this study can be found in the European Nucleotide Archive (ENA) under the study accession number: PRJEB24540.

P. viticola Transcriptome Analysis
The analysis of P. viticola transcriptome on Mgaloblishvili (PvM), Pinot noir (PvP), and Bianca (PvB) showed that 5251 genes were expressed by the pathogen at 3 dai. Further analysis highlighted the presence of 1519 DEGs when infecting Mgaloblishvili in comparison with Pinot noir (PvM_vs_PvP), Mgaloblishvili in comparison with Bianca (PvM_vs_PvB), and Pinot noir in comparison with Bianca (PvP_vs_PvB; Table 1). Only two genes were differentially expressed by P. viticola during infection on Mgaloblishvili with respect to Pinot noir (PvM_vs_PvP), while 533 and 984 genes were differentially expressed by the pathogen when comparing the inoculation of Mgaloblishvili and Bianca (PvM_vs_PvB) and Pinot noir and Bianca (PvP_vs_PvB), respectively (Table 1). P. viticola DEGs were mostly up-regulated (96%) in Mgaloblishvili and Pinot noir with respect to Bianca (Table 1; Table S1). To identify how the P. viticola transcriptome is modulated upon infection of grapevine accessions with different susceptibility levels, we mainly focused on P. viticola DEGs specific for each contrast (PvM_vs_PvP, PvM_vs_PvB, and PvP_vs_PvB) or shared among contrasts (PvM_vs_PvP with PvP_vs_PvB; PvM_vs_PvB with PvP_vs_PvB; PvM_vs_PvP with PvM_vs_PvB; and PvM_vs_PvP with PvP_vs_PvB and PvM_vs_PvB). The number of DEGs within and among contrasts is reported in Figure 1. Two P. viticola DEGs were found in the comparison between Mgaloblishvili and Pinot noir (PvM_vs_PvP), one specific and one shared in the contrast Pinot noir/Bianca (PvP_vs_PvB). The comparison of P. viticola transcriptome on Mgaloblishvili and Bianca (PvM_vs_PvB) and Pinot noir and Bianca (PvP_vs_PvB) resulted in the presence of 13 and 463 specific DEGs, respectively. Five hundred and twenty DEGs were common in the comparisons Mgaloblishvili/Bianca and Pinot noir/Bianca (PvM_vs_PvB with PvP_vs_PvB), while no DEGs were shared among the three contrasts (PvM_vs_PvP with PvP_vs_PvB and PvM_vs_PvB) nor between PvM_vs_PvP with PvM_vs_PvB.
Differences between biological processes over-represented in the two contrasts are mainly related to "response to stress" (GO:0006950) in PvM_vs_PvB, which was not represented in PvP_vs_PvB.

P. viticola Effectors
The transcripts encoding for proteins with known or putative P. viticola effector function were identified among the P. viticola DEGs at 3 dai. Thirty-five cytoplasmic and apoplastic effectors and seventy-six P. viticola genes encoding proteins possessing a signal peptide for secretion were differentially expressed by P. viticola when comparing the pathogen transcriptomes on Mgaloblishvili and Pinot noir (PvM_vs_PvP), Mgaloblishvili and Bianca (PvM_vs_PvB), and Pinot noir and Bianca (PvP_vs_PvB; Table S3; Figure 2; Figure 3). A single P. viticola gene (PVITv1009470), encoding a trypsin-like apoplastic effector, was differentially expressed only during the interaction with Mgaloblishvili compared to Pinot noir (PvM_vs_PvP; Figure 2). Four P. viticola genes encoding apoplastic effectors (the elicitin PVITv1020900 and trypsin PVITv1006768, PVITv1031227, and PVITv1029052; Figure 2) and two genes encoding proteins possessing a signal peptide for secretion (PVITv1033624, and PVITv1005898; Figure 3) were up-regulated only in Mgaloblishvili compared to Bianca (PvM_vs_PvB). In general, none of the transcripts associated with elicitin and NPP1 effectors showed any known function. PVITv1033624 and PVITv1005898 encoded enzymes involved in β-glucan hydrolysis (glucanendo-1,6-beta-glucosidase) and pectin de-esterification (pectinesterase).
Genes 2020, 11, 261 8 of 24 Differences between biological processes over-represented in the two contrasts are mainly related to "response to stress" (GO:0006950) in PvM_vs_PvB, which was not represented in PvP_vs_PvB.
Forty-five transcripts, corresponding to four cytoplasmic effectors (three RxLR: PVITv1002967, Forty-five transcripts, corresponding to four cytoplasmic effectors (three RxLR: PVITv1002967, PVITv1021251, PVITv1008294; and one YxSLK: PVITv1028596), seven apoplastic effectors (from PVITv1011483 to PVITv1028723), and 34 transcripts coding for secreted proteins (from PVITv1025751 to PVITv1038141), were specifically associated with the pathogen when infecting the V. vinifera cv. Pinot noir in the comparison with Bianca (PvP_vs_PvB). The YxSLK effector encodes a spermidine synthase, while the others are completely uncharacterized, as none of them have significant homology with characterized proteins. Among the genes encoding for secreted proteins, some were linked to oxidative stress (6-phosphogluconolactonase, monodehydroascorbate reductase NADH) and nutrition (acid phosphatase).
The above cited genes were all up-regulated with a LFC (log2 Fold Change) greater than 3. Only two effectors (PVITv1021785 and PVITv1036189), encoding uncharacterized proteins, were down-regulated in the two V. vinifera varieties compared to Bianca (PvM_vs_PvB and PvP_vs_PvB).

Grapevine Transcriptome Network Analysis
The network analysis was performed on 3000 genes with the highest gene expression variation across samples. Genes with similar expression pattern were grouped into modules. A total of 10 coexpression modules were identified via hierarchical clustering (Figure 4).
Pinot noir in the comparison with Bianca (PvP_vs_PvB). The YxSLK effector encodes a spermidine synthase, while the others are completely uncharacterized, as none of them have significant homology with characterized proteins. Among the genes encoding for secreted proteins, some were linked to oxidative stress (6-phosphogluconolactonase, monodehydroascorbate reductase NADH) and nutrition (acid phosphatase).
The above cited genes were all up-regulated with a LFC (log2 Fold Change) greater than 3. Only two effectors (PVITv1021785 and PVITv1036189), encoding uncharacterized proteins, were downregulated in the two V. vinifera varieties compared to Bianca (PvM_vs_PvB and PvP_vs_PvB).

Grapevine Transcriptome Network Analysis
The network analysis was performed on 3000 genes with the highest gene expression variation across samples. Genes with similar expression pattern were grouped into modules. A total of 10 coexpression modules were identified via hierarchical clustering (Figure 4). The number of genes per module (module size) ranged from 40 (purple) to 977 (turquoise). Gene Set Enrichment analyses were performed on the modules to identify enriched GO terms or pathways (Table S4), showing that all the modules, except the purple one, had significant enrichments, providing evidence of functional coherence within modules. Some of the most highly significant GO terms and KEGG pathways were: "monoterpenoid biosynthesis" and "glutathione metabolism" (green module); "alpha-linolenic acid metabolism", "biosynthesis of secondary metabolites", and "flavonoid biosynthesis" (black module); "glucosinolate biosynthesis", "cyanoamino acid metabolism", and "plant-pathogen interaction" (blue module); "flavone and flavonol biosynthesis" and "phagosome and cutin, suberine, and wax biosynthesis" (brown module); "linoleic acid metabolism" (pink module); "plant hormone signal transduction" (red module); "phenylpropanoid biosynthesis" and "flavonoid biosynthesis" (turquoise module); and "plant-pathogen interaction" and "monoterpenoid biosynthesis" (yellow module).
All genes of each module were screened to identify transcription factors as main regulators of modules: turquoise module includes the most TFs (Table S5). In addition, the analysis of the conserved Transcription Factor Binding Sites (TFBSs) on promoter regions of genes included in each module revealed that a total of 299 TFBSs (p-value < 0.05) were over-represented across modules, except in the purple and magenta modules. Each module was enriched in minimum 14 (green The number of genes per module (module size) ranged from 40 (purple) to 977 (turquoise). Gene Set Enrichment analyses were performed on the modules to identify enriched GO terms or pathways (Table S4), showing that all the modules, except the purple one, had significant enrichments, providing evidence of functional coherence within modules. Some of the most highly significant GO terms and KEGG pathways were: "monoterpenoid biosynthesis" and "glutathione metabolism" (green module); "alpha-linolenic acid metabolism", "biosynthesis of secondary metabolites", and "flavonoid biosynthesis" (black module); "glucosinolate biosynthesis", "cyanoamino acid metabolism", and "plant-pathogen interaction" (blue module); "flavone and flavonol biosynthesis" and "phagosome and cutin, suberine, and wax biosynthesis" (brown module); "linoleic acid metabolism" (pink module); "plant hormone signal transduction" (red module); "phenylpropanoid biosynthesis" and "flavonoid biosynthesis" (turquoise module); and "plant-pathogen interaction" and "monoterpenoid biosynthesis" (yellow module).
All genes of each module were screened to identify transcription factors as main regulators of modules: turquoise module includes the most TFs (Table S5). In addition, the analysis of the conserved Transcription Factor Binding Sites (TFBSs) on promoter regions of genes included in each module revealed that a total of 299 TFBSs (p-value < 0.05) were over-represented across modules, except in the purple and magenta modules. Each module was enriched in minimum 14 (green module) and maximum 89 (blue module) TFBSs (Table S6). These findings indicated that WGCNA software could efficiently exploit gene expression levels to build modules able to recover coregulation.
Gene expression pattern for each module was extracted from the network analysis to determine the genes discriminating for resistance/susceptibility to P. viticola. Only the yellow and the blue module genes were selected for patterning since they included most of the genes involved in the plant-pathogen interaction pathways (Table S4): the yellow module mainly associated with resistance to P. viticola in Bianca, while the blue module related to the response of Mgaloblishvili and Bianca to the pathogen. In the yellow module, five main clusters have been identified and three of them showed a different expression pattern among cultivars. Mgaloblishvili and Pinot noir showed similar expression patterns, different from those of Bianca (where the majority of genes were up-regulated; Figure S1). This module included key genes that mimic Bianca response to the pathogen, such as genes encoding for receptor kinase proteins, PR proteins, or genes related to the signaling pathway, production of ROS (Reactive Oxygen Species) and antimicrobial compounds. In the blue module, four main clusters have been identified. Most of them showed a similar expression pattern among the three cultivars (mainly up-regulated), other showed a similar expression pattern between Bianca and Pinot noir (mainly down-regulated), different from the expression pattern of Mgaloblishvili ( Figure S2). This module included genes related to unique Mgaloblishvili resistant behavior, such as genes encoding for protein ubiquitination, mitogen-activated protein kinases, wall-associated receptor kinases, cellulose, and lignin biosynthesis proteins. Table 2 summarizes the characteristics of the libraries used for each genotype (M, Mgaloblishvili; P, Pinot Noir; B, Bianca) and treatment (i = inoculated; ni = non-inoculated). The raw reads ranged from 5,042,556 to 35,478,918 with an average per triplicate ranging from 20,481,001 (Mni at 3 dai) to 29,590,959 (Mi at 3 dai). The reads available for the de novo transcriptome assembly that passed the trimmomatic filtering step ranged from 59,610,375 (Mni at 3 dai) to 86,910,171 (Mi at 3 dai). In total, 21 de novo transcriptomes were assembled ( Table 2). The number of transcripts was comprised between 38,749 (Pni at 0 dai) and 60,459 (Bi at 3 dai) whereas the number of supertranscripts ranged between 28,818 (Pni at 0 dai) and 37,891 (Bi at 3 dai).

De novo Grapevine Transcriptome Assembly
The transcriptomes were merged and clustered using cd-hit-est producing 171,743 clusters. A selection on the clusters was carried out to retrieve only those formed by transcripts belonging to a single genotype at the three time-points together.
The selection resulted in 30,201 clusters for Mi and/or Mni, 26,756 clusters for Pi and/or Pni, and 50,515 clusters for Bi and/or Bni. In addition, clusters containing transcripts belonging to Mi, Mni, Bi, and Bni but not Pi and/or Pni were also selected resulting in 47,020 clusters. Subsequently, the sequences of each selected cluster subset were filtered retrieving only those with a length ≥300 bp to be used in the blastx analysis. The transcripts obtained in each subset were: 16,225 for Mi/Mni, 14,000 for Pi/Pni, and 28,185 for the Bi/Bni. Moreover, 12,579 and 17,258 transcripts were retrieved for Mi/Mni and Bi/Bni, respectively, from clusters formed by transcripts of Mi/Mni and Bi/Bni in any combination.
This analysis focused on differences at nucleotide level and length between the transcripts assembled for each genotype. The aim was to highlight those transcripts differing among the transcriptomes that could be considered as typical and characteristic (unique) for each genotype investigated. However, transcripts differing at nucleotide level with the parameters used in our analysis could have higher homology at protein level. Therefore, a blastx analysis was performed and the Uniprot IDs (IDs) of the hits with pident ≥70.00 and a ratio alignment/subject length ≥0.8 were extracted separately for each genotype and treatment (inoculated or non-inoculated).  Gene Ontology and enrichment analyses were carried out to identify the most enriched GO categories in inoculated vs non-inoculated samples per cultivar (Table S7, Figures S3-S5). Based on the biological process GO categories, the most enriched terms were "phosphorylation" and "oxidation-reduction" for Mgaloblishvili (16% and 14% respectively) and Bianca (15% and 16% respectively); and "cellular protein modification" and "oxidation-reduction" for Pinot noir (16% and 12% respectively). The most enriched cellular component GO category was "integral component of membrane" per cultivar. While, about the molecular function GO categories, Mgaloblishvili showed an enrichment of "nucleic acid binding" terms (14%), Pinot noir an enrichment of "protein binding" terms (22%), and Bianca inoculated samples "hydrolase activity" (18%).

P. viticola Transcriptome
In a previous study [17], the authors showed that Mgaloblishvili, a Georgian (Caucasus) cultivar of V. vinifera, is able to reduce growth and sporulation of P. viticola, compared to the susceptible V. vinifera control (Pinot noir), by four and ten times, respectively. In contrast, the interspecific hybrid Bianca undergoes hypersensitive response upon pathogen inoculation. Here, the transcriptome data (RNA-seq) previously obtained [17] were used to explore candidate pathways and genes that are differentially expressed by P. viticola during the interaction with two V. vinifera varieties characterized by susceptible (Pinot noir) and resistant (Mgaloblishvili) phenotypes and a resistant hybrid (Bianca).
Transcript abundance levels within an organism follow a log-normal distribution, meaning that a small percentage of all transcripts cover a large fraction of all the transcripts present in the cell at any given moment [21]. For this reason, transcriptomic studies require deeper sequencing in order to identify transcripts with lower than average expression levels. It should be noted that in dual-omics experiments, in which the transcriptomes of two organisms are screened at once, the problem is exacerbated. Thus, when mapping RNA reads on the P. viticola genome, only 5251 genes appear to be expressed, and most of them at very low levels. Nonetheless, the samples relative to the 3 dai time point contain enough reads to run differential expression analysis. Almost all of the P. viticola DEGs (99%) concerned the comparison of its transcriptome in Mgaloblishvili and Pinot noir versus Bianca, since the samples concerning P. viticola in Mgaloblishvili and Pinot noir are much more similar. The great similarity between the genes expressed by P. viticola in Mgaloblishvili and Pinot noir could be due to the fact that these are both V. vinifera varieties, while Bianca is characterized by the presence of 20% genes of American grapevine species [54]. Therefore, the comparison of P. viticola transcriptome on V. vinifera with a non-vinifera variety was mainly used to elucidate the genes involved in the plant-pathogen interaction. The biological functions associated with P. viticola DEGs are almost identical in Mgaloblishvili and Pinot noir and are mostly related to biosynthetic processes. Based on PvM_vs_PvB and PvP_vs_PvB contrasts, the P. viticola DE transcriptome on Mgaloblishvili differs from the one on Pinot noir mainly for the presence of GO terms related to response to stress, which could be associated with the action of the plant on the pathogen: at 3 dai P. viticola growth was, in fact, clearly hampered in Mgaloblishvili, while no alteration of the vegetative structures of the pathogen was observed in Pinot noir [17].

P. viticola Effectors
When analyzing the genes encoding effectors expressed by P. viticola in the comparison among grapevine accessions at 3 dai, it was possible to see that no CRN-encoding genes were differentially expressed by the pathogen, while eight transcripts with RXLR and YxSLK motif were up-regulated in both Mgaloblishvili and Pinot noir, and in Pinot noir compared to Bianca (PvP_vs_PvB). Among genes encoding apoplastic effectors, we found few NPP1 (three) and elicitin (seven) DEGs and several trypsin DEGs (17), associated with protein hydrolysis, up-regulated in Mgaloblishvili and/or Pinot noir. Finally, 76 DEGs encoding proteins with a signal peptide for secretion were found. Most of the differentially expressed genes encoding for effectors were up-regulated by the pathogen in both Mgaloblishvili and Pinot noir (57) or in PvP_vs_PvB (45) and only a few were up-regulated exclusively in Mgaloblishvili (seven) or Bianca (two).
Of the 15 apoplastic effectors and 38 secreted proteins encoded by genes with a known function up-regulated by P. viticola in both V. vinifera cultivars, many are proteins involved in protein hydrolysis and modification (glycosylation and folding), as previously found by other authors [21,27]. Another important category is formed by cell wall degrading enzymes (CWDEs) targeting cellulose, hemicellulose, and pectins [55], which constitute the primary defense of the plant cell, and β-1,3-glucans such as callose, which is produced by the plant in response to the pathogen as a physical barrier to penetration [56,57]. CWDEs are extracellular effectors that degrade a wide range of polysaccharides and glycoproteins not only during plant penetration but also to allow the release of nutrients [55].
Five genes encoding apoplastic effectors and two genes encoding secreted proteins were up-regulated by P. viticola exclusively during the interaction with Mgaloblishvili. These effectors are putatively involved in the hydrolysis of antifungal proteins produced by the host cell and cell wall degradation during colonization [58]. While β-(1,6)-glucanase activity has been detected during fungal pathogenesis, its role remains unclear as host plants lack its β-(1,6)-glucan substrate [59]. An explanation of the up-regulation of this gene could be related to the growth of P. viticola mycelium: β-glucanases play a role in cell wall growth and extension by partially hydrolyzing localized areas and enabling insertion of new cell wall material without disrupting its overall integrity [59]. Three days after inoculation on Mgaloblishvili, P. viticola indeed showed contorted, hyper-branched hyphae that could indicate a high cell wall synthesis [17]. These candidate effectors should be better investigated since they are specifically related to the immune response of the resistant V. vinifera cultivar.
Based on the comparative analysis of the P. viticola transcriptome, it was possible to identify four candidate cytoplasmic effectors with a specific role in pathogenesis on Pinot noir: three RXLR and an YxSLK. Of these effectors, only the YxSLK had a function assigned to it: it encodes a spermidine synthase, a protein catalyzing the biosynthesis of a polyamine involved in several fungal functions, such as dimorphism, spore germination and appressorium formation, and in virulence [60,61]. The transcripts encoding non-cytoplasmic effectors associated with susceptibility in Pinot noir showed functions analogous to those previously described and included also genes encoding enzymes involved in response to oxidative stress (monodehydroascorbate reductase (NADH), responsible for scavenging ROSs in plants and reducing phenoxyl radicals [62,63]), metabolism (6-phosphogluconolactonase, involved in the pentose phosphate pathway, a pathway highly represented in P. infestans hyphae [64]), and phosphate nutrition (acid phosphatase [65]).

Networking Analysis Reveals A Putative Gene of Susceptibility in V. vinifera
Network analysis clustered genes in 10 coexpression modules (Figure 4). These modules showed enrichment of the most common plant-pathogen interaction pathways, such as plant hormone signal transduction and biosynthesis of secondary metabolites (monoterpenes and flavonoids) [66,67], confirming the results of the previous study [17].
The availability of a resistant V. vinifera accession provided us the opportunity to explore the genes associated with susceptibility within the species that, to the best of our knowledge, have not been discovered yet [68]. In plant-pathogen interaction, apart from R (resistance) genes, that can detect effectors and activate ETI, there are S (susceptibility) genes, which are required for successful pathogen infection and are considered essential for compatible plant-pathogen interaction [8]. S genes can be involved in early pathogen establishment, modulation of host defenses, and pathogen sustenance [68], and their disruption can confer disease resistance [8]. S genes are genes coding for effector targets that function as negative defense regulators or susceptibility factors [69]. One of the methods used for individuating S genes is based on gene expression studies, where inoculation with the pathogen results in the gene up-regulation during compatible interaction and transgenic overexpression of the gene leads to susceptibility [70]. It must be pointed out that S genes can show a constitutive expression and an increased expression upon pathogen inoculation [71]. To identify putative S genes in V. vinifera, the following hypothesis has been made: the putative S genes should have been down-regulated in Mgaloblishvili (resistant V. vinifera) inoculated samples in comparison with non-inoculated samples, and up-regulated in Pinot noir (susceptible V. vinifera) inoculated versus non-inoculated samples at 1 dai, the most informative time point in P. viticola-grapevine interaction [17]. Among the genes examined, only a gene encoding a LOB domain-containing (LBD) protein (XM_010660137.1) showed these characteristics. In Arabidopsis, the root-specific LBD transcription factor AtLBD20 acts as a repressor of a subset of jasmonate-mediated defense mechanisms during infection of roots with Fusarium oxysporum Schltdl., and the knockout of this gene resulted in increased resistance to F. oxysporum [72]. The down-regulation of this gene in Mgaloblishvili and its up-regulation in Pinot noir make it a candidate S gene associated with a negative regulation of immunity signals [8,69]. Further studies are needed to confirm this hypothesis.

De novo Transcriptome Assembly Reveals New Grapevine Traits Related to P. viticola Response
V. vinifera is one of the most highly heterozygous species and has high genetic diversity [73] such that the commonly used [17,[74][75][76][77] reference genome PN40024 12X v2 [78] is not comprehensive of the entire complexity of modern cultivars. For this reason, de novo transcriptome assembly of Mgaloblishvili, Pinot noir, and Bianca was performed to characterize the different responses to P. viticola infection. The largest number of unique transcripts was detected in Bianca samples, as expected due to its non-vinifera background [54], while the lowest was in Pinot noir (Table S7).
The GO analysis showed that most of the unique transcripts of Mgaloblishvili inoculated samples were involved in the processes of "phosphorylation", "integral component of membrane", and "nucleic acid binding" (Figure S3), that can be related to signaling, recognition, and transcriptional regulation in the plant-pathogen interaction [79][80][81][82][83].
In Pinot noir, most of the unique transcripts of inoculated samples were involved in "cellular protein modification" and "protein binding" ( Figure S4) with a significant role in the plant-pathogen interaction [83].
Unlike Mgaloblishvili, Bianca showed a higher enrichment of GO categories related to the "oxidation-reduction biological process" and "hydrolase activity" ( Figure S5) due to the ability of Bianca to trigger a localized HR and secrete hydrolytic enzymes (such as PR proteins, glucanases, chitinases, and proteases) soon after the onset of infection to counteract the pathogen invasion [56,[84][85][86].

Transcriptomic Changes Occurring In Grapevine During The First Stages of P. viticola Infection
As already found in similar studies [12,87], 1 dai was the most informative time point in our experimental conditions (Table 3; Tables S8-S10) [17]. This result is strictly related to the infection process, when the haustorium and the plant cell membrane get in contact for the first time [17,56,88]. The number of DEGs was larger for Bianca samples, followed by Mgaloblishvili and Pinot noir samples. This trend is correlated to the genotype resistance level against the pathogen. Bianca is a resistant variety showing HR [56,84,85] and promptly (at 1 dai) activates the molecular mechanisms preventing mycelial growth [89]. Mgaloblishvili showed a down regulation of P. viticola growth and sporulation associated with PAMP, DAMP, and effector recognition [17]. Whereas Pinot noir is a susceptible variety, not able to counteract pathogen growth [17,87].

Unique DEGs of Mgaloblishvili in Response to P. viticola Infection
To elucidate the Mgaloblishvili response to P. viticola, unique DEGs having a LFC value above the threshold were discussed (Table S8). Among up-regulated DEGs, five genes, encoding two receptor-like kinase (RLKs) proteins, two 18.3 kDa heat shock proteins, and valencene synthase, were strictly linked to the plant-pathogen interaction mechanism [17]. RLKs (receptor-like kinases) appear to play a central role in the plant-pathogen interaction, signaling during pathogen recognition, and activation of the plant defense mechanisms [90]. Synthesis of small heat shock proteins (sHSPs), with a molecular mass of 15-42 kD, can be induced by heat shock and various abiotic or biotic stresses [91,92] and are involved in the initial stages of fungal recognition in HR-independent defense mechanism [93]. Valencene synthase (already identified in previous study [17]) catalyzes the conversion of farnesyl diphosphate to valencene, a sesquiterpene with antimicrobial activity [94][95][96].
Among the Mgaloblishvili down-regulated DEGs, only one transcript appeared to be strictly linked to the plant-pathogen interaction mechanism, and not already identified in the previous study [17]. It encodes a probable xyloglucan endotransglucosylase/hydrolase (XTH) protein 32. XTHs are a class of cell wall-modifying proteins involved in cell growth, having the ability to loosen cell walls [97,98]. The overexpression of these genes makes the membrane more permeable [99]. Membrane fluidity has a role in the P. viticola-grapevine interaction, where a more rigid membrane is associated with resistance to the pathogen [100]. Therefore, the down-regulation of a XTH gene could likely be a part of the resistance mechanism of Mgaloblishvili by increasing the compactness of cells and membrane permeability.

Conclusions
In conclusion, we identified the dynamics of gene and pathway regulation in P. viticola and susceptible/resistant grapevines during the battle for infection. On the one hand, during the interaction with Mgaloblishvili, P. viticola specifically up-regulated genes encoding apoplastic effectors and putative proteases that were not differentially expressed during the interaction with the susceptible (Pinot noir) and resistant (Bianca) grapevine cultivars. On the other hand, Mgaloblishvili confirmed its unique response to P. viticola infection based on the up-regulation of resistance genes involved in signaling, induction of basal immune response, and metabolism of terpenes, as previously described [17]. Moreover, the de novo assembly of the grapevine accessions allowed us to identify new genes that can be at the basis of V. vinifera resistance to P. viticola.
Apart from the detection of genes putatively involved in resistance, this study also highlighted the presence of a possible susceptibility gene, a transcription factor that was down-regulated in the resistant V. vinifera cv. Mgaloblishvili, up-regulated in the susceptible V. vinifera cv. Pinot noir, and not differentially expressed in the resistant hybrid Bianca. New perspectives in genetic improvement of grapevine for resistance to P. viticola are needed, and in this regard the exploitation of susceptibility genes could be remarkable. Recently, disrupting of S genes to interfere with the host-pathogen compatibility and obtaining disease resistance has become of great interest, due to the accomplishment of this task in a transgene-free way via new genome editing tools [101].
The results obtained in the present study contribute to elucidate the mechanisms of interaction between P. viticola and V. vinifera during infection. Further analyses are needed to investigate the role of the candidate effectors in pathogenicity through functional analysis and to confirm the role of the putative resistance/susceptibility genes of V. vinifera.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/11/3/261/s1, Table S1. List of transcripts differentially expressed by P. viticola upon inoculation of Mgaloblishvili in comparison with Pinot noir (PvM_vs_PvP), Mgaloblishvili in comparison with Bianca (PvM_vs_PvB), and Pinot noir in comparison with Bianca (PvP_vs_PvB), three days after inoculation, Table S2. List of the GO terms related to upand down-regulated differentially expressed genes expressed by P. viticola upon inoculation of Mgaloblishvili in comparison with Bianca (PvM_vs_PvB) and Pinot noir in comparison with Bianca (PvP_vs_PvB), three days after inoculation, Table S3. List of the differentially expressed genes encoding cytoplasmic and apoplastic effectors and proteins with signal peptide for secretion up-and down-regulated by P. viticola upon inoculation of Mgaloblishvili in comparison with Pinot noir (PvM_vs_PvP), Mgaloblishvili in comparison with Bianca (PvM_vs_PvB), and Pinot noir in comparison with Bianca (PvP_vs_PvB), three days after inoculation, Table S4. Enriched pathway terms in modules detected via the network analysis represented in Figure 4, Table S5. List of transcription factors (TFs) detected in each module of network analysis represented in Figure 4, Table S6. Transcription factors (TFs) enrichment of module genes detected via network analysis represented in Figure 4, Table S7. List of transcripts of three grapevine cultivars (Mgaloblishvili, Pinot noir, and Bianca), inoculated and non-inoculated with P. viticola used for GO analysis, Table S8. List of unique differentially expressed genes identified in Mgaloblishvili inoculated with P. viticola at 1, 2, and 3 days after inoculation. dai = day after inoculation, Table S9. List of unique differentially expressed genes identified in Pinot noir inoculated with P. viticola at 1, 2, and 3 days after inoculation. dai = day after inoculation, Table S10. List of unique differentially expressed genes identified in Bianca inoculated with P. viticola at 1, 2, and 3 days after inoculation. dai = day after inoculation, Figure S1. Heatmap of normalized expression values of genes of yellow module detected by network analysis (Figure 4) for Mgaloblishvili, Pinot noir, and Bianca inoculated and non-inoculated with P. viticola at 1, 2, and 3 days after inoculation, Figure S2. Heatmap of normalized expression values of genes of blue module detected by network analysis (Figure 4) for Mgaloblishvili, Pinot noir, and Bianca inoculated and non-inoculated with P. viticola at 1, 2, and 3 days after inoculation, Figure S3. Overview of the most enriched GO categories in inoculated vs. non-inoculated samples of Mgaloblishvili. A: Biological process GO categories; B: cellular component GO categories; C: molecular function GO categories, Figure S4. Overview of the most enriched GO categories in inoculated vs. non-inoculated samples of Pinot noir. A: Biological process GO categories; B: cellular component GO categories; C: molecular function GO categories, Figure S5. Overview of the most enriched GO categories in inoculated vs. non-inoculated samples of Bianca. A: Biological process GO categories; B: cellular component GO categories; C: molecular function GO categories, Figure S6. Comparison between differentially expressed genes (DEGs) of Mgaloblishvili samples inoculated with P. viticola, at 1, 2, and 3 days after inoculation. M = Mgaloblishvili; i = inoculated; dai = day after inoculation, Figure S7. Comparison between differentially expressed genes (DEGs) of Pinot noir samples inoculated with P. viticola, at 1, 2, and 3 days after inoculation. P = Pinot noir; i = inoculated; dai = day after inoculation, Figure S8. Comparison between differentially expressed genes (DEGs) of Bianca samples inoculated with P. viticola, at 1, 2, and 3 days after inoculation. B = Bianca; i = inoculated; dai=day after inoculation.