Genome-wide Identification of the POD Gene Family and their Expression Profiling in Grapevine (Vitis

The purpose of this study is by extending the analysis of class III peroxidases (PODs) in grapevine and provide further insights into the organ-specific developmental role in transcriptional dynamics and gene duplication analysis of this economically important fruit crop species. Herein, we comprehensively identified 47 PODs in the grapevine genome and are further classified into 7 subgroups based on their phylogenetic analysis. Results of motif composition and gene structure organization analysis revealed that PODs in the same subgroup shared similar conjunction while the protein sequences were highly conserved. Intriguingly, the integrated analysis of chromosomal mapping and gene collinearity analysis proposed that both dispersed and tandem duplication events contributed to the expansion of PODs in grapevine. Also, the gene duplication analysis suggested that most of the genes (20) were dispersed followed by (15) tandem, (9) segmental or whole-genome duplication, and (3) proximal, respectively. The evolutionary analysis of PODs, such as Ka/Ks ratio of the 15 duplicated gene pairs were less than 1.00, indicated that most of the gene pairs exhibiting purifying selection and 7 pairs underwent positive selection with value greater than 1.00. The Gene Ontology Enrichment (GO), Kyoto Encyclopedia of Genes Genomics (KEGG) analysis, and cis-elements prediction also revealed the positive functions of PODs in plant growth and developmental activities, and response to stress stimuli. Further, based on the publically available RNA-sequence data, the expression patterns of PODs in tissue-specific response during several developmental stages revealed diverged expression patterns. Subsequently, 30 genes were selected for RT-PCR validation in response to (NaCl, drought, and ABA), which showed their critical role in grapevine. In conclusion, we predict that these results will lead to novel insights regarding genetic improvement of grapevine.


Background
Peroxidases (EC 1.11.1.X) is a group of well-known large multi-gene family and that are broadly dispersed in living organisms. Based on their structure variations, the peroxidase (PODs) can be characterized into two main groups such as either heme PODs and nonheme PODs [1]. Meanwhile, the heme PODs can be ordered into two more sub-families like animal PODs and non-animal PODs [2]. The non-animal superfamily contains three major sub-distinct classes namely, class I, II, and III [3]. The class III peroxidases (EC 1.11.1.7) are abbreviated in various ways in previous studies (POX, POD, Px, PER, and Prx) and act as plant-specific oxidoreductases [3,4]. In this study, we will use the abbreviation for class III peroxidase as POD. The PODs are involved in a broad range of physiological activities, such as the formation of lignin, cell wall components, defense against pathogenicity or herbivore, and abiotic stress tolerance [5][6][7][8][9].
In recent times, due to the results of transcriptomic data, a large number of PODs have been accompanying numerous biological processes [10][11][12]. However, the direct role of this multi-gene family is still elusive and only a few studies have demonstrated their functional role [13][14][15][16]. For example, Arabidopsis thaliana and Populus trichocarpa PODs (AtPrx72 and PtrPO21), play a significant role in lignification of the leaves [17,18]. The overexpression of the POD genes in A. thaliana ( AtPrx22, AtPrx39, and AtPrx69) improve cold tolerance [19]. Moreover, the cotton GhPOX1 has been studied for a higher production level of reactive oxygen species [20]. Several, POD genes in roots of Zea mays are known to regulate by methyl jasmonate, salicylic acid and pathogen elicitors [21]. Taken together, based on these results, the PODs play an important role in biological, physiological, and in response to stress stimuli, therefore, its comprehensive analysis is necessary to further explore its role in plant growth and development.
In addition, the POD family members have been well-studied and characterized by bioinformatics analysis in several plant species including, Arabidopsis thaliana [22] [26], Zea mays [27], Pyrus bretschneideri [28]. Nevertheless, to date, no previous genome-wide analysis has been carried out of this gene family in grapevine. While, a large number of genes in this family suggested their functional diversity among each individual proteins [11]. Grapevine (Vitis vinifera L) is one of the widely popular and important fruit crops in the world [26]. A common goal of current plant genomics research is to create an expandable platform for global classification and analysis of plant gene family. Hence, it's necessary to provide a foundation for future research. In the meantime, the availability of the grapevine genome (Version 2.1) facilitate the research in grapevine momentously for its genetic studies by improvement in the quality of berry.
In the present study, we performed a wide-ranging bioinformatics analysis of POD gene family and verified their role against various stress responses (i-e., NaCl, drought, and ABA) in grapevine. In total, 47 genes were identified in the grapevine genome and were systematically analyzed by genome-wide approaches. Thus, the study including, physicochemical properties, phylogenetic relationships, chromosomal mapping, collinear correlation, gene duplication events, rate of substitution rates, motif composition and gene structure, promoter sequence analysis, GO and KEGG enrichment analysis, and expression profiling using RNA-seq data and RT-PCR analysis in response to salt, drought and ABA. In general, the results of our study will undoubtedly be helpful for future research on fruits crop species and pay the base for functional characterization of the PODs gene family.

Characterization of POD Gene Family in Grapevine
In this study, a total of 47 POD genes were identified from the grapevine genome and for simplicity, we denominated as VvPOD1-VvIQD47 based on their orthologous position with Arabidopsis thaliana. We also studied some useful information of PODs including, the protein identifier, chromosomal localization, coding sequence (CDS) length (bp), and various physicochemical properties such as, protein length (aa), molecular weight (MW) kDa, isoelectric point (PIs), and grand average of hydropathicity (GRAVY). While, the gene duplication types (i.e., dispersed, tandem, proximal, and segmental or whole-genome duplication) and subcellular localization analysis were also briefly studied for each of POD proteins (Supplementary Table S1 of GRAVY ranged from -0.37 (VvPOD10) to 0.03 (VvPOD25). Intriguingly, the variability was observed in most of the genes for GRAVY, indicating mostly hydrophilic properties and only a few of them (VvPOD46, VvPOD43, VvPOD31, and VvPOD25) are hydrophobic in nature by showing positive values. Additionally, the gene duplication analysis intimated that most of the genes (20) were dispersed followed by tandem (15), segmental or wholegenome duplication (9), and proximal (3), respectively.

Phylogenetic Relationships, Gene Structure Organization of POD Gene Family in Grapevine
To investigate the evolutionary relationships, we used the 47 POD gene grapevine and 73 Arabidopsis thaliana to construct a maximum likelihood approach tree by using MEGA 7.0.
The phylogenetic tree reveals that PODs can be further subcategorized into 7 subgroups ( Figure 1). The results exhibited that there is an uneven distribution of VvPOD genes 6 compared with AtPODs. For instance, we observed that subgroup 7 contains the most number of genes (15 and 17) as compared to other subgroups in grapevine and Arabidopsis. The phylogenetic tree also revealed the relatively close genetic relationships with Arabidopsis.
The 10 conserved motifs ranging from (motif 1-10) of VvPODs were explored by the MEME program. Markedly, motifs 1-4 were most common among the members of PODs, suggesting unique features among subgroups ( Figure 2a). Also, the LOGOS for these motifs were obtained by MEME, the higher number (100) consensus sequences were observed in motif-2 while a less number (50) were recorded in motif 4, and motifs 8-10, respectively (Supplementary Figure S1). The gene structure organization was analyzed based on CDS and untranslated regions (UTRs) by using TBtools. The result reveals that VvPOD members are highly conserved within each other and displayed a similarity among subgroups ( Figure 2b). Further, these findings indicated the structural diversification among VvPOD gene family.

Chromosomal Localization, Gene Collinearity, and Ka/Ks Analysis of POD
To illustrate the chromosomal localization among 47 POD members and the gene collinearity analysis between grapevine and Arabidopsis were drawn with the help of TBtools software. The results for PODs chromosomal localization unveiled the irregular distribution patterns ranging from 1-9 proteins per chromosome except (chr5, chr9, and chr15) across 19 different chromosomes (i.e., Chr01-Chr19) in the grapevine genome. Also, the number of genes on each chromosome were distinct such as the high number of genes (9) were observed on Chr12, followed by Chr1 and Chr18 each with 5 genes, chr6 has 4 genes, while 3 genes were allocated on the Chr7, Chr10 and unknown chromosome (ChrUn), respectively, as described in Figure 3. Thus, among POD members high variation patterns were observed in the grapevine genome. Furthermore, the gene collinearity 7 relationships between V. vinifera ( VvPOD) and Arabidopsis (AtPOD) was also illustrated by using circos plot with the help of TBtools software. As a consequence, high conservation was observed between VvPOD and AtPOD genes ( Figure 3). Table.1. The POD genes in grapevine with outlier Ka/Ks and various types of duplications (i.e., Dispersed, proximal, tandem, and segmental).
The selection pressure among various types of duplications (i.e., dispersed, tandem, proximal, segmental or WGD), also intended by calculating the rates of synonymous substitution (Ks) and non-synonymous substitution (Ka). During evolutionary processes, the genes are typically exposed to various types of selection pressure, such as purifying selection (Ka/Ks<1), positive selection (Ka/Ks>1), and neutral selection (Ka/Ks=1) [29].
Among 47 VvPOD members, we selected 22 pairs (i.e., 10 pair dispersed, 1 pair proximal, 7 pair tandem, and 4 pair segmental or WGD) as presented in Table 1. Results showed that most of the gene pairs having less than 1.00 Ka/Ks ratio suggested purifying selection, thus revealed limited divergence after gene duplications. Though, 7 pairs were observed with higher than 1.00 values, implicating positive selection.

Cis-Regulatory Elements Analysis in Grapevine
The GO enrichment analysis for POD genes was performed to understand their functional regulatory mechanism by using the orthologous pairs of Arabidopsis thaliana. The three common subgroups were observed such as molecular functions (MF), cellular component (CC), and biological process (BP). In the MF processes, "oxidoreductase and catalytic activity" (GO:0016491 and GO:0003824), are highly enriched GO terms. Similarly, for CC processes and BP most of the GO terms are responsive to "cell wall, plasmodesma, symplast, cell-cell junction, plant-type cell wall" (GO:0005618, GO:0009506, GO:0055044, 8 GO:0005911, and GO:0009505), and "response to toxic substance, cellular response to stimulus, oxidation-reduction process, metabolic and cellular process" (GO:0009636, GO:0051716, GO:0055114, GO:0008152, and GO:0009987), and are briefly summarized in Supplementary Table S2. As results, the GO terms for MF, CC, and BP, suggested the crucial role of PODs in various activities of grapevine.
Additionally, the KEGG enrichment analysis indicated the three major pathways among PODs in grapevine such as "Biosynthesis of other secondary metabolites, phenylpropanoid biosynthesis, and metabolism (Supplementary Table S3).
Moreover, the cis-acting elements in the promoter region of POD members were performed by using the PlantCARE database. In brief, most of the genes were largely participating in light regulation with key regulatory elements (GT1-motif, G-Box, GATA-motif, and AE-Box), followed by hormones (CGTCA-motif, TGACG-motif, ABRE, and GARE-motif), stress and other regulatory factors (LTR, ARE, CCAAT-Box, CAT-BOX, o2-site,), and circadian, respectively. Thus, we observed the diversified role of POD members and their indirect involvement in several biotic-abiotic/hormone signaling processes (Supplementary Table   S4).

Expression Profiling of POD Genes in Different Organs and Developmental Stages in Grapevine
The expression profiling of PODs in grapevine was based on 19 various tissue and organs during their developmental stages, the RNA-seq data were retrieved from NCBI (GSE36128) according to the previously reported study [30]. To represent the spatiotemporal expression, a heatmap was generated on FPKM-based (Log 2 ) values of 47 VvPOD genes ( Figure 4). Results revealed that 9 genes (VvPOD1, VvPOD2, VvPOD6, VvPOD10, VvPOD12, VvPOD27, VvPOD32, VvPOD37, and VvPOD46) displayed a striking expression among all tissues and organs, implicating their major participation in tissue-specific 9 response and development. Moreover, the rest of the genes showed either moderate or weak expression abundance in all the selected tissues and organs, speculating their limited response in grapevine.

qRT-PCR Analysis of POD Genes in Response to (NaCl, drought, and ABA)
To investigate the role of VvPOD genes under diverse abiotic stress conditions, we performed qRT-PCR analysis of randomly selected 30 candidate genes and that were subjected to NaCl, drought, and ABA stress treatment. The results directed that all the genes responded variably and showed higher, moderate or low expression level compared to the controls. In response to salt stress, approximately 52% of the total genes showed higher expression level, whereas the rest of the genes showed either moderate or low expression. Interestingly, in the case of ABA and drought stress, about 78% and 72% genes were observed to be down-regulated (Figure 5a). Moreover, the correlation analysis based on Pearson's Correlation Coefficient (PCC) of the relative expression indicated largely a highly positive correlation and some of them were found with inverse correlation ( Figure 5b). Taken together, these results of POD genes based on expression level respond to multiple stresses and might play an important role in the maintenance of plant growth.

Discussion
The PODs multi-gene family are involved in the various biological process by regulating plant growth and developmental processes. While, POD family members have been comprehensively analyzed by genome-wide approaches in several species including,  In this study, we utilized the publically available RNA-seq data of 19 different tissues and organs to validate the tissue-specific response of 47 VvPODs. Among them, 9 genes (VvPOD1, VvPOD2, VvPOD6, VvPOD10, VvPOD12, VvPOD27, VvPOD32, VvPOD37, and VvPOD46) exhibited significantly higher expression level among all the tissue and organs.
The rest of the other genes showed down-regulation with similar tendency throughout grapevine tissues and organs developmental phases.
To better understand the gene functions, expression profiling mainly provides a valuable clue to its regulatory role against specific treatment. Abiotic stresses considerably harbor plant growth and productivity, thus lead to economic losses. The drought severity hampered grapevine productivity and has deleterious impacts on viticulture worldwide [41,42]. Salt stress is one of the severe abiotic stress, while grapevine plants generally modify their physiology to combat salt stress severity [43]. In regulating plant tolerance to several abiotic stress, an eminent hormone, such as ABA plays an essential role. Many researchers have indicated that plant peroxidases are involved in various cellular processes during plant growth and development, and its response to abiotic and biotic stress have been reported over the years [27,44]. For instance, membrane-bounded peroxidase (pmPOX1, pmPOX2a, pmPOX2b, and pmPOX3) in Zea mays roots are regulated by methyl jasmonate, salicylic acid and pathogen elicitors [21]. Similarly, the overexpression of CpPrx01 altered the growth pattern of plants and condensed the level of a well-known hormone, indole-3-acetic acid (IAA) [45]. Thus, in this study, the qRT-PCR validation of 30 POD genes in response to NaCl, drought, and ABA displayed a range of differential expression. Among them, salt stress was found to be highly regulated compared to drought and ABA stress.
To improve and comprehend the possible transcriptional regulation functions of PODs, we also analyzed the GO, KEGG enrichment, and cis-elements in grapevine. Thus, these results inferred their diversified functions and indirect involvement in several biological processes, such as biotic-abiotic/hormone signaling.

Conclusions
In conclusion, we systematically identified a total of 47 POD genes in grapevine and were categorized into 7 subgroups, as supported by phylogenetic analysis. The GO, KEGG, and cis-elements analysis also extended our repositories on the diverse functions of PODs in plant developments during various stress-related activities in grapevine. While, the transcriptional profiling of various organs during several developmental stages and RT-PCR analysis, provide a stand-point of the PODs in improving the plant growth. Thus, the results of our study increase our understanding of POD genes in grapevine and laying the solid foundation for genetic improvement in other fruit crops.

Identification of POD Genes Family in Grapevine
In order to identify the POD genes, we used the 73 reference sequences of Arabidopsis against grapevine (genome version 2.1) with the help of BioEdit tools. In general, we retrieved the sequences of grapevine and Arabidopsis from online sources such as ensembl (https://plants.ensembl.org/index.html) and TAIR (http://www.arabidopsis.org/).
The domain composition was verified by using the NCBI-Conserved Domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART databases (http://smart.embl-heidelberg.de/) [46]. Those sequences with absent of POD domains and sequence with obvious error in length (>100aa length) were eliminated from the study before analysis.

Calculation of Synonymous (Ks) and Non-synonymous (Ka) for Duplicated Genes
The rate of Ka/Ks was carried out for various types of duplicated pairs (i.e., dispersed, proximal, tandem, and segmental) by using MEGA 7.0 [48]. The Ks and Ka ratio was

Gene Structure, Motifs Composition, and Physicochemical Analysis of POD Protein
For the gene structure illustration, we utilized the GFF3 file of the grapevine genome and images were implemented by TBtools software [49]. The motifs analysis of POD protein were performed by the Multiple Em of Motif Elicitation (MEME Suite) version 5.0.5 and then demonstrated by TBtools software. The following parameter was calibrated for this purpose as follows: the maximum number of motifs 10, with a minimum width of 50 and a maximum of 100 and the other parameter was set as default [50]. While, for each gene of POD, the several physicochemical properties (i.e., molecular weight (MW), isoelectronic points (PIs), and GRAVY) were intended by ExPASY PROTPARAM tools (http://web.expasy.org/protparam/). The subcellular localization was further predicted by using the WOLF PSORT online server (https://wolfpsort.hgc.jp/).

Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomics (KEGG) and Cis-Elements Predictions of Class III Peroxidase Gene Family
For the GO and KEGG enrichment analysis, the online panther server (http://pantherdb.org/) and the genome server was implemented (https://www.genome.jp/kegg/pathway.html) and their enriched pathways were further explored by TBtools software [49]. The POD promoter sequences (i.e., selected as 1500bp) was initially imported in Generic File Format (GFF) from the grapevine genome. After that various cis-regulatory elements for each promoter sequence was identified by the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [51].

Chromosomal Mapping and Gene Collinearity Analysis of Grapevine and Arabidopsis
The grapevine genomic database (CRIBI. Available online: http://genomes.cribi.unipd.it/grape/, V2.1), was utilized for the chromosomal locations of POD genes and were mapped based on information available. Similarly, for gene collinearity analysis, the grapevine and Arabidopsis relationship was verified and demonstrated with the help of Circos (TBtools software) program was applied [49].

Plant Material and Methods
In the present study, we used the two-year-old potted grapevine plants (V. vinifera cv. Summer Black), selected from greenhouse condition (25 ± 5 °C) under 16-h light/8-h dark photoperiod and 65% relative humidity (RH) at the Nanjing Agricultural University Nanjing-China. The grapevine plants were subjected to abiotic stress including, NaCl (100 mM), drought (irrigated with 15% (w/v) polyethylene glycol), and ABA (100 μM/100 mL of water) with an interval of 1, 6, 24 h against control (CK), with three biological replicates and the three samples were mixed to make one composite sample while using the 4th unfolded leaf from cv. Summer Black grapevine. At last, all the samples were quickly frozen in liquid nitrogen and stored at −80 °C for further use.

RNA Isolation and Transcriptional Profiling of POD Gene Family in Grapevine
Total RNA was extracted from the grapevine leaves with the help of Trizol (Invitrogen) according to the manufacturer's instructions. Then, the Primer Script RT reagent kit (TAKARA, Dalian, China) was used to convert the RNA into reverse-transcribed into cDNA.
The gene-specific primers (Supplementary Table S5) was designed by the help of Becan Designer 7.9 and their specificity was confirmed using the BLAST tool against the grapevine genome. Also, we performed the RT-PCR by following the previously reported studies [52,53], and the relative fold expression was calculated using the comparative Ctmethod. The expression patterns of all POD gene were analyzed based on a previous study [38,54], and for reference gene in qRT-PCR, we used the housekeeping actin gene   Figure 1 Phylogenetic relationship of POD genes between grapevine and Arabidopsis. The