Genome-Wide Identification, Characterization, and Expression Analysis of the Grapevine Superoxide Dismutase (SOD) Family

Abstract

Superoxide dismutase (SOD) is an essential enzyme of the plant antioxidant system that responds to oxidative damage caused by adverse conditions. However, little is known about the SOD gene family in Vitis vinifera (Vv). In the present study, ten SOD genes, including 6 copper/zinc SODs, 2 iron SODs, and 2 manganese SODs, were identified in the grapevine genome where they were unevenly distributed on 12 chromosomes. Ten VvSOD genes were divided into three main groups based on phylogenetic analysis, subcellular localization, and the distribution of conserved protein motifs. Additionally, many cis-elements related to different stresses were found in the promoters of the 10 VvSOD genes. Syntenic analysis revealed that VvMSD1 and VvMSD2 were derived from segmental duplication, and VvCSD4 and VvCSD5 belong to a pair of tandemly duplicated genes. Gene expression analysis based on microarray data showed that the 10 VvSOD genes were expressed in all the tested tissues. Interestingly, the segmentally duplicated gene pair (VvMSD1 and VvMSD2) exhibited differential expression patterns in various organs. In contrast, the tandemly duplicated gene pair (VvCSD4 and VvCSD5) displayed similar expression patterns in the tested organs. Our results provide a basis for further functional research on the SOD gene family in grapevine.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Grapevine (Vitis vinifera) is one of the most cultivated and economically valuable fruit crops in the world. However, abiotic stresses, particularly drought and salt stress, threaten the growth of grapevine globally, thereby affecting fruit yield and quality [1–3]. The most common result of such stress is the generation of toxic reactive oxygen species (ROS). Membrane damage, protein oxidation, DNA lesions, and even irreparable metabolic dysfunction and cell death can occur as a result of excessive ROS (such as superoxide anions, hydroxyl radicals, hydrogen peroxide, and singlet oxygen) [4, 5].

To cope with ROS toxicity, plants have developed efficient and sophisticated antioxidative response systems, including the synthesis of low relative molecular mass antioxidant molecules (e.g., l-ascorbic acid and glutathione) and various enzymes. Among these enzymatic components, the superoxide dismutases (SODs) constitute the first line of defense against ROS by catalyzing the dismutation of the superoxide O₂^-to O₂ and H₂O₂ [6, 7]. In plants, SODs have been detected in the roots, leaves, fruits, and seeds where they play important roles in protecting cells from oxidative damage [8].

Multiple SOD isozymes exist in plants and are classified into three types based on their metal cofactor, protein folds, and subcellular distribution: copper/zinc- (Cu/Zn-) SOD, manganese- (Mn-) SOD, and iron- (Fe-) SOD, and these SODs are located in different compartments of the cell [7, 9, 10]. Cu/Zn-SODs are mainly distributed in the cytosol, chloroplasts, peroxisomes, and/or the extracellular space, while Fe-SODs mainly occur in the chloroplasts, as well as in peroxisomes and mitochondria, and Mn-SODs are mainly localized not only in the mitochondria but also in different types of peroxisomes [11]. In addition, a new type of SOD, nickel (Ni)SOD, was first discovered, cloned, and characterized in Streptomyces [12]. However, no evidence for NiSOD has been found in plants [13].

In recent years, some studies have reported that SODs contribute to the response to various environmental stimuli in plants, such as cold, drought, salinity, auxin, and ethylene [14–16]. Owing to their crucial roles in the antioxidant system, a considerable number of SOD genes have been cloned from various monocot and dicot plants [17–19]. With the development of high-throughput sequencing technologies, the genome-wide identification of SOD gene families has been performed in many plant species, including Arabidopsis thaliana [20], Dimocarpus longan [17], Sorghum bicolor [21], Populus trichocarpa [9], Musa acuminate [22], Gossypium raimondii, Gossypium arboreum [23], and Gossypium hirsutum [24].

In the present study, we identified the SOD gene family in grapevine on the genome-wide scale. Genomic organization, gene structure, motif composition, subcellular localization, syntenic analysis, and phylogenetic relationships were analyzed using bioinformatics. Then, the gene ontology and putative promoters of the grapevine SODs were also investigated, and the cis-elements involved in stress responses were analyzed. Putative transcription factors (TFs) which may regulate the VvSOD genes were predicted. Furthermore, we studied the expression patterns of the grapevine SOD gene family under abiotic stress (salt, drought, cold, and heat) using microarray data and a real-time quantitative- (qRT-) PCR detection system. Finally, coexpression networks of VvSOD genes and TF genes under the four abiotic stresses were generated based on the qRT-PCR data. This was the first comprehensive study of the SOD gene family in grapevine and may provide valuable information for understanding the classification, evolution, and putative functions of the grapevine SOD family on the genome-wide scale.

2. Materials and Methods

