Revisiting Vitis vinifera Subtilase Gene Family: A Possible Role in Grapevine Resistance against Plasmopara viticola

Subtilisin-like proteases, also known as subtilases, are a very diverse family of serine peptidases present in many organisms. In grapevine, there are hints of the involvement of subtilases in defense mechanisms, but their role is not yet understood. The first characterization of the subtilase gene family was performed in 2014. However, simultaneously, the grapevine genome was re-annotated and several sequences were re-annotated or retrieved. We have performed a re-characterization of this family in grapevine and identified 82 genes coding for 97 putative proteins, as result of alternative splicing. All the subtilases identified present the characteristic S8 peptidase domain and the majority of them also have a pro-domain I9 inhibitor, a protease-associated (PA) domain, and a signal peptide for targeting to the secretory pathway. Phylogenetic studies revealed six subtilase groups denominated VvSBT1 to VvSBT6. As several evidences have highlighted the participation of plant subtilases in response to biotic stimulus, we have investigated subtilase participation in grapevine resistance to Plasmopara viticola, the causative agent of downy mildew. Fourteen grapevine subtilases presenting either high homology to P69C from tomato, SBT3.3 from Arabidopsis thaliana or located near the Resistance to P. viticola (RPV) locus were selected. Expression studies were conducted in the grapevine-P. viticola pathosystem with resistant and susceptible cultivars. Our results may indicate that some of grapevine subtilisins are potentially participating in the defense response against this biotrophic oomycete.


