The Class III Peroxidase (POD) Gene Family in Cassava: Identification, Phylogeny, Duplication, and Expression

The class III peroxidase (POD) enzymes participate in plant development, hormone signaling, and stress responses. However, little is known about the POD family in cassava. Here, we identified 91 cassava POD genes (MePODs) and classified them into six subgroups using phylogenetic analysis. Conserved motif analysis demonstrated that all MePOD proteins have typical peroxidase domains, and gene structure analysis showed that MePOD genes have between one and nine exons. Duplication pattern analysis suggests that tandem duplication has played a role in MePOD gene expansion. Comprehensive transcriptomic analysis revealed that MePOD genes in cassava are involved in the drought response and postharvest physiological deterioration. Several MePODs underwent transcriptional changes after various stresses and related signaling treatments were applied. In sum, we characterized the POD family in cassava and uncovered the transcriptional control of POD genes in response to various stresses and postharvest physiological deterioration conditions. These results can be used to identify potential target genes for improving the stress tolerance of cassava crops.


Introduction
Peroxidases (EC 1.11.1.X) form a large family of enzymes that are widely distributed in living organisms and catalyze the oxidoreduction reaction between hydrogen peroxide (H 2 O 2 ) as an electron acceptor and diverse electron donors, such as auxin, phenolic compounds, or secondary metabolites [1,2]. According to their protein sequences and structure, peroxidases are classified as either non-heme peroxidases or heme peroxidases [3]. The majority of heme peroxidase members are divided into animal and non-animal groups [4]. On the basis of their amino acid sequences and catalytic properties, non-animal heme peroxidases are assigned to one of three large families: class I, II, or III [3,5]. The class

Phylogenetic and Comparative Analyses of PODs in Cassava
The homology and similarity of the POD genes in cassava were determined by performing multiple sequence alignments. Then, the radiation phylogenetic tree of the 91 POD genes was constructed using the neighbor-joining (NJ) method with a bootstrap value of 1000 using MEGA 5.1 (University College Dublin, Dublin, Ireland). The phylogenetic analyses indicated that MePOD genes can be divided into six subgroups on the basis of the observed genetic distance and bootstrap support ( Figure 1). The large subgroups A and D consist of 23 and 24 MePOD members, respectively, whereas the small subgroups C and F contain 9 and 8 MePOD members, respectively. Subgroups B and E are composed of 15 and 12 MePODs members, respectively. These results show that a diversified POD family exists in cassava.

Conserved Motif and Gene Structure Analysis of POD Families in Cassava
The structural features of MePODs were investigated by identifying 10 conserved motifs using the MEME database in accordance with the phylogenetic relationship. Then, the conserved motifs were submitted in their entirety to the InterProScan database for annotation. Eight domains (domains 1, 2, 3, 4, 5, 6, 7, and 9) were noted as POD protein motifs, which are an essential feature of the peroxidase family. On the basis of the motif analyses, 83 MePODs were assigned to one of five subgroups (A-E). Each of these 83 MePODs contains at least nine POD motifs, except for MePOD57 (in subgroup D) and MePOD69 (in subgroup B), which have seven and five motifs, respectively. The presence of these motifs suggests that the identified proteins are characteristic of the POD family ( Figure 2). Subgroup F is distinct from the others: its members contain domains 2, 4, 5, 7, and 8. These results indicate that the proteins assigned to the same subfamilies share similar POD motif characteristics, further supporting their phylogenetic classification as PODs in cassava.

120
Next, the exon-intron structures of cassava POD genes were analyzed. Subgroup F is exon-rich,

121
with five to nine exons, whereas other subfamilies have fewer (between one and four) exons, except 122 for MePOD02 (in subgroup B), which has five exons ( Figure 3). High proportions of POD genes 123 contain four exons in subgroups A, B, C, and D, with four-exon POD genes forming 84%, 67%, 57%, 124 and 76% of the genes in these subgroups, respectively, whereas only 50% of the subgroup E genes 125 have four exons. Generally, POD genes in the same subgroup show similar exon-intron features, Figure 2. The motif analyses of POD family members in cassava according to their evolutionary relationship. The POD motifs were identified by the MEME database. The 10 different colors of the boxes on the right represent diverse conserved motifs, while the gray lines indicate non-conserved sequences.
Next, the exon-intron structures of cassava POD genes were analyzed. Subgroup F is exon-rich, with five to nine exons, whereas other subfamilies have fewer (between one and four) exons, except for MePOD02 (in subgroup B), which has five exons ( Figure 3). High proportions of POD genes contain four exons in subgroups A, B, C, and D, with four-exon POD genes forming 84%, 67%, 57%, and 76% of the genes in these subgroups, respectively, whereas only 50% of the subgroup E genes have four exons. Generally, POD genes in the same subgroup show similar exon-intron features, providing further evidence of their phylogenetic relationship.