We downloaded the most recent version (V 2.1) of a 12X assembly of the grapevine (V. vinifera) genome from CRIBI (http://genomes.cribi.unipd.it/grapevine/) for identifying SOD genes. To identify members of the SOD gene family in grapevine, the full-length protein sequences from A. thaliana (TAIR locus number: AT1G08830.1, AT2G28190, AT5G18100, AT4G25100, AT5G51100, AT5G23310, AT3G56350, and AT3G10920) were used as BLASTp queries against all the grapevine protein sequences with the threshold expectation value set to 1.0. All the hits were further submitted to Pfam analysis (http://pfam.xfam.org/) to verify the presence of the SOD domain. Physicochemical characteristics of the SOD proteins were calculated using the ProtParam tool (http://web.expasy.org/protparam/), including the number of amino acids, molecular weight, theoretical isoelectric point (pI), aliphatic index, and grand average of hydropathicity (GRAVY) [25].

ProtComp9.0 server (http://linux1.softberry.com/) was used to predict the subcellular localizations of the SOD proteins [22]. Conserved protein motifs were predicted by MEME Suite (http://meme-suite.org/tools/meme) with the default settings, except that the number of motifs was set to 8, and the minimum and maximum motif widths were set to 20 and 150 amino acids [26]. The web-based Gene Structure Display Server (GSDS; http://gsds.cbi.pku.edu.cn/index.php) program was utilized to illustrate exon and intron organization for the grapevine SOD genes.

2.1. Phylogenetic Analysis

To investigate the phylogenetic relationships of the SODs among grapevine and A. thaliana, multiple SOD protein sequences were aligned, and an unrooted phylogenetic tree was constructed in MEGA 6.06 [27]. The phylogenetic tree was constructed using the neighbor-joining (NJ) method and Jones-Taylor-Thornton (JTT) model. In the phylogenetic tree, the degree of support for a grouping pattern was evaluated using bootstraps (1000 replicates).

2.2. Chromosomal Locations and Syntenic Analysis

The chromosomal locations of the grapevine SOD genes were verified from the CRIBI database (http://genomes.cribi.unipd.it/grapevine/), and chromosomal images were drawn using MapInspect software [28].

Duplication patterns of the VvSOD genes were assigned based on their locations. Genes located in a 200 kb region of a chromosome or a scaffold were defined as a gene cluster derived from tandem duplications [29]. Segmentally duplicated genes were identified as genes located on duplicated chromosomal blocks and were determined using MCScanX software (http://chibba.pgml.uga.edu/mcscan2/), which detects gene duplication events using an $E$-value threshold of 10⁻⁵ [30].

To detect the selection mode of the segmentally duplicated gene pairs, the ratio of nonsynonymous substitutions to synonymous substitutions was evaluated. Firstly, the CDSs (nucleotide coding sequences) of the VvSOD genes were aligned using Muscle in MEGA6.06. Then, nonsynonymous substitutions ($\text{Ka}$), synonymous substitutions ($\text{Ks}$), and the ratio between them ($\text{Ka}/\text{Ks}$) were calculated using MEGA 6.06.

2.3. GO Analysis and Promoter Sequence Analysis

The gene ontology (GO) term IDs of each VvSOD gene member in grapevine were obtained from the CRIBI database (http://genomes.cribi.unipd.it/DATA/V2/annotation/bl2go.annot.txt). The annotations of the GO term IDs were collected from the Gene Ontology Consortium (http://www.geneontology.org). Ten grapevine SOD proteins were assessed based on their molecular functions, biological processes, and cellular localizations.

The 2000 bp upstream sequences of the coding region of each VvSOD gene were downloaded from the CRIBI database. Then, the cis-elements in the promoters of each VvSOD gene were predicted using the PlantCARE server (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [31].

2.4. Prediction of Potential Regulatory Interactions between TFs and VvSODs

The TFs which possess overrepresented targets in the VvSOD gene sets were detected from the Plant Transcriptional Regulatory Map (PlantRegMap) with a regulation prediction tool (http://plantregmap.cbi.pku.edu.cn/regulation_prediction.php) [32]. 2000 bp upstream sequences of the coding region of each VvSOD gene were used as the input.

2.5. Microarray Data Analysis

Microarray gene expression profiles covering most of the grapevine organs at different stages were downloaded from Gene Expression Omnibus (GEO, available online: https://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE36128 [33]. The collected plant organs were as follows: bud, inflorescence, tendril, leaf, stem, root, developing berry, withering berry, seed, rachis, anther, carpel, petal, pollen, and seedling.

Expression analyses of SOD genes and TF genes in response to abiotic stresses were based on microarray data (including accession number GSE31594 and GSE31675) downloaded from the NCBI GEO datasets. Four abiotic stresses, including cold, salt, heat, and polyethylene glycol (PEG), were analyzed. The fold changes compared with the corresponding control in each experiment were used to generate heatmaps using the R package (https://www.r-project.org/) pheatmap.

2.6. Coexpression Network and Interaction of VvSODs and TF Genes

For the four stress treatments, the Pearson correlation coefficient (PCC) value (Table S4) was calculated between each pair of VvSODs and TF genes using gene expression values from microarray data in SPSS v 20.0 (IBM Corp., Armonk, NY, US) [28]. Coexpressed gene pairs were filtered with a PCC cutoff of 0.95. The coexpression network map was generated by Cytoscape 3.3 (http://www.cytoscape.org).

2.7. Plant Material and Experimental Treatments

For the cold, heat, and salt treatments, six-week-old “Pinot Noir” grapevine (V. vinifera, the sequenced genotype PN40024) plantlets were used in this study. The plantlets were maintained at a photoperiod of 16 h light/8 h dark and a temperature of 23°C. Low temperature and heat stress treatments were applied by placing the plants at 4°C (drug storage box, HYC-360, Haier) and 42°C (intelligent artificial climate box, RXZ-380C, Jiangnan Instrument Factory, Ningbo) with a photoperiod of 16 h light/8 h dark, respectively [28]. Leaf samples were harvested from both treatments at 0, 6, 12, and 24 h after initiation of the treatments. Plants subjected to 200 mmol/L NaCl were sampled at 0, 6, 12, and 24 h after treatment. Three independent biological replications were performed for each treatment. Tissue samples were immediately frozen in liquid nitrogen and stored at −70°C until analysis.

For drought stress treatment, in vitro-grown grapevine plants (PN40024) were maintained on half-strength Murashige and Skoog medium containing 0.3 mg/L indole 3-butyric acid (IBA) and placed at 25°C in a culture room under a photoperiod of 16/8 h and a light intensity of 100 μmol·m^-2·s^-1. Five-week-old tissue-cultured grapevine plants were transplanted into pots and acclimatized in the growth chamber with 16 h light at 24°C/8 h dark at 22°C and 70%–80% relative humidity [34]. Plants grown in pots were fully watered for the first four weeks, after which water was withheld to impose the drought treatment. The grapevine leaf samples were collected at 0 (control), 2, 4, and 8 days after the treatment, respectively.

2.8. Total RNA Extraction and cDNA Synthesis

Total RNA was extracted from the leaves using the Plant Total RNA Isolation Kit Plus. The total RNA was eluted in 30 μL of RNase-free water and stored at −70°C. The purity and density of the extracted RNA was assessed on a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, USA) using RNase-free water as a blank.

Single-stranded cDNA was synthetized using a Takara PrimeScript RT reagent kit with a gDNA eraser (Takara, Dalian, China) and oligo dT primers as described by the manufacturer’s protocol. Approximately 2 μg of total RNA in a single 20 μL reaction was converted to single-stranded cDNA using standard thermal cycling conditions.

2.9. qRT-PCR Analysis of VvSOD Genes

qRT-PCR analyses were conducted on a Roche LightCycler 480 PCR system using Lightcycler 480 SYBR Green I Master. The oligonucleotide primers (Table S3) for amplifying specific VvSODs were designed using the PrimerQuest Tool (http://sg.idtdna.com/Primerquest/Home/Index). The 10 μL reaction volume contained 5 μL SYBR Green I Master, 1 μL of diluted cDNA, 1 μL of each primer (10 μM), and the addition of ddH₂O to bring the total volume to 10 μL. The thermocycle protocol was as follows: denaturation at 95°C for 10 min, followed by 40 cycles of 95°C for 10 s, 58°C for 20 s, and 72°C for 10 s. At the end, a melting curve was generated from 65 to 97°C. Three replicates were used for each sample. Expression levels were calculated using the 2^-ΔΔT method and the housekeeping gene (actin-101-like, VIT_012s0178g00200) was used as a reference gene for normalizing the expression of the VvSOD genes [35]. Fold differences were visualized by normalizing all of the data based on setting the expression level at 0 h as a value of 1 (values above 1 and below 1 were considered as up- and downregulated, respectively) [36].

3. Results

3.1. Genome-Wide Identification of the SOD Gene Family in Grapevine

Based on a BLAST search using the known A. thaliana SOD protein sequences and the identification by the Pfam database, a total of 10 SOD genes were identified in the grapevine genome. Based on the domain analysis, the 10 grapevine SOD genes were classified into 3 groups (Cu/ZnSODs, FeSODs, and MnSODs), including 6 VvSOD proteins with a copper-zinc domain (Pfam: 00080) and 4 SOD proteins with an iron/manganese superoxide dismutase alpha-hairpin domain (Pfam: 00081) or an iron/manganese superoxide dismutase, C-terminal domain (Pfam: 02777). Based on the domain analysis and chromosome location, the 10 grapevine genes were named as follows: VvCSD1, VvCSD2, VvCSD3, VvCSD4, VvCSD5, VvCSD6, VvMSD1, VvMSD2, VvFSD1, and VvFSD2 (Table 1).

Table 1

Chromosome location, subcellular prediction, and physicochemical properties of grape SOD proteins.

Gene name	Gene ID	ORF length	Chromosomal localization	Molecular weight	Theoretical pI	Aliphatic index	Grand average of hydropathicity (GRAVY)	Subcellular prediction
VvCSD1	VIT_202s0025g04830.1	969	chr2(4357308 4364558)	33844.25	5.66	88.66	0.007	Cytoplasmic
VvCSD2	VIT_206s0061g00750.1	666	chr6(18276905 18281206)	22537.5	5.87	92.22	0.067	Cytoplasmic
VvCSD3	VIT_208s0007g07280.1	471	chr8(20888649 20892112)	15786.74	7.19	83.78	-0.164	Cytoplasmic
VvCSD4	VIT_214s0030g00830.1	381	chr14(5076284 5079112)	12926.4	5.33	81.19	-0.202	Cytoplasmic
VvCSD5	VIT_214s0030g00950.1	351	chr14(5248285 5250067)	12149.69	4.93	84.83	0.157	Cytoplasmic
VvCSD6	VIT_214s0036g01320.1	240	chr14(12383065 12384913)	8594.7	5.81	83.67	-0.394	Cytoplasmic
VvMSD1	VIT_206s0004g07950.1	687	chr6(8686646 8687855)	25615.29	7.14	94.56	-0.358	Mitochondrial
VvMSD2	VIT_213s0067g02990.1	852	chr13(1602822 1606910)	31456.18	9.29	86.08	-0.364	Mitochondrial
VvFSD1	VIT_210s0042g00100.1	990	chr10(12756538 12763591)	37779.23	8.62	77.11	-0.522	Chloroplast
VvFSD2	VIT_216s0013g00260.1	681	chr16(5200339 5206132)	25268.79	8.48	83.76	-0.375	Chloroplast

The physicochemical analysis indicated that the length of the amino acid sequences, MWs, pIs, aliphatic indexes, and GRAVY values varied among these VvSOD proteins. Considerable variations in the number of amino acids among these VvSOD proteins were observed, ranging from 79 to 329 aa. The predicted molecular weight of the 10 VvSOD proteins varied from 8594.7 to 37779.23 Da. The results revealed that four Fe-MnSODs were basic and five Cu/ZnSODs were acidic, except for the slightly basic VvCSD3. The GRAVY numeric values of the three VvSOD proteins (VvCSD1, VvCSD2, and VvCSD5) were positive, indicating that they are hydrophobic proteins and the rest of the VvSODs are hydrophilic proteins (Table 1).

The ProtComp9.0 program was used to predict the subcellular localizations of the VvSOD proteins. Among them, six Cu/ZnSODs were cytoplasmic, two MnSODs were in the mitochondria, and two FeSODs were located on the chloroplast (Table 1). These results were in accordance with a previous study [8, 37].

The MEME server was used for conserved motif analysis, and a total of eight conserved motifs were identified (Figure 1). Among them, motif 3 was observed in all six Cu/ZnSODs, while motifs 5, 6, and 7 were found in four Fe-MnSODs (Figure 2). The exon numbers varied among the 10 VvSODs, ranging from 4 to 9 (Figure 2). It is worth noting that the tandemly duplicated gene pair (VvCSD4 and VvCSD5) and segmentally duplicated gene pair (VvMSD1 and VvMSD2) both exhibited the same exon numbers.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

3.2. Phylogenetic Analysis of SOD Genes

Twenty-eight SOD protein sequences from grapevine and A. thaliana (Table S1) were used to construct an unrooted phylogenetic tree (Figure 3). The SODs were divided into three major groups, namely, Cu/ZnSODs, MnSODs, and FeSODs. The FeSODs and MnSODs were clustered within a large clade with a high bootstrap value (99%), implying that the FeSODs and MnSODs originated from a common ancestor. Additionally, species-specific duplications occurred between VvCSD4 and VvCSD5 and VvMSD1 and VvMSD2.

[figure omitted; refer to PDF]

3.3. Chromosomal Locations and Syntenic Analysis

The 10 VvSOD genes were located on seven chromosomes, including VvCSD1 on chromosome 2; VvCSD2 and VvMSD1 on chromosome 6; VvCSD4, VvCSD5, and VvCSD6 on chromosome 14; and VvCSD3, VvFSD1, VvMSD2, and VvFSD2 on chromosomes 8, 10, 13, and 16, respectively (Figure 4). One segmental duplication event was detected between gene VvMSD1 and VvMSD2 on chromosomes 6 and 13. Moreover, VvCSD4 and VvCSD5 were identified as a pair of tandemly duplicated genes. The $\text{Ka}/\text{Ks}$ ratios were calculated to better understand the duplication and functional divergence of the duplicated VvSOD genes during their evolutionary course. As shown in Table S2, the $\text{Ka}/\text{Ks}$ ratio of the segmentally duplicated gene pair (VvMSD1 and VvMSD2) was 0.29, indicating that they underwent purifying selection. The $\text{Ka}/\text{Ks}$ ratio of the tandemly duplicated gene pair (VvCSD4 and VvCSD5) was 1.29, which indicated that they were driven by positive selection during their evolutionary process.

[figure omitted; refer to PDF]

3.4. GO Analysis of the VvSOD Proteins

Biological processes, molecular functions, and cellular components of genes are characteristics of genes or gene products that elucidate the diverse molecular functions of proteins [38]. GO analysis was performed to further characterize the predicted functions of the 10 VvSOD proteins (Figure 5). The “molecular function” data showed that all VvSOD genes were involved in the “superoxide dismutase activity” (GO:0004784). According to the “biological process” results, all of the VvSOD genes participated in the “oxidation-reduction process” (GO:0055114). Moreover, the grapevine SOD genes may be involved in the biological processes responding to biotic stimulus and abiotic stimulus, such as strong light, ozone, salt, and UV stresses. In addition, they may also participate in some biological metabolic processes, such as superoxide metabolic processes, tRNA metabolic processes, and rRNA processing.

[figure omitted; refer to PDF]

3.5. Analysis of cis-Elements in Putative VvSOD Gene Promoters

To further determine the regulatory roles of VvSODs under various stresses, the cis-elements in the VvSOD gene promoter sequences were researched. The cis-elements were divided into four major subgroups: stress-responsive, hormone-responsive, light-responsive, and MYB binding site (Figure 6). Among 10 VvSODs, nine of them possessed several MYB binding sites. Furthermore, we found that abundant light-responsive cis-elements existed in the VvSODs, particularly in the VvMSD2 gene with a number of 19.

[figure omitted; refer to PDF]

3.6. Prediction of Potential Regulatory Interactions between TFs and VvSODs

To investigate the potential regulatory interactions between TFs and VvSODs, regulation prediction was done. A total of 21 TFs which possess overrepresented targets in the input gene set under a cutoff $p$ value ≤ 0.05 were found. Among them, 10, 2, 9, 10, 9, 5, 1, 13, 3, and 8 TFs targeted VvCSD1, VvCSD2, VvCSD3, VvCSD4, VvCSD5, VvCSD6, VvFSD1, VvFSD2, VvMSD1, and VvMSD2, respectively (Table S6).

3.7. Expression Analysis of the VvSOD Genes Based on Microarray Data

To gain more insight into the role that VvSODs play in plant growth and development, we analyzed the expression patterns of the VvSOD genes in different grapevine tissues and organs based on the microarray data from 54 grapevine samples [33]. Interestingly, VvMSD2 displayed consistently high expression in all 54 tissues, whereas VvMSD1 exhibited lower expression levels compared with the others (Figure 7). VvCSD4 and VvCSD5 exhibited similar expression patterns in all the tested tissues. This similarity was also evident in gene VvFSD1 and VvFSD2. Additionally, two genes (VvMSD1 and VvMSD2) demonstrated distinct tissue-specific expression patterns.

[figure omitted; refer to PDF]

Furthermore, four Affymetrix microarray datasets were specifically analyzed for determining the expression profiles of VvSODs and TF genes in response to abiotic stresses, which included PEG, salt, heat, and cold. Six out of 10 probes corresponding to VvSOD genes were found (VvCSD1, VvCSD2, VvCSD3, VvCSD4, VvMSD2, and VvFSD1). Eight out of 21 probes corresponding to TF genes were found (VIT_18s0122g00380, VIT_17s0000g01260, VIT_08s0007g06270, VIT_04s0008g01470, VIT_05s0020g03880, VIT_19s0014g01680, VIT_15s0048g02870, and VIT_10s0003g00140). In the cold, salt, and PEG stress, the expression levels of VvCSD1 and VvCSD2 were higher than those of the other four genes (Figure 8(a)). VvFSD1 was expressed at an extremely low level in these three treatments. However, VvCSD3 and VvFSD1 were highly induced under the heat stress treatment. In each abiotic stress, there existed a single TF gene which responded strongly. However, the other was expressed at the similar level (Figure 8(b)). For example, listed genes, VIT_10s0003g00140 in cold stress, VIT_15s0048g02870 in salt and PEG stress, and VIT_08s0007g06270 in heat stress, were all altered strongly by 8 times compared with the control. VIT_10s0003g00140, VIT_15s0048g02870 and VIT_08s0007g06270 belong to ERF, HD-ZIP, and SBP TF family, respectively (Table S6).

[figures omitted; refer to PDF]

3.8. Expression Analysis of VvSOD Genes Based on qRT-PCR

Previous studies have demonstrated the important roles that SOD genes play in the plant response to abiotic stresses. The analysis of microarray data also indicated that the VvSOD genes were induced by different stresses. To gain more insight into the relative expressional diversity of eight members of the VvSOD gene family to cold, heat, salt, and drought, we detected the expression levels using qRT-PCR (Figure 9). Our results demonstrated that all eight VvSOD genes were induced noticeably by the four stresses, exhibiting different expression levels. All eight members were upregulated in the cold treatment, except for VvCSD3 and VvFSD2 (Figure 9). Among the six upregulated genes, VvCSD1, VvCSD2, and VvFSD1 were altered strongly by 3–7.7 times compared with the control. Upon heat stress, expression of VvCSD1 increased by 12-fold, 11-fold, and 9-fold at 6 h, 12 h, and 24 h, respectively. VvCSD2, VvCSD3, VvFSD2, and VvMSD2 were also upregulated at 6 h with a lower expression level compared with that in VvCSD1. However, the other three genes (VvCSD4, VvCSD6, and VvFSD1) were downregulated during the treatment periods (Figure 9). In the salt treatment, the eight genes showed the lowest expression levels compared with the other three stresses. The expression levels of most of the genes exhibited slight variations, except VvCSD3 and VvCSD6, which were upregulated more than 2-fold. Upon drought stress, the expression level of VvCSD1, VvCSD2, VvCSD4, and VvCSD6 peaked at 12 h, but six genes (VvCSD1, VvCSD2, VvCSD3, VvFSD1, VvFSD2, and VvMSD2) were downregulated at 24 h. Only two genes (VvCSD3 and VvMSD2) were downregulated at 6 h and maintained similar levels at the following time points.

[figures omitted; refer to PDF]

3.9. Coexpression Network of the VvSOD Gene Family and TF Genes

To determine the potential coexpression relationships between the VvSOD genes and TF genes, we conducted one coexpression network based on the microarray data under four stress treatments ($p$ value < 0.05, and $\text{PCC}<-0.95$ or >0.95). It showed that VvSOD genes and TF genes corresponded both positively and negatively (Figure 10). A total of 11 gene pairs under three treatments (heat, salt, and drought) were correlated with each other, including five VvSODs (VvCSD2, VvCSD3, VvCSD4, VvMSD2, and VvFSD1) and eight TF genes (VIT_05s0020g03880, VIT_19s0014g01680, VIT_10s0003g00140, VIT_17s0000g01260, VIT_08s0007g06270, VIT_04s0008g01470, VIT_15s0048g02870, VIT_18s0122g00380). Under heat treatment, three gene pairs showed a negative correlation (Figure 10(a)). Only one gene pair was positively correlated with one another under salt stress (VvFSD1 and VIT_08s0007g06270) (Figure 10(b)). Interestingly, correlated gene pair (VvCSD3 and VIT_17s0000g01260, noted by yellow node) was also predicted to have regulation interaction (Table S6). Additionally, the result inferred that there was a negative correlation between VIT_04s0008g01470 with two VvSOD genes (VvCSD3 and VvMSD2) under PEG treatment (Figure 10(c)).

[figures omitted; refer to PDF]

4. Discussion

With the development of second-generation sequencing technologies and the elucidation of the multiple functions of SOD genes, research into the identification of SOD genes has progressed rapidly. The SOD gene family has been identified in many species, such as Arabidopsis [20], longan [17], rice [39], poplar [9], banana [22], sorghum [40], tomato [8], cotton [23, 41], and cucumber [37].

In the present study, we identified 10 SOD genes in the grapevine genome that could be classified into three SOD types (six Cu/ZnSODs, two FeSODs, and two MnSODs) based on the domain and phylogenetic tree analyses. The 10 SOD genes possessed different intron numbers. A previous study reported that Cu/Zn-SODs containing different intron numbers and positions showed no exon-intron structural similarities in related species [42]. It is worth noting that one tandem pair of VvSOD genes (VvMSD1 and VvMSD2) exhibited similar intron/exon organization patterns.

The promoters of VvSODs contain large amounts of light-responsive cis-elements that may contribute to light stress resistance. Gupta et al. detected increased resistance to oxidative stress in transgenic plants with an overexpressed Cu/Zn-SOD gene [43]. Transgenic A. thaliana plants overexpressing a miR398-resistant form of CSD2 accumulated more CSD2 mRNA than plants overexpressing a regular CSD2 and as a result were much more tolerant of intensive light [44]. Similar observations have also been reported in Nicotiana tabacum and N. plumbaginifolia [45, 46]. All VvSOD genes except for VvFSD2 possessed the MYB binding site. It had been reported that MYB transcription factors played regulatory roles in developmental processes and defense responses in plants [47]. The existence of abundant MYB-related cis-elements indicated that it may involve in the regulation of the expression of VvSOD genes.

By predicting TFs which may regulate VvSODs, several TFs were obtained (Table S6). These TFs belonged to a different TF protein family, including SBP, MYB, WRKY, ERF, CPP, LBD, TCP, HD-ZIP, AP2, and Dof family. Many of them had been found involve in a variety of biological functions like flower development, biotic and abiotic stresses, cell proliferation, hormonal response, and carbohydrate metabolism [48–54]. It has been reported that TFs regulate the target genes by both activate way and repress way [55, 56]. The coexpression results of VvSODs and TF genes also showed evidence for this conclusion. TF regulation prediction showed that VIT_17s0000g01260 may interact with VvCSD3. Coexpression result identified this prediction with a negative relationship. However, the roles of TFs in VvSOD regulatory networks, as well as their interaction mechanism, remain to be fully elucidated through experimental verification and bioinformatic analysis.

Gene duplication provides new genetic material and novel genes and occurs via tandem, segmental, or genome-wide duplication [57–60]. In our study, only one segmental duplication event (VvMSD1 and VvMSD2) and one tandem duplication event (VvCSD4 and VvCSD5) were identified.

The expression profiles of the 10 VvSOD genes in the different tissues indicated that most of the VvSOD genes were expressed, and a few genes displayed tissue-specific expression patterns. VvCSD4 and VvCSD5, identified as tandemly duplicated gene pairs, exhibited similar gene expression patterns. Moreover, VvMSD1 and VvMSD2, the segmentally duplicated gene pair, showed different expression patterns in all the tested tissues. Li et al. reported that the attribution of duplication mode to the expression divergence implies a different evolutionary course of duplicated genes [61]. Adams et al. found that organ-specific and biased expression or silencing of duplicated gene pairs existed in Gossypium hirsutum [62]. A similar phenomenon was also reported in wheat [63].

The qRT-PCR results revealed that every VvSOD gene responded to at least one abiotic stress treatment performed in our study (Figure 9). Notably, only VvCSD1 was activated strongly under heat stress. Moreover, VvCSD1 also expressed highly under cold stress. This suggest that VvCSD1 may play a vital role in the low/high temperature defense. VvCSD6 showed the highest expression level compared with other genes under drought stress. Promoter analysis revealed that VvCSD6 harbored 3 MYB binding sites which involved in drought inducibility, which could explain the significant expression of VvCSD6 under drought stress. Under salt stress, only the genes VvCSD3 and VvCSD6 altered obviously. Taken together, the expression pattern under various stress conditions suggested different VvSOD genes could perform its own function under different stresses. This research constitutes the first comprehensive study of the SOD gene family in grapevine and may provide valuable information for elucidating the classification, evolution, and putative functions of the grapevine SOD family on the whole-genome scale.

5. Conclusion

In the present study, we identified 10 SOD gene grapevine, including three types (6 Cu/Zn SODs, 2 FeSODs, and 2 MnSODs), which were unevenly distributed on 12 chromosomes. Syntenic analysis revealed that VvMSD1 and VvMSD2 were derived from segmental duplication, and VvCSD4 and VvCSD5 belong to a pair of tandemly duplicated genes. Promoter analysis indicated that VvSODs contain large amounts of light-responsive cis-elements which may contribute to light stress resistance. Furthermore, the expression profiles of the 10 VvSOD genes in the different tissues and abiotic stresses indicate that a few genes displayed tissue-specific expression patterns and VvSODs play roles in different aspects of abiotic stress. Finally, coexpression networks of VvSODs and TF genes under the four abiotic stresses were generated based on the microarray data. VIT_17s0000g01260 predicted to regulate VvCSD3, coexpression analysis verified the corelationship of these two genes. This suggested that VIT_17s0000g01260 may participate in the regulation of VvCSD3. However, the regulation mechanism needs further investigation.

This was the first comprehensive study of the SOD gene family in grapevine and may provide valuable information for understanding the classification, evolution, and putative functions of the grapevine SOD family on the genome-wide scale.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

References

[1] Y. Ma, J. Wang, Y. Zhong, F. Geng, G. R. Cramer, Z. M. Cheng, "Subfunctionalization of cation/proton antiporter 1 genes in grapevine in response to salt stress in different organs," Horticulture Research, vol. 2 no. 1, article 15031,DOI: 10.1038/hortres.2015.31, 2015.

[2] M. Tu, X. Wang, Y. Zhu, D. Wang, X. Zhang, Y. Cui, Y. Li, M. Gao, Z. Li, Y. Wang, X. Wang, "VlbZIP30 of grapevine functions in dehydration tolerance via the abscisic acid core signaling pathway," Horticulture Research, vol. 5 no. 1,DOI: 10.1038/s41438-018-0054-x, 2018.

[3] J. P. Londo, A. P. Kovaleski, J. A. Lillis, "Divergence in the transcriptional landscape between low temperature and freeze shock in cultivated grapevine (Vitis vinifera)," Horticulture Research, vol. 5 no. 1,DOI: 10.1038/s41438-018-0020-7, 2018.

[4] T. Karuppanapandian, H. W. Wang, N. Prabakaran, K. Jeyalakshmi, M. Kwon, K. Manoharan, W. Kim, "2,4-dichlorophenoxyacetic acid-induced leaf senescence in mung bean (Vigna radiata L. Wilczek) and senescence inhibition by co-treatment with silver nanoparticles," Plant Physiology and Biochemistry, vol. 49 no. 2, pp. 168-177, DOI: 10.1016/j.plaphy.2010.11.007, 2011.

[5] S. H. Lee, N. Ahsan, K. W. Lee, D. H. Kim, D. G. Lee, S. S. Kwak, S. Y. Kwon, T. H. Kim, B. H. Lee, "Simultaneous overexpression of both CuZn superoxide dismutase and ascorbate peroxidase in transgenic tall fescue plants confers increased tolerance to a wide range of abiotic stresses," Journal of Plant Physiology, vol. 164 no. 12, pp. 1626-1638, DOI: 10.1016/j.jplph.2007.01.003, 2007.

[6] J. M. Tepperman, P. Dunsmuir, "Transformed plants with elevated levels of chloroplastic SOD are not more resistant to superoxide toxicity," Plant Molecular Biology, vol. 14 no. 4, pp. 501-511, DOI: 10.1007/BF00027496, 1990.

[7] R. G. Alscher, N. Erturk, L. S. Heath, "Role of superoxide dismutases (SODs) in controlling oxidative stress in plants," Journal of Experimental Botany, vol. 53 no. 372, pp. 1331-1341, DOI: 10.1093/jexbot/53.372.1331, 2002.

[8] K. Feng, J. Yu, Y. Cheng, M. Ruan, R. Wang, Q. Ye, G. Zhou, Z. Li, Z. Yao, Y. Yang, Q. Zheng, H. Wan, "The SOD gene family in tomato: identification, phylogenetic relationships, and expression patterns," Frontiers in Plant Science, vol. 7,DOI: 10.3389/fpls.2016.01279, 2016.

[9] J. J. Molina-Rueda, C. J. Tsai, E. G. Kirby, "The Populus superoxide dismutase gene family and its responses to drought stress in transgenic poplar overexpressing a pine cytosolic glutamine synthetase (GS1a)," PLoS One, vol. 8 no. 2, article e56421,DOI: 10.1371/journal.pone.0056421, 2013.

[10] J. J. P. Perry, D. S. Shin, E. D. Getzoff, J. A. Tainer, "The structural biochemistry of the superoxide dismutases," Biochimica et Biophysica Acta (BBA) - Proteins and Proteomics, vol. 1804 no. 2, pp. 245-262, DOI: 10.1016/j.bbapap.2009.11.004, 2010.

[11] W. Wang, M. X. Xia, J. Chen, R. Yuan, F. N. Deng, F. F. Shen, "Gene expression characteristics and regulation mechanisms of superoxide dismutase and its physiological roles in plants under stress," Biochemistry, vol. 81 no. 5, pp. 465-480, DOI: 10.1134/S0006297916050047, 2016.

[12] H.-D. Youn, E.-J. Kim, J.-H. Roe, Y. C. Hah, S.-O. Kang, "A novel nickel-containing superoxide dismutase from Streptomyces spp," Biochemical Journal, vol. 318 no. 3, pp. 889-896, DOI: 10.1042/bj3180889, 1996.

[13] F. Wolfe-Simon, D. Grzebyk, O. Schofield, P. G. Falkowski, "The role and evolution of superoxide dismutases in algae1," Journal of Phycology, vol. 41 no. 3, pp. 453-465, DOI: 10.1111/j.1529-8817.2005.00086.x, 2005.

[14] B. Wang, U. Lüttge, R. Ratajczak, "Specific regulation of SOD isoforms by NaCl and osmotic stress in leaves of the C3 halophyte Suaeda salsa L," Journal of Plant Physiology, vol. 161 no. 3, pp. 285-293, DOI: 10.1078/0176-1617-01123, 2004.

[15] M. Pilon, K. Ravet, W. Tapken, "The biogenesis and physiological function of chloroplast superoxide dismutases," Biochimica et Biophysica Acta-Bioenergetics, vol. 1807 no. 8, pp. 989-998, DOI: 10.1016/j.bbabio.2010.11.002, 2011.

[16] A. C. Asensio, M. Gil-Monreal, L. Pires, Y. Gogorcena, P. M. Aparicio-Tejo, J. F. Moran, "Two Fe-superoxide dismutase families respond differently to stress and senescence in legumes," Journal of Plant Physiology, vol. 169 no. 13, pp. 1253-1260, DOI: 10.1016/j.jplph.2012.04.019, 2012.

[17] Y. L. Lin, Z. X. Lai, "Superoxide dismutase multigene family in longan somatic embryos: a comparison of CuZn-SOD, Fe-SOD, and Mn-SOD gene structure, splicing, phylogeny, and expression," Molecular Breeding, vol. 32 no. 3, pp. 595-615, DOI: 10.1007/s11032-013-9892-2, 2013.

[18] D. Zhu, J. G. Scandalios, "Maize mitochondrial manganese superoxide dismutases are encoded by a differentially expressed multigene family," Proceedings of the National Academy of Sciences, vol. 90 no. 20, pp. 9310-9314, DOI: 10.1073/pnas.90.20.9310, 1993.

[19] K. H. Baek, D. Z. Skinner, P. Ling, X. Chen, "Molecular structure and organization of the wheat genomic manganese superoxide dismutase gene," Genome, vol. 49 no. 3, pp. 209-218, DOI: 10.1139/g05-102, 2006.

[20] D. J. Kliebenstein, R. A. Monde, R. L. Last, "Superoxide dismutase in Arabidopsis: an eclectic enzyme family with disparate regulation and protein localization," Plant Physiology, vol. 118 no. 2, pp. 637-650, DOI: 10.1104/pp.118.2.637, 1998.

[21] E. Filiz, H. Tombuloğlu, "Genome-wide distribution of superoxide dismutase (SOD) gene families in Sorghum bicolor," Turkish Journal of Biology, vol. 39 no. 1, pp. 49-59, DOI: 10.3906/biy-1403-9, 2015.

[22] X. Feng, Z. Lai, Y. Lin, G. Lai, C. Lian, "Genome-wide identification and characterization of the superoxide dismutase gene family in Musa acuminata cv. Tianbaojiao (AAA group)," BMC Genomics, vol. 16 no. 1,DOI: 10.1186/s12864-015-2046-7, 2015.

[23] W. Wang, M. Xia, J. Chen, F. Deng, R. Yuan, X. Zhang, F. Shen, "Genome-wide analysis of superoxide dismutase gene family in Gossypium raimondii and G. arboreum," Plant Gene, vol. 6 no. C, pp. 18-29, DOI: 10.1016/j.plgene.2016.02.002, 2016.

[24] W. Wang, X. Zhang, F. Deng, R. Yuan, F. Shen, "Genome-wide characterization and expression analyses of superoxide dismutase (SOD) genes in Gossypium hirsutum," BMC Genomics, vol. 18 no. 1,DOI: 10.1186/s12864-017-3768-5, 2017.

[25] E. Gasteiger, A. Gattiker, C. Hoogland, I. Ivanyi, R. D. Appel, A. Bairoch, "ExPASy: the proteomics server for in-depth protein knowledge and analysis," Nucleic Acids Research, vol. 31 no. 13, pp. 3784-3788, DOI: 10.1093/nar/gkg563, 2003.

[26] T. L. Bailey, J. Johnson, C. E. Grant, W. S. Noble, Nucleic Acids Research, vol. 43 no. W1, pp. W39-W49, DOI: 10.1093/nar/gkv416, 2015.

[27] K. Tamura, G. Stecher, D. Peterson, A. Filipski, S. Kumar, "MEGA6: Molecular Evolutionary Genetics Analysis version 6.0," Molecular Biology and Evolution, vol. 30 no. 12, pp. 2725-2729, DOI: 10.1093/molbev/mst197, 2013.

[28] Y. Xi, J. Liu, C. Dong, Z. M. (. M.). Cheng, "The CBL and CIPK gene family in grapevine (Vitis vinifera): genome-wide analysis and expression profiles in response to various abiotic stresses," Frontiers in Plant Science, vol. 8,DOI: 10.3389/fpls.2017.00978, 2017.

[29] Y. Zhong, Z. M. Cheng, "A unique RPW8-encoding class of genes that originated in early land plants and evolved through domain fission, fusion, and duplication," Scientific Reports, vol. 6 no. 1, article 32923,DOI: 10.1038/srep32923, 2016.

[30] K. Zhu, F. Chen, J. Liu, X. Chen, T. Hewezi, Z. M. Cheng, "Evolution of an intron-poor cluster of the CIPK gene family and expression in response to drought stress in soybean," Scientific Reports, vol. 6 no. 1, article 28225,DOI: 10.1038/srep28225, 2016.

[31] M. Lescot, P. Déhais, G. Thijs, K. Marchal, Y. Moreau, Y. van de Peer, P. Rouzé, S. Rombauts, "PlantCARE, a database of plant cis-acting regulatory elements and a portal to tools for in silico analysis of promoter sequences," Nucleic Acids Research, vol. 30 no. 1, pp. 325-327, DOI: 10.1093/nar/30.1.325, 2002.

[32] J. Jin, F. Tian, D. C. Yang, Y. Q. Meng, L. Kong, J. Luo, G. Gao, "PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants," Nucleic Acids Research, vol. 45 no. D1, pp. D1040-D1045, DOI: 10.1093/nar/gkw982, 2017.

[33] M. Fasoli, S. Dal Santo, S. Zenoni, G. B. Tornielli, L. Farina, A. Zamboni, A. Porceddu, L. Venturini, M. Bicego, V. Murino, A. Ferrarini, M. Delledonne, M. Pezzotti, "The grapevine expression atlas reveals a deep transcriptome shift driving the entire plant into a maturation program," Plant Cell, vol. 24 no. 9, pp. 3489-3505, DOI: 10.1105/tpc.112.100230, 2012.

[34] K. Zhu, X. Wang, J. Liu, J. Tang, Q. Cheng, J. G. Chen, Z. M. Cheng, "The grapevine kinome: annotation, classification and expression patterns in developmental processes and stress responses," Horticulture Research, vol. 5 no. 1,DOI: 10.1038/s41438-018-0027-0, 2018.

[35] M. Wang, A. Vannozzi, G. Wang, Y. H. Liang, G. B. Tornielli, S. Zenoni, E. Cavallini, M. Pezzotti, Z. M. Cheng, "Genome and transcriptome analysis of the grapevine (Vitis vinifera L.) WRKY gene family," Horticulture Research, vol. 1 no. 1, article 14016,DOI: 10.1038/hortres.2014.16, 2014.

[36] J. Liu, N. Chen, F. Chen, B. Cai, S. Dal Santo, G. Tornielli, M. Pezzotti, Z. M. Cheng, "Genome-wide analysis and expression profile of the bZIP transcription factor gene family in grapevine (Vitis vinifera)," BMC Genomics, vol. 15 no. 1,DOI: 10.1186/1471-2164-15-281, 2014.

[37] Y. Zhou, L. Hu, H. Wu, L. Jiang, S. Liu, "Genome-wide identification and transcriptional expression analysis of cucumber superoxide dismutase (SOD) family in response to various abiotic stresses," International Journal of Genomics, vol. 2017,DOI: 10.1155/2017/7243973, 2017.

[38] M. Ashburner, C. A. Ball, J. A. Blake, D. Botstein, H. Butler, J. M. Cherry, A. P. Davis, K. Dolinski, S. S. Dwight, J. T. Eppig, M. A. Harris, D. P. Hill, L. Issel-Tarver, A. Kasarskis, S. Lewis, J. C. Matese, J. E. Richardson, M. Ringwald, G. M. Rubin, G. Sherlock, "Gene ontology: tool for the unification of biology. The gene ontology consortium," Nature Genetics, vol. 25 no. 1, pp. 25-29, DOI: 10.1038/75556, 2000.

[39] B. Dehury, K. Sarma, R. Sarmah, J. Sahu, S. Sahoo, M. Sahu, P. Sen, M. K. Modi, G. D. Sharma, M. D. Choudhury, M. Barooah, "In silico analyses of superoxide dismutases (SODs) of rice (Oryza sativa L.)," Journal of Plant Biochemistry and Biotechnology, vol. 22 no. 1, pp. 150-156, DOI: 10.1007/s13562-012-0121-6, 2013.

[40] E. Filiz, H. Tombuloğlu, "Genome-wide distribution of superoxide dismutase (SOD) gene families inSorghum bicolor," Turkish Journal of Biology, vol. 39, pp. 49-59, DOI: 10.3906/biy-1403-9, 2015.

[41] J. Zhang, B. Li, Y. Yang, W. Hu, F. Chen, L. Xie, L. Fan, "Genome-wide characterization and expression profiles of the superoxide dismutase gene family in Gossypium," International Journal of Genomics, vol. 2016,DOI: 10.1155/2016/8740901, 2016.

[42] R. C. Fink, J. G. Scandalios, "Molecular evolution and structure–function relationships of the superoxide dismutase gene families in angiosperms and their relationship to other eukaryotic and prokaryotic superoxide dismutases," Archives of Biochemistry and Biophysics, vol. 399 no. 1, pp. 19-36, DOI: 10.1006/abbi.2001.2739, 2002.

[43] A. S. Gupta, J. L. Heinen, A. S. Holaday, J. J. Burke, R. D. Allen, "Increased resistance to oxidative stress in transgenic plants that overexpress chloroplastic Cu/Zn superoxide dismutase," Proceedings of the National Academy of Sciences, vol. 90 no. 4, pp. 1629-1633, DOI: 10.1073/pnas.90.4.1629, 1993.

[44] R. Sunkar, A. Kapoor, J. K. Zhu, "Posttranscriptional induction of two Cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance," Plant Cell, vol. 18 no. 8, pp. 2051-2065, DOI: 10.1105/tpc.106.041673, 2006.

[45] J. Kurepa, M. V. Montagu, D. Inzé, "Expression of sodCp and sodB genes in Nicotiana tabacum: effects of light and copper excess," Journal of Experimental Botany, vol. 48 no. 12, pp. 2007-2014, DOI: 10.1093/jxb/48.12.2007, 1997.

[46] E. W. Tsang, C. Bowler, D. Hérouart, W. Van Camp, R. Villarroel, C. Genetello, M. Van Montagu, D. Inze, "Differential regulation of superoxide dismutases in plants exposed to environmental stress," Plant Cell, vol. 3 no. 8, pp. 783-792, DOI: 10.1105/tpc.3.8.783, 1991.

[47] C. Yanhui, Y. Xiaoyuan, H. Kun, L. Meihua, L. Jigang, G. Zhaofeng, L. Zhiqiang, Z. Yunfei, W. Xiaoxiao, Q. Xiaoming, S. Yunping, Z. Li, D. Xiaohui, L. Jingchu, D. Xing-Wang, C. Zhangliang, G. Hongya, Q. Li-Jia, "The MYB transcription factor superfamily of Arabidopsis: expression analysis and phylogenetic comparison with the rice MYB family," Plant Molecular Biology, vol. 60 no. 1, pp. 107-124, DOI: 10.1007/s11103-005-2910-y, 2006.

[48] K. Yamasaki, T. Kigawa, M. Inoue, M. Tateno, T. Yamasaki, T. Yabuki, M. Aoki, E. Seki, T. Matsuda, E. Nunokawa, Y. Ishizuka, T. Terada, M. Shirouzu, T. Osanai, A. Tanaka, M. Seki, K. Shinozaki, S. Yokoyama, "A novel zinc-binding motif revealed by solution structures of DNA-binding domains of Arabidopsis SBP-family transcription factors," Journal of Molecular Biology, vol. 337 no. 1, pp. 49-63, DOI: 10.1016/j.jmb.2004.01.015, 2004.

[49] S. P. Pandey, I. E. Somssich, "The role of WRKY transcription factors in plant immunity," Plant Physiology, vol. 150 no. 4, pp. 1648-1655, DOI: 10.1104/pp.109.138990, 2009.

[50] M. Martín-Trillo, P. Cubas, "TCP genes: a family snapshot ten years later," Trends in Plant Science, vol. 15 no. 1, pp. 31-39, DOI: 10.1016/j.tplants.2009.11.003, 2010.

[51] T. Nakano, K. Suzuki, T. Fujimura, H. Shinshi, "Genome-wide analysis of the ERF gene family in Arabidopsis and rice," Plant Physiology, vol. 140 no. 2, pp. 411-432, DOI: 10.1104/pp.105.073783, 2006.

[52] D. Lijavetzky, P. Carbonero, J. Vicente-Carbajosa, "Genome-wide comparative phylogenetic analysis of the rice and Arabidopsis Dof gene families," BMC Evolutionary Biology, vol. 3 no. 1,DOI: 10.1186/1471-2148-3-17, 2003.

[53] Z. Yang, S. Gu, X. Wang, W. Li, Z. Tang, C. Xu, "Molecular evolution of the CPP-like gene family in plants: insights from comparative genomics of Arabidopsis and rice," Journal of Molecular Evolution, vol. 67 no. 3, pp. 266-277, DOI: 10.1007/s00239-008-9143-z, 2008.

[54] R. Stracke, M. Werber, B. Weisshaar, "The R2R3-MYB gene family in Arabidopsis thaliana," Current Opinion in Plant Biology, vol. 4 no. 5, pp. 447-456, DOI: 10.1016/S1369-5266(00)00199-0, 2001.

[55] T. Eulgem, P. J. Rushton, S. Robatzek, I. E. Somssich, "The WRKY superfamily of plant transcription factors," Trends in Plant Science, vol. 5 no. 5, pp. 199-206, DOI: 10.1016/S1360-1385(00)01600-9, 2000.

[56] D. S. Latchman, "Transcription factors: bound to activate or repress," Trends in Biochemical Sciences, vol. 26 no. 4, pp. 211-213, DOI: 10.1016/S0968-0004(01)01812-6, 2001.

[57] A. R. Kersting, E. Bornberg-Bauer, A. D. Moore, S. Grath, "Dynamics and adaptive benefits of protein domain emergence and arrangements during plant genome evolution," Genome Biology and Evolution, vol. 4 no. 3, pp. 316-329, DOI: 10.1093/gbe/evs004, 2012.

[58] M. Buljan, A. Bateman, "The evolution of protein domain families," Biochemical Society Transactions, vol. 37 no. 4, pp. 751-755, DOI: 10.1042/BST0370751, 2009.

[59] X. Yang, G. A. Tuskan, Z.-M. Cheng, "Divergence of the Dof gene families in poplar, Arabidopsis, and rice suggests multiple modes of gene evolution after duplication," Plant Physiology, vol. 142 no. 3, pp. 820-830, DOI: 10.1104/pp.106.083642, 2006.

[60] Y. Zhong, Y. Jia, Y. Gao, D. Tian, S. Yang, X. Zhang, "Functional requirements driving the gene duplication in 12 Drosophila species," BMC Genomics, vol. 14 no. 1,DOI: 10.1186/1471-2164-14-555, 2013.

[61] Z. Li, H. Zhang, S. Ge, X. Gu, G. Gao, J. Luo, "Expression pattern divergence of duplicated genes in rice," BMC Bioinformatics, vol. 10 no. S6,DOI: 10.1186/1471-2105-10-S6-S8, 2009.

[62] K. L. Adams, R. Cronn, R. Percifield, J. F. Wendel, "Genes duplicated by polyploidy show unequal contributions to the transcriptome and organ-specific reciprocal silencing," Proceedings of the National Academy of Sciences, vol. 100 no. 8, pp. 4649-4654, DOI: 10.1073/pnas.0630618100, 2003.

[63] A. Bottley, G. M. Xia, R. M. D. Koebner, "Homoeologous gene silencing in hexaploid wheat," The Plant Journal, vol. 47 no. 6, pp. 897-906, DOI: 10.1111/j.1365-313X.2006.02841.x, 2006.