INTRODUCTION
Subtilisin-like proteases (SBTs) are the second largest family of serine peptidases present in archaea, bacteria, eukarya, fungi, and yeast (Siezen et al., 1991). They belong to the S8 family within the SB clan of serine proteases, according to the classification of the peptidase database MEROPS (Rawlings et al., 2014;http://merops.sanger.ac.uk). The majority of plant subtilases are synthesized as an inactive pre-proprotein precursor. Their structure usually presents a signal peptide, a pro-domain (also known as I9 inhibitor domain), a subtilase domain (also known as S8 peptidase domain), and a proteaseassociated domain (PA), although some of them may have only one or even additional domains (Siezen and Leunissen, 1997;Dodson and Wlodawer, 1998;Antão and Malcata, 2005;Siezen et al., 2007;Vartapetian et al., 2011). The presence of a highly conserved catalytic triad within the S8 peptidase domain, composed by aspartate (Asp), histidine (His), and serine (Ser) amino acid residues is characteristic of the subtilase family (Dodson and Wlodawer, 1998). Additionally, certain subtilases may also have a conserved catalytic asparagine (Asn) residue in the same S8 peptidase domain (Siezen and Leunissen, 1997;Dodson and Wlodawer, 1998;Jordá et al., 1999). Furthermore, several plant subtilases also contain a fibronectin (Fn) III-like domain, required for their activation (Rawlings and Salvesen, 2013).
In opposition to mammals on which only nine subtilases have been identified, subtilases from plants are especially abundant, with 63 known genes in Oryza sativa (Tripathi and Sowdhamini, 2006), 56 genes in Arabidopsis thaliana (Rautengarten et al., 2005), 15 genes in Lycopersicon esculentum (Meichtry et al., 1999), 23 genes in the moss Physcomitrella patens, 90 genes in Populus trichocarpa (Schaller et al., 2012), and 82 genes in Solanum tuberosum (Norero et al., 2016). The expansion of the SBT family in plants was accompanied by functional diversification, and novel, plant-specific physiological roles were acquired in the course of evolution. In addition to their contribution to general protein turnover, plant SBTs are involved in the development of seeds and fruits, cell wall modification, processing of peptide growth factors, epidermal development, and pattern formation and in biotic and abiotic stress responses (reviewed in Schaller et al., 2012). In plant-pathogen interactions, the first evidences of subtilase participation were reported by Granell et al. (1987). This subtilase, named pathogenesisrelated protein 69 (P69), was later associated to the response of tomato leaves to Phytophothora infestans Fischer et al., 1989) and characterized as an alkaline proteinase located in the vacuole and intercellular spaces of leaf parenchyma cells (Vera and Conejero, 1988;Vera et al., 1989;Tornero et al., 1997;Jordá et al., 1999Jordá et al., , 2000Meichtry et al., 1999). P69 was also the first plant subtilase for which two protein substrates were identified, systemin (Schaller and Ryan, 1994) and the leucine-rich repeat protein (LRP; Tornero et al., 1996), although the consequences of these substrates processing events for plant pathogen interaction still remain unknown. More recently, Ramírez and co-workers have identified a SBT3.3 gene from A. thaliana as encoding a serine protease homolog to the P69C subtilase from tomato and associated its function in immune priming responses (Ramírez et al., 2013). Also, in S. tuberosum, expression profile analysis of detached potato leaves after P. infestans infection or after BABA or BTH treatment highlighted an expression increase of several subtilases genes (Norero et al., 2016). Moreover, the subtilase St SBTc-3 has been found as a major protein in apoplast of detached potato leaves after P. infestans infection (Fernández et al., 2012) and it was shown that this subtilase evidences DEVDasa activity and is related to programmed cell death functions (Fernández et al., 2015).
In grapevine, the first clues highlighting subtilase participation defense mechanisms were reported by Figueiredo et al. (2008), when comparing resistant and susceptible genotypes prior and post-inoculation with Plasmopara viticola. A subtilisinlike protease (XM_010660203.1), identified as a cucumisin was constitutively expressed in resistant genotype and increased its expression after P. viticola inoculation (Figueiredo et al., 2008(Figueiredo et al., , 2012Monteiro et al., 2013). Also in this pathosystem, it was shown that after treatment with serine protease inhibitors, plants became more sensitive to P. viticola (Gindro et al., 2012). It was hypothesized that some components of P. viticola secretome could inhibit the endogenous subtilases of susceptible varieties, thereby inhibiting the plant's normal defense reaction, while resistant or immune varieties may possess endogenous subtilases that are not recognized due to slight structural modifications of the protein patterns of these cultivars. In this case, plant defense mechanisms would continue to operate, with fatal consequences for the pathogen and restricting its development (Gindro et al., 2012).
In 2014, a first attempt to characterize the grapevine subtilase family was made by Cao and co-workers where the subtilase sequences were identified based on the presence of the PA domain (Cao et al., 2014). In parallel, a re-annotation of the grapevine genome was conducted (Vitulo et al., 2014) and several subtilase sequences were either re-annotated or retrieved. The aim of this work was to identify subtilisin-like proteases in the grapevine genome and characterize them based on phylogenetic analyses, gene and protein primary structure. Additionally, expression analysis of selected subtilase genes was conducted to identify subtilases potentially involved in downy mildew resistance.
Gene exon/intron structure information was collected from grapevine genome annotation at NCBI.

Phylogenetic Analysis
Two phylogenetic analyses were carried out. The first one with all the identified grapevine subtilase amino acid sequences (97) and a second one combining these amino acid sequences with 56 SBTs amino acid sequences from A. thaliana (Rautengarten et al., 2005) and 14 from Solanum lycopersicum described in Meichtry et al. (1999).

Selection of Grapevine Subtilase Sequences Putatively Involved in Pathogen Resistance
Previous studies in plants associated some subtilases with the defense response to pathogen attack, like the subtilase SBT3.3 in A. thaliana (Ramírez et al., 2013), the P69 in S. lycopersicum (Tornero et al., 1996;Jordá et al., 1999) and the cucumisin in grapevine (Figueiredo et al., 2008(Figueiredo et al., , 2012. The subtilase genes from V. vinifera were blasted against the A. thaliana genome (TAIR database, https://www.arabidopsis.org/) and the tomato genome (SolGenomics database, https://solgenomics.net/) to retrieve the grapevine sequences presenting higher sequence similarity to A. thaliana SBT3.3 and tomato P69 genes. SolGenomics results were corroborated in NCBI BLAST tool, restricting to S. lycopersicum organism, and was assumed the NCBI accession for further studies. Moreover, subtilase sequences with a chromosomal location near the RPV locus on grapevine genome were also selected for further studies. Multiple alignment of the grapevine subtilases selected as putatively involved in plant resistance was performed in DNASTAR software (version 13, Burland, 1999; http://www.dnastar.com/).
The protein interaction network of the selected subtilases was obtained (STRING, version 10.0, Szklarczyk et al., 2014; http://string-db.org/). The gene accessions for all proteins that interact with the selected grapevine subtilases were queried at the NCBI database. The gene ontology (GO) terms for all the interacting proteins were also obtained with the Blast2GO tool (version 3.3, Conesa et al., 2005; https://www.blast2go.com/).

Experimental Design for Expression Analysis: Plant Material and Inoculation Experiments
Two expression analysis experiments were conducted: (1) using two grapevine cultivars (resistant and susceptible) inoculated with P. viticola to access subtilase expression during inoculation; (2) using several grapevine accessions with different degrees of resistance toward P. viticola to access subtilase constitutive expression.
For the first analysis, two Vitis vinifera cultivars were selected to access subtilase expression during interaction with P. viticola. The cultivar Regent, bread by multiple introgressions from resistant wild genotypes (Welter et al., 2007), presenting a high degree of resistance to downy and powdery mildews (Anonymous, 2000), and Trincadeira, a Portuguese traditional grapevine cultivar widely used for quality wine production and highly sensitive to P. viticola (Figueiredo et al., 2008). Both cultivars were propagated under identical greenhouse conditions, briefly grapevine wood cuttings were grown in 12 cm diameter pots in Fruhstorfer Erde (soil) Type P for 10 weeks, under natural day/night rhythm with temperatures ranging between 5 and 28 • C, according to Figueiredo et al. (2012). For plant inoculation, P. viticola sporangia were collected after an overnight incubation of symptomatic leaves from greenhouse infected plants in a moist chamber at room temperature. Sporangia were carefully recovered by brushing, dried, stored at −25 • C and checked for their vitality by microscopy as in Kortekamp et al. (2008). A suspension containing 10 4 sporangia ml −1 was used to spray the abaxial leaf surface in order to challenge the plants. Mock inoculations with water were also made. After inoculation, plants were kept in a moist chamber (100% humidity) for 8 h and then under greenhouse conditions at 25 • C during the inoculation time course. The third to fifth fully expanded leaves beneath the shoot apex were harvested at 6, 12, and 24 hpi, immediately frozen in liquid nitrogen and stored at −80 • C. For each genotype and condition (inoculated and mock inoculated), three independent biological replicates were collected, being each biological replicate a pool of three leaves from three different plants.
For the second analysis, to access if subtilases are constitutively expressed, young leaves from several resistant (V. labrusca, V. rupestris, V. rotundifolia, V. riparia, and V. candicans) and tolerant (V. sylvestris) Vitis species and V. vinifera cultivars [Regent (resistant) and Trincadeira (susceptible)] were harvested from five different plants (per biological replicate), at the Portuguese Grapevine Germplasm Bank at INIA-Estação Vitivinícola Nacional (Dois Portos), (Supplementary Data 1). Leaves were immediately frozen in liquid nitrogen and stored at −80 • C. Three biological replicates were collected from each accession.

RNA Extraction and cDNA Synthesis
Total RNA was isolated from frozen leaves with the Spectrum TM Plant Total RNA Kit (Sigma-Aldrich, USA), according to manufacturer's instructions. Residual genomic DNA was digested with DNase I (On-Column DNase I Digestion Set, Sigma-Aldrich, USA). RNA purity and concentration were measured at 260/280 nm using a spectrophotometer (NanoDrop-1000, Thermo Scientific) while RNA integrity was verified by agarose gel electrophoresis (1.2% agarose in TBE buffer). Genomic DNA (gDNA) contamination was checked by qPCR analysis of a target on the crude RNA (Vandesompele et al., 2002). Complementary DNA (cDNA) was synthesized from 2.5 µg of total RNA using RevertAid R H Minus Reverse Transcriptase (Fermentas, Ontario, Canada) anchored with Oligo(dT) 23 primer (Fermentas, Ontario, Canada), according to manufacturer's instructions.

Quantitative Real Time PCR
Quantitative RT-PCR (qPCR) experiments were carried out using Maxima TM SYBR Green qPCR Master Mix (2×) kit (Fermentas, Ontario, Canada) in a StepOne TM Real-Time PCR system (Applied Biosystems, Sourceforge, USA). A final concentration of 2.5 mM MgCl 2 and 0.2 µM of each primer were used in 25 µL volume reactions, together with 4 µL of cDNA as template. Primer sequences and reaction details are provided in Supplementary Data 2.
Thermal cycling for all genes started with a denaturation step at 95 • C for 10 min followed by 40 cycles of denaturation at 95 • C for 15 s and annealing (Supplementary Data 2) for 30 s. Each set of reactions included a control without cDNA template. Dissociation curves were used to analyse non-specific PCR products. Three biological replicates and two technical replicates were used for each sample. Gene expression (fold change) was calculated by the Hellemans et al. (2007). The reference genes used for the normalization were the previously described in Monteiro et al. (2013). Statistical significance (p < 0.05) of gene expression was determined by the Mann-Whitney U test using IBM R SPSS R Statistics version 23.0 software (SPSS Inc., USA).

Identification and Characterization of Subtilisin-Like Serine Protease Genes in Grapevine
The first characterization of the grapevine subtilase family was made by Cao and co-workers in 2014, where 80 subtilase genes were identified (Cao et al., 2014). Subtilase search was restricted to the subtilase conserved PA domain, although subtilases are usually characterized by three conserved domains (PA, S8 peptidase and I9 inhibitor). In parallel, the grapevine genome was re-annotated (Vitulo et al., 2014), nine of the previously identified genes were completely removed from the databases and eight were re-annotated. Thus, we have performed a new characterization of this family in grapevine using the subtilase PA, S8 peptidase and I9 inhibitor domains as query in the new grapevine genome annotation version. Eighty-two V. vinifera subtilase genes were identified, from which it is predicted to obtain 97 subtilase proteins (Supplementary Data 3). This search resulted in the introduction of 17 new subtilase genes and the re-annotation of 8 from the subtilase genes previously identified (Cao et al., 2014). The number of genes members found in V. vinifera subtilase family is similar to the one from P. trichocarpa (90 subtilase genes; Schaller et al., 2012) and potato (82 subtilase genes; Norero et al., 2016) and higher than those reported in other plant species, like Arabidopsis or tomato, which were detected 56 and 15 subtilase genes, respectively (Meichtry et al., 1999;Rautengarten et al., 2005).
The 82 identified genes identified were mapped in V. vinifera chromosomes. These genes were unevenly distributed among 15 of the 19 grapevine chromosomes (Figure 1). No subtilase genes were detected on chromosomes 1, 5, 14, and 17, and the specific location of 8 of the 82 subtilase genes is still unknown. The majority of the subtilase genes were located on chromosomes 6, 13, and 16, with 9 genes in chromosome 6 and 10 in chromosomes 13 and 16 (Figure 1).

Gene Structure Analysis
Grapevine STB genes were checked for exon-intron structure. Details of the exon-intron structures are shown in Supplementary Data 3. The number of introns varied between 0 and 18, around 24% of the grapevine subtilase genes are intronless and 5% present a high number of introns (17-18 introns). Intronless subtilase genes have been reported in Arabidopsis and potato (Rautengarten et al., 2005;Norero et al., 2016) being the highest number of introless genes reported in potato (63% of the StSBT genes; Norero et al., 2016). Intronless genes can serve as beacons in analyses of gene function and evolution, they have been found in large gene families and related to gene duplications, inheritance from ancient prokaryotes, retroposition or other mechanisms (Yan et al., 2016). In general, most closely related members in the same group shared a similar exon-intron structure.

Protein Structure and Domain Analysis
Both molecular weight (Mw) and isoelectric point (pI) were predicted for the grapevine subtilase proteins (Supplementary Data 3). V. vinifera subtilase proteins have a wide range of molecular weights, between 19 and 164 kDa, slightly higher than the already described for other plant serine proteases (between 19 and 110 kDa; Antão and Malcata, 2005). The majority (60%) present a theoretical molecular weight between 80 and 90 kDa, and only 26% between 60 and 80 kDa as previously described for other plant serine proteases (Antão and Malcata, 2005). The remaining 14% presented a molecular weight lower than 60 kDa or higher than 90 kDa. Grapevine subtilases have a theoretical pI between 4.69 and 9.57. This pI range is comparable to other subtilase proteins, like, for example the STB3.3 (AT1G32960.1) from A. thaliana and P69C (CAA76726) from tomato with a theoretical pI of 6.27 and 5.27, respectively (predicted from protein sequence with Compute pI/Mw tool from ExPASy, http://web.expasy.org/compute_pi/) or the CpSUB1 from papaya with a pI of 8.97 (Othman and Nuraziyan, 2010).
Subtilases are characterized by a multidomain structure comprising a signal peptide, a propeptide, a protease domain (S8) and a protease-associated (PA) domain (Siezen et al., 2007). The protease domain (S8 domain) that defines the subfamily S8A (PF00082), includes the catalytic triad and a proteaseassociated domain (PA) (PF02225) which is an insertion of about 120-160 amino acids long between the His and Ser active site residues that cause a displacement of the reactive Ser from the catalytic triad to the C-terminal (Siezen and Leunissen, 1997). This domain has been implicated in protein-protein interactions or substrate specificity (Siezen and Leunissen, 1997 Data 3). This I9 inhibitor domain (PF05922), also known as pro-domain N terminus or propeptide, is involved after removal, on the pro-enzyme activation, working as a molecular chaperone in the folding of the mature peptidase. Thus, the I9 inhibitor prevents the access of the substrate to the active site and activates the peptidase when it is removed either by autocatalytic cleavage or by interaction with a secondary peptidase (Siezen, 1996). The PA domain was detected in all of the subtilase sequences, but only in 46 sequences the Evalue was considered significant (E ≤ −5). On the other 39 sequences the presence of the PA domain presents low E values or is listed as unintegrated signature. However, the shift on the reactive Ser from the catalytic triad to the C-terminal is present in all of the grapevine subtilase sequences, which may suggest sequence divergence for this domain.
Not all the grapevine subtilases exhibited the three domains simultaneously. Despite being conserved in plant subtilases, the simultaneous presence of I9 inhibitor, S8 peptidase, and PA domain is not obligatory requisite. Moreover, it is yet to be confirmed if the simultaneous presence or not of this set of domains has some effects in subtilase functions. An example of the non-simultaneous existence of the three conserved domains is the P69C subtilase from tomato (Tornero et al., 1997).
The presence of subtilases with domain repeats can be a result of the evolution and a way to improve the subtilase features and its functions. Gene duplication and mutation processes in biological evolution have been largely recognized since the 1930s (Bridges, 1936;Brown and Doolittle, 1995;Zhang, 2003). Gene duplication may result in domain repeats in protein structure. These repeats have a rich variety of functional properties involving protein-protein interactions as well as binding to other molecules like DNA or RNA. Furthermore, long tandem of repeats can play an important role in the folding of three dimensional structures of multi-domain proteins. Structural studies in proteomics have shown that the abundance of domain repeats in organisms of higher complexity is highly correlated with domain families involved in complexassembly, cell-adhesion, and signaling processes (Han et al., 2007).
Six VvSBT presented an additional domain, the fibronectin (Fn) III-like domain (PF06280), (Supplementary Data 3). This domain of unknown function is required by some plant subtilases for their activity (Rawlings and Salvesen, 2013). The SBT3.3 subtilase structure from A. thaliana showed the presence of the three conserved domains and also a fibronectin (Fn) III-like domain (Rose et al., 2010).

Subcellular Location Prediction
Predictions of the subcellular location of a gene product can provide additional information for its functional involvement. Different subcellular locations of plant subtilases have been found to correlate with their different physiological functions (Rautengarten et al., 2005;Cao et al., 2014). For example, the CpSUB1 subtilase from papaya is secreted to extracellular space, where it plays a role in the early stage of fruit development and ripening by degrading cell wall matrix (Othman and Nuraziyan, 2010). Rice subtilase RSP1 is only present in the reproductive organ and absent in leaves, roots, embryos, or rice panicles (Yoshida and Kuboyama, 2001). This suggests that the role for each plant subtilase is related to its location event in spite of analogous structural features (Othman and Nuraziyan, 2010). Sixty VvSBTs possess signal sequences for targeting to the secretory pathway (S), 26 subtilases do not contain any known targeting motif, 8 family members are predicted to be targeted to mitochondria (M) and 3 to the chloroplast (C), (Supplementary Data 4). Hence, grapevine subtilases may have a diversity of functions, most of them probably with roles in the extracellular space or matrix.

Phylogenetic Analysis of Grape Subtilases
A phylogenetic analysis of the 97 grapevine subtilase proteins was carried out and the consensus phylogeny obtained is shown in Figure 2. Based on the phylogenetic relationships of the grapevine subtilases proteins found, an outgroup was identified and 6 SBT groups were established and named VvSBT1 to VvSBT6 (Figure 2). VvSBT1 comprise 23 subtilases and include all the subtilases annotated as cucumisin (degradative subtilases from Cucumis melo). VvSBT1 coding genes are unevenly distributed among grapevine chromosomes and present a structure with variable number of introns (Supplementary Data 3). VvSBT2 is the smallest group containing 4 subtilase proteins annotated as xylem serine proteins. All of the coding genes present a 10 intron structure and are distributed between chromosome 8 and 15 (Supplementary Data 3). The VvSBT3 group comprises 14 subtilases proteins, most of them presenting similarity with A. thaliana SBT3.3/SBT3.5 proteins. All of the coding genes present more than 9 introns and an uneven chromosomal location (Supplementary Data 3). VvSBT4 includes 13 subtilase proteins, most of them showing similarity to A. thaliana SBT5.3/SBT5.4. Within this group subtilase coding genes present between 5-10 introns and for most of the genes the chromosomal location is unknown (Supplementary Data 3). VvSBT5 is the largest group including 34 subtilase proteins that are annotated as subtilisin-like proteins. Within this group, most of the coding genes are intronless or present 1-2 intron. VvSBT6 is comprised by 6 subtilase proteins all containing an additional fibronectin III-like domain, the coding genes are all located in chromosome 11 and present 13 introns (Supplementary Data 3). Rautengarten and co-workers have equally performed a phylogenetic analysis of the predicted 56 A. thaliana (AtSBT) subtilase sequences that showed a division of the subtilases into 6 groups (Rautengarten et al., 2005). In potato, 5 subtilase groups were considered (Norero et al., 2016).

Phylogenetic Analysis of Grapevine, Tomato, and Arabidopsis Subtilases
Biological functions of plant subtilases remain largely unknown. The phylogenetic analysis of 168 amino acid sequences including grapevine, tomato and Arabidopsis subtilases evidenced 8 clusters named from I to VIII (Figure 3).
Cluster II includes the grapevine VvSBT6 group, the three subtilase sequences considered as an outgroup (Figures 2, 3) and the AtSBT2 group (which included the ALE1 protease necessary for cuticle formation and epidermal differentiation during embryo development in A. thaliana (Tanaka et al., 2001).
Cluster III comprises the AtSBT3 group including the AtSBT3.3 recently described as being involved in immune priming events in Arabidopsis (Ramírez et al., 2013) and several proteins from the VvSBT3 group (Figure 2).
Cluster V includes two grapevine subtilases from SBT3 group, both named CO(2)-response secreted protease and two AtSBT subtilases, AtSBT5.1 and 5.2, suggesting that these subtilases may share common functions.
Cluster VI groups the VvSBT2 group and the AtSBT4.14 and SBT4.15 involved in the xylem differentiation (Zhao et al., 2000).
Cluster VII is defined by the remaining members of the AtSBT4 group and seven grapevine subtilases from the VvSBT1 group, all presenting the coding genes located in chromosome 6 (Supplementary Data 3).
Clusters VIII includes all of the remaining sequences from the VvSBT1 group, being located in chromosome 16 and 13, respectively. These VvSBT could have evolved separated from Arabidopsis subtilases.
The 14 selected grapevine subtilases potentially linked to V. vinifera immunity were further analyzed. A prediction of glycosylated sites was performed in the selected grapevine proteins ( Table 3) as it has been shown that glycosylated plant subtilases are secreted to plant extracellular matrix (ECM; Bykova et al., 2006;Cedzich et al., 2009). Since the ECM is where the first host-pathogen interaction, recognition and signaling events take place (Dixon and Lamb, 1990), the accumulation of subtilases in plant ECM may account for an important role during pathogenesis. The most important protein glycosylation form is N-linked, formed by the covalent attachment of asparaginelinked carbohydrates to the protein (Gupta and Brunak, 2002;Bykova et al., 2006). Protein N-glycosylation was previously described in subtilases P69B from tomato (Bykova et al., 2006). From the 14 protein sequences analyzed, only two may not contain a signal peptide (XM_002275345.2, XM_010659200.1), and thus may not be glycosylated in vivo, even though they contain potential motifs. The remaining 12 proteins seem to contain a signal peptide and N-glycosylation was predicted in several Asp residues ( Table 3).
The protein-protein interaction network for the selected subtilases putatively involved in grapevine immunity was performed. By understanding the protein environment where these proteins are likely involved, it is feasible to obtain relevant information about their function and the biological processes. For this analysis, the STRING database was used (Szklarczyk et al., 2014). The top 50 proteins that interact with the 14 grapevine subtilases were analyzed individually in UniProt to access the biological processes to which they are associated. Five of these proteins were predicted to interact with all of the selected grapevine subtilases and are involved in biological processes associated to defense responses, namely fatty acid beta-oxidation, protein kinase activity, ER-associated ubiquitin-dependent protein catabolic process, defense response and protein serine/threonine kinase activity (Supplementary Data 5). Lipid peroxidation and lipidassociated signaling have been recently associated to grapevine resistance to P. viticola (Figueiredo et al., 2015;Guerreiro et al., 2016). Also protein kinases are known to regulate the majority of cellular pathways, especially those involved in signal transduction (Dhanasekaran and Premkumar Reddy, 1998). This proteinprotein interaction network result reinforces the hypothesis that these 14 subtilases may have some involvement in the grapevine immunity.

Expression Analysis
Expression profiles of the selected grapevine subtilases were first analyzed in two V. vinifera cultivars, Regent and Trincadeira (resistant and susceptible to P. viticola, respectively) at 6, 12 and 24 h after inoculation (hpi) with P. viticola. Early inoculation time points were chosen to access the signaling events during pathogen recognition: between 6 and 12 hpi stomatal penetration and development of stomatal vesicles with primary hyphae occurs and at 24 hpi elongated hyphae invade the intercellular space of the mesophyll progressing to the branching stage in susceptible plants and stopping the development in resistant plants (Kortekamp and Zyprian, 2003;Unger et al., 2007).
Of the 14 grapevine subtilases analyzed, two presented no amplification (XM_002275393.2: sequence similarity with P69C and XM_002277863.3: located in rpv9) in both cultivars and were retrieve from our study. Three subtilases that present either sequence homology with AtSBT3.3 (XM_003634104.2 and XM_002273159.3) or are located in the rpv loci (XM_010659200.1) presented the same expression pattern during inoculation time-course in the resistant cultivar Regent, being up-regulated at 6 hpi, decreasing their expression at 12 hpi and increasing again at 24 hpi (Figure 4). In the susceptible genotype, these subtilases were more expressed at 24 hpi (Figure 4). Also, the subtilase XM_003634105.2 (presenting sequence homology with AtSBT3.3) showed an up-regulation at the early time-points analyzed (6 and 12 hpi) in the resistance cultivar Regent, being down-regulated at 24 hpi. In the susceptible cultivar Trincadeira, this gene was down-regulated at 6 hpi and increased its expression from 12 hpi (Figure 4). Accordingly, in Vitis pseudoreticulata leaves infected with the biotrophic ascomycete Erysiphe necator (Schw.) Burr., presented an up-regulation of all these subtilase genes after infection (Weng et al., 2014).
The grapevine subtilase gene XM_010649370.1, showing homology with AtSBT3.3, exhibited an increase of expression  Comparison between chromosomal location of grapevine subtilase genes and chromosomal location of "Resistance to Plasmopara viticola-RPV" loci. Asterisk (*) indicate subtilases selected for further studies.
at both 6 and 24 hpi, being down-regulated at 12 hpi, in the resistant cultivar, while in the susceptible cultivar this gene is over-expressed during all the inoculation time-course (Figure 4). The expression of this subtilase has been previously analyzed in different grapevine tissues and abiotic stimuli (Cao et al., 2014), exhibiting a constitutive high-level of expression in roots, leaves, and stem and the expression was supressed in abiotic (salt, heat, cold, and drought) stress conditions (Cao et al., 2014), which may suggest that it could have a participation in response to biotic stimulus instead. Recently, a Gossypium babardense subtilase gene GbSBT1, that show a high sequence homology with the grapevine subtilase XM_010649370.1, was associated with the defense response to extracellular stimulations like Verticillium dahlia infection, the cause of the Verticillium wilt disease. The corresponding GbSBT1 protein showed to be mainly localized at the cell membrane and moves to the cytoplasm following treatment with jasmonic acid and ethylene, which supports the hypothesis of some grapevine subtilases are located near the place where occurs the plant-pathogen interaction. Moreover, Duan and co-workers, observed a reduction of the tolerance of a cotton resistant genotype, when the GbSBT1 gene was silenced, and an activation of the expression levels of defense-related genes (Duan et al., 2016). Two subtilases presenting sequence homology with AtSBT3.3 (XM_010663620.1) and located at the rpv9 (XM_002284065.3) showed a similar expression pattern, being up-regulated at the later time-point analyzed (24 hpi) in the resistant cultivar Regent, however in the susceptible cultivar Trincadeira both genes have an earlier and higher expression starting at 12 hpi. The expression of the subtilase XM_002284065.3 was also analyzed by Cao and co-workers and an increase of expression during abiotic stress was shown (Cao et al., 2014). Both results suggest an involvement of this subtilase in response to abiotic and biotic environmental stimulus.
The expression of XM_002278414.3 (presenting homology with AtSBT3.3) and XM_002275435.2 (presenting homology with tomato P69C) was either not altered or down-regulated during inoculation time-course in both cultivars (Figure 4). Accordingly this subtilase presented low when submitted to abiotic stress conditions (heat and drought; Cao et al., 2014).   Two subtilases showing sequence homology with tomato P69C (XM_002275345.2 and XM_002275374.2) were upregulated at 12 hpi and down-regulated at the other timepoints (6 and 24 hpi), in the resistant cultivar Regent, and not amplified in the susceptible cultivar Trincadeira. In S. tuberosum leaves, was observed the up-regulation of the expression of three genes, homologous to these two grapevine subtilase genes, after inoculation with Phytophthora infestans or elicitation with DLβ-aminobutiric acid (BABA), a resistance gene inductor. This may suggest that these two genes could be pathogen induced and associated to the defense responses of resistant genotypes (Norero et al., 2016).
When comparing both incompatible (Regent) and compatible (Trincadeira) interactions it is clear that the increase in subtilase expression in Trincadeira presents a delay when compared to the resistant cultivar in which several subtilases are highly expressed 6 hpi (Figure 4). An early increase of expression of some subtilases in the incompatible interaction may be related to the successful establishment of a defense strategy against the invading pathogen.
None of the subtilases analyzed exhibited an up-regulated constitutive expression in comparison to Trincadeira (Figure 5). FIGURE 4 | Heatmap of the 14 grapevine subtilase expression in V. vinifera cv Regent and V. vinifera cv Trincadeira at 6, 12, and 24 hpi with P. viticola. Each column indicates a time-point (6, 12, and 24 hpi) and each row represents a subtilase gene in the resistant grapevine genotype (Regent) or in the susceptible grapevine genotype (Trincadeira) and was colored according to the log2 ratio of expression. Green indicates lower expression, red indicates higher expression, black indicates no expression (see the color scale) and indicates no amplification. Asterisks (*) represent significant difference (p ≤ 0.05) between target and control samples (Mann-Whitney U-test; SPSS Inc., USA, V20).

CONCLUSIONS
This work presents the identification and characterization of the grapevine subtilase family, comprising 82 genes (20 of each are intronless), after the grapevine genome reannotation (Vitulo et al., 2014). These genes present an uneven distribution along 15 of the 19 grapevine chromosomes and encode 97 putative VvSBT proteins, due to alternative splicing. Phylogenetical analysis allowed the characterization of six groups (VvSBT1-VvSBT6) based on amino acid similarity.
Our results suggest that grapevine subtilases may exert several functions and that several grapevine subtilases may be potentially involved in pathogen defense, particularly to P. viticola. We have shown that several grapevine subtilases present sequence similarity with Arabidopsis SBT3.3 and tomato P69, some are located in the P. viticola resistance associated locis (Rpv) or are up-regulated during P. viticola infection. Also, the majority of subtilases are predicted to be located in the extracellular space which reinforces their putative role in the defense mechanisms against pathogens. XM_010660203.1, the cucumisin previously identified as possibly involved in the grapevine defense mechanisms (Figueiredo et al., 2008(Figueiredo et al., , 2012Monteiro et al., 2013), presented the higher increase of expression after inoculation and it is constitutively expressed in several resistant grapevine genotypes, thus it may be considered a strong resistance-associated candidate. More studies must be conducted to define subtilase functions and their role on plant-pathogen interactions particularly in the grapevine resistance against P. viticola.

AUTHOR CONTRIBUTIONS
AF conceived the study and planned the experiment. JF and MM performed the experiments; JF has done the bioinformatic analysis. GC and OP performed the phylogenetical analysis. AF, MS, RM, and JF performed data analysis. AF, JF, and MS wrote the manuscript. All authors have read and approved the manuscript.