Expression Profiles of POD Genes in Different Tissues of Two Cassava Genotypes
The expression levels of MePOD genes in different tissues were investigated by performing RNA-Seq analysis on the storage roots, stems, and leaves of a cultivated variety (Arg7) and wild subspecies (W14). The resulting expression data covered 59 and 56 MePOD genes in the transcriptome dataset of Arg7 and W14, respectively ( Figure 6A; Table S3). Of these genes in Arg7, 15 (25%), 9 (15%), and 8 (14%) MePODs had high transcriptional levels (log2-based > 4) in stems, leaves, and storage roots, respectively. The number of MePODs with high expression (log2-based > 4) in the stems, leaves, and storage roots of W14 was 9 (16%), 7 (13%), and 10 (18%), respectively. Notably, MePOD5 in subgroup E and MePOD89 in subgroup F were strongly expressed (log2-based fold change > 4) in the three diverse tissues of Arg7 and W14. These POD genes may play a molecular role in the development and function of different cassava tissues. Log2-based fold changes were applied to build the heat map using Mev4.9.0 software.

Expression Profiles of Cassava PODs During PPD
Transcriptome analyses were conducted at different postharvest periods of the storage roots of SC124 to examine the possible function of MePOD genes during PPD ( Figure 6C; Table S5)

239
The expansion of a gene family primarily occurs via three kinds of modes: segmental duplication 240 of multiple genes, tandem duplication of individual genes, and whole-genome duplication [36,37].

241
To analyze the duplication modes of the POD genes in cassava, we first identified the chromosomal 242 locations of the MePOD genes. Chromosomal mapping revealed that these genes are widely 243 distributed among 17 chromosomes and one scaffold (cassava has 18 chromosomes in total) ( Figure  Figure 7. Expression profiles of cassava POD genes in the leaves of Arg7 after exposure to MeJA, SA, ABA, H 2 O 2 , salt, osmotic stress (mannitol treatment), cold stress, and Xam. Log2-based qRT-PCR fold changes were used to build the heat map with Mev4.9.0 software. The changes in color represent the relative gene expression level.

Discussion
Given the significant role of PODs in various physiological processes, including responses to biotic and abiotic stresses, it was necessary to scientifically investigate the potential functions of POD genes in cassava, which is an important crop. In this study, we identified 91 PODs in the cassava genome ( Figure 1); thus, cassava has more POD members than Arabidopsis but fewer than rice, Populus trichocarpa, Medicago sativa, maize, and Pyrus bretschneideri [9,[24][25][26][27][28][29][30]. We found that 92% (84/91) of the MePODs have a molecular mass in the range of 30 to 45 kDa, which is in accordance with previous studies [2,7]. Most of the POD genes (89/91) in cassava harbor more than one exon (Figure 3), which is similar to the proportion of single-exon POD genes in Pyrus bretschneideri and Zea mays (PbPRX (90/94) and ZmPRX (89/107), respectively) [29,30]. The similarities in gene structure and motif composition among the members in each MePOD subgroup support the phylogenetic classification presented here.
The expansion of a gene family primarily occurs via three kinds of modes: segmental duplication of multiple genes, tandem duplication of individual genes, and whole-genome duplication [36,37]. To analyze the duplication modes of the POD genes in cassava, we first identified the chromosomal locations of the MePOD genes. Chromosomal mapping revealed that these genes are widely distributed among 17 chromosomes and one scaffold (cassava has 18 chromosomes in total) (Figure 4), which is in accordance with the wide chromosomal distribution of PODs in Arabidopsis, rice, Populus trichocarpa, maize, and Pyrus bretschneideri [9,24,27,29,30]. Secondly, 16 paralogous POD genes were characterized in the cassava genome, indicating that tandem duplication contributed to MePOD expansion. Accumulated evidence has demonstrated that duplication events have been important for gene expansion in the POD family. A total of 37 PRX genes in Populus and 24 POD genes in maize were identified as tandem duplications, further supporting that tandem duplication has been a significant means of POD gene expansion [27,29]. Almost all these paralogous MePODs had low or no expression after drought treatment, but 63% (10 out of 16) from the postharvest transcriptome were expressed (Table S6), of which MePOD-2, -30, -32, -33, -39, and -44 were significantly upregulated and MePOD-34, -56, and -62 were repressed at some time point during the PPD process ( Figure 6C). These results indicate that most of the MePOD genes resulting from tandem duplication-driven expansion are involved in the PPD process of cassava storage roots. Paralogous PRX genes were also found to be involved in other biological processes. In maize, paralogous genes ZmPRX-26, -42, and -75 were induced after NaCl, PEG, SA, or H 2 O 2 treatment [29]. In Chinese pear, the expression of paralogous genes PbPRX-42 and -64 increased during fruit development [30].
The POD family is positively related to the reduced production of hydrogen peroxide and the decreased formation of reactive oxygen species, and the suppression of these species increases plant resistance to stresses [4,5,11,12,19,20]. In this study, the total number of POD genes responding to drought (log2-based fold change > 1) was greater in both the roots and leaves of W14 than that in Arg7 and SC124, suggesting the comprehensive activation of PODs in response to drought in W14 ( Figure 6B). The wild ancestor W14 has been previously confirmed to be more resistant to drought than the two cultivars SC124 and Arg7 [38,39]. Accumulated evidence suggests that the overexpression of POD genes results in increased plant tolerance to drought and osmotic stresses [19,20,40,41]. The activity of POD enzyme was significantly enhanced under drought stress [42]. Consequently, we conclude that the high ratio of MePODs induced by drought in W14 might contribute to its strong drought tolerance.
Previous studies have suggested that ROS production results in the deterioration process in cassava during the postharvest period, and a reduction in ROS accumulation could delay the PPD process [35,43]. The POD family mainly participates in the peroxidative cycle and hydroxylic cycle, resulting in the reduced production of H 2 O 2 and the decreased formation of ROS [4,5,11,12]. Some PRXs have been shown to change in expression during the fruit storage process [44,45]. The activity of POD enzyme significantly increased during cassava PPD process, suggesting their possible role during the postharvest period of cassava [35,44]. In this study, we found that 78% (71 out of 91) of PODs (log2-based fold change > 1) were upregulated in the storage roots of SC124 ( Figure 6C). Interestingly, 13% (12 out of 91) of PODs (log2-based fold change > 1) were induced at all points. Collectively, these results indicate that MePOD genes are involved in the PPD process in cassava storage roots.
Previous research has indicated that PODs can extensively participate in plants' responsse to biotic and abiotic stresses [18][19][20][21][22][23]. Here, we selected nine genes (MePOD-13, -16, -17, -19, -23, -68, -74, -85, and -86) to further examine their expression levels after various treatments (Figure 7; Table S6). These genes are located on different regions of chr7, 13, 3, 17, 8, 15, 10, 9, and 18, respectively (Figure 4). Phylogenetic analysis indicates that MePOD-16, -68, -74, and -85 belong to subgroup A; MePOD-17 belongs to subgroup D; and MePOD-13, -19, -23, and -86 belong to subgroup E (Figure 1). The results show that all nine of the analyzed MePODs were upregulated in response to at least two types of treatments. MePOD17 and MePOD85 (log2-based fold change > 1) were induced by six treatments (MeJA, salt, cold stress, osmotic stress, ABA, and Xam and MeJA, osmotic stress, SA, ABA, H 2 O 2 , and Xam, respectively); MePOD13 was upregulated by five treatments (MeJA, osmotic stress, ABA, SA, H 2 O 2 , and Xam); and MePOD23 and MePOD86 were upregulated by four treatments (MeJA, osmotic stress, SA, and H 2 O 2 and MeJA, SA, H 2 O 2 , and Xam, respectively). Of these, MePOD13 and MePOD23 were induced after H 2 O 2 treatment in cassava leaves (Figure 7) but exhibited the opposite trend of expression during the PPD process in storage roots ( Figure 6C), suggesting their differential roles in diverse tissues. MePOD-13, -19, -23, -68, and -86 (belonging to subgroup E, except for MePOD68) were upregulated by MeJA and SA treatments but downregulated by ABA treatment (Figures 1 and 7). The expression of some PODs has been induced by MeJA and SA treatments in other plant species [2,46]. The opposite direction of expression of these POD genes between MeJA and SA treatments and ABA treatment may be due to the antagonism between MeJA/SA and ABA [47,48]. Whereas ABA plays a prominent role in plants' tolerance to drought stress [38], MeJA-and SA-mediated signaling pathways are also activated under drought stress [49,50]. The induction of these genes by MeJA, SA, and drought suggests their possible involvement in MeJA-and SA-mediated drought responses in cassava. The responses of POD genes to multiple treatments have been observed in other plants. In Arabidopsis, AtPrx33 and AtPrx34 were upregulated after H 2 O 2 and flg22 treatments [51]. In maize, ZmPRX-26, -42, and -71 were induced by H 2 O 2 , salt, and PEG treatments [29]. Phylogenetic analysis of MePODs with AtPrx-33 and -34 and ZmPRX-26, -42, and -71 found that MePOD86 shares a close phylogenetic relationship with ZmPRX71 ( Figure S1), suggesting their functional conservation in multiple treatments. Multiple stresses, such as cold, salt, or PEG, induced the activity of POD enzyme, demonstrating the response of POD genes to environmental stress [52][53][54]. These results suggest that MePODs participate in the response to multiple stresses or related signals and are candidate targets for the genetic improvement of cassava.

Plant Materials and Treatments
Three cassava genotypes, W14, SC124, and Arg7, were planted in the greenhouse of the Chinese Academy of Tropical Agricultural Science (Haikou, China). Their characteristics were described in our previous studies [38,55]. Stem segments containing three nodes were cut from eight-month-old cassava plants and planted in pots, as described in Hu's study [51]. The transcripts of W14 and Arg7 MePOD genes in stems and leaves after being planted for 90 days and middle storage roots after being planted for 150 days were examined by RNA-Seq. After W14, SC124, and Arg7 were cultured for 90 days, they were subjected to drought stress by withholding water for 12 days, after which their leaves and roots were sampled to study the transcriptional responses by RNA-Seq. To examine the expression profiles of MePOD genes after the plants were exposed to stress and related signaling treatments, the 60-day-old Arg7 variety was treated with 100 µM MeJA for 0, 2, 6, 10, and 24 h; 300 mM NaCl for 0 h, 2 h, 6 h, 3 days, and 14 days; a low temperature (4 • C) for 0, 2, 5, 15, and 48 h; 100 µM SA for 0, 2, 6, 10, and 24 h; 200 mM mannitol (to induce osmotic stress) for 0 h, 2 h, 6 h, 3 days, and 14 days; 100 µM ABA for 0, 2, 6, 10, and 24 h; 10% H 2 O 2 for 0, 2, 6, 10, and 24 h; or Xam for 0, 2, 6, 12, and 24 h. Ten-month-old cassava storage roots (CSR) were cut into 5-mm-thick slices and placed into Petri dishes containing wet filter paper for 0, 6, 12, and 48 h to study the expression changes in MePOD genes via RNA-Seq during CSR deterioration. All samples were frozen immediately in liquid nitrogen and stored at −80 • C for RNA-Seq and qRT-PCR.
All the predicted cassava POD genes identified from HMMER and BLAST were confirmed only if they included the POD special domains examined by SMART (available online: http://smart.emblheidelberg.de/) [61]. Multiple sequence alignment of all predicted MePOD protein sequences was performed with Clustal W in BioEdit software [62]. The phylogenetic tree of the full-length MePOD protein sequences was created using MEGA 5.0 (University College Dublin, Dublin, Ireland) with the neighbor-joining (NJ) method, and bootstrap analysis was conducted with 1000 replicates [63].

Protein Properties and Structure Analyses of PODs in Cassava
The ProtParam database (available online: http://web.expasy.org/protparam/) was used to predict the properties, including amino acid numbers, molecular weights (MW), and isoelectric points (pI), of all presumed POD proteins [64]. The motifs were analyzed using the MEME program (available online: http://meme-suite.org/tools/meme), in which the maximum number of motifs was set to 10, the optimum width of motifs was set to 15-50 amino acid residues, and the other settings were kept at default values [65]. Subsequently, these 10 motifs were annotated in InterProScan (available online: http://www.ebi.ac.uk/Tools/pfa/iprscan/) [66]. The gene structures of each MePOD were investigated using GSDS software (available online: http://gsds.cbi.pku.edu.cn/) using each MePOD's genomic DNA sequence and its corresponding CDS sequence, which were retrieved from the cassava genome database [67].

Chromosomal Location and Duplication Pattern Analyses
According to the results of BLASTN in the Phytozome 12.0 cassava database, MePOD genes were mapped to different chromosomes. On the basis of the calculated value of nucleotide sequence similarity and the phylogenetic relationship of cassava POD genes, paralogous genes were identified. The gene duplication pattern of paralogous MePOD genes was determined by the following two criteria: (1) the identity of the aligned region was >90% and (2) the alignment covered >90% of the longer gene. Circos software (Canada's Michael Smith Genome Sciences Center, Vancouver, Canada) was used to draw the duplication events of MePOD genes [68,69]. The values of nonsynonymous substitution (Ka) and synonymous substitution (Ks) were calculated suing DnaSP 5.0 software [70]. Ka/Ks rate > 1 indicates positive evolution, Ka/Ks rate = 1 indicates neutral evolution, and Ka/Ks rate < 1 indicates negative evolution [71].

Transcriptome Analyses of PODs in Cassava
RNA-Seq was used to determine the expression of cassava MePOD genes. Total RNA was isolated from frozen stems, leaves, roots, and storage roots using the plant RNeasy extraction kit (TIANGEN, Beijing, China) and quantified with a NanoDrop 2000c (Thermo Scientific Inc., Waltham, MA, USA). Total RNA (3 µg) of each sample was used to construct the cDNA library according to the Illumina instructions and then sequenced using an Illumina GA II (Illumina Inc., San Diego, USA). The original data processing and analysis methods were described in our previous study [72].

Quantitative Real-Time PCR Analyses
Leaf samples of the Arg7 variety subjected to MeJA, SA, ABA, NaCl, low temperature, mannitol, H 2 O 2 , or Xam treatments were collected to perform qRT-PCR analysis. Total RNA (1 µg) from each sample was used to synthesize the first-strand cDNA using oligo-dT primer by SuperScript reverse transcriptase (Takara, Dalian, China). The cDNA product was diluted to 50 ng·µL −1 , and 1 µL was used for qRT-PCR. The qRT-PCR reaction mixtures (20 µL) contained 0.6 µL of each gene-specific primer (300 nmol·µL −1 ), 10 µL of 2× FastFire qPCR PreMix (Tiangen, Beijing, China), and 7.8 µL of RNase-free water. The qRT-qPCR thermal cycling included cDNA denaturation at 95 • C for 1 min, with 40 cycles of 95 • C for 5 s and 60 • C for 15 s in the Mx3005P Real-Time PCR System (Agilent Inc., Palo Alto, CA, USA) with the SYBR green method. The β-tubulin gene of cassava was chosen as an internal control. All qRT-qPCR experiments were performed in triplicate, and the gene-specific primers used in expression analysis are listed in Table S7. The data obtained from the qRT-qPCR were analyzed with Tukey's post-hoc ANOVA in SPSS 22.0 (SPSS Inc., Chicago, USA) (P < 0.05) after fold treatment with the 2 − Ct method.

Conclusions
In this study, we identified 91 PODs from the cassava genome and studied their basic classification, protein motif, gene structure, chromosomal distribution, and duplication pattern. Comprehensive transcriptional level analyses revealed the involvement of MePODs in biotic and abiotic stress responses, hormone responses, and storage root deterioration. Several MePOD genes (MePOD-13, -17, -85, and -86) were found to be transcriptionally upregulated after multiple different treatments, suggesting that these genes are good candidates to target for cassava improvement. These findings increase our understanding of POD-mediated stress and hormone responses and storage root deterioration in cassava, laying a foundation for the genetic improvement of cassava.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/11/2730/ s1. Table S1. Characteristics of PODs in cassava. Table S2. The Ka/Ks ratios of duplicated POD genes in cassava. Table S3. The expression profiles (log2-based values) of the cassava POD genes in different tissues. Table S4. The expression profiles (log2-based fold changes) of the cassava POD genes after drought treatment. Table S5. The expression profiles (log2-based fold changes) of the cassava POD genes after harvest. Table S6. The expression profiles (log2-based values) of the cassava POD genes after various stress treatments. Table S7. Primers used in qRT-PCR analysis. Figure S1. Phylogenetic analyses of 91 MePODs, 2 AtPRXs, 3 OsPRXs, and 3 ZmPRXs.