Systematic Analysis of Cinnamyl Alcohol Dehydrogenase Family in Cassava and Validation of MeCAD13 and MeCAD28 in Lignin Synthesis and Postharvest Physiological Deterioration

1. Introduction

Cassava (Manihot esculenta Crantz) is a vital crop that plays a significant role in global food security because of its carbohydrates and essential nutrients, particularly in regions where it serves as a staple food source for millions of people [1]. Its global production was 330 million tons in 2022 (FAOSTAT) [2]. However, cassava faces challenges during postharvest handling and storage, with postharvest physiological deterioration (PPD) being a major concern that affects the quality and market value of tuberous roots [3]. Mechanical damages caused by collision, vibration, and extrusion during harvesting and storage are the main causes of PPD [4]. Then, a series of biochemical and physiological changes occurs, leading to typical phenomenon such as rapidly tissue deterioration, discoloration, and loss of firmness. Ultimately, nutritional quality decreases and susceptibility to microbial decay increases [5]. The main research on PPD primarily focuses on reactive oxygen species (ROS) elimination, programmed cell death [6,7,8,9,10]. Furthermore, wound healing can maintain the shelf life of tuberous roots by preventing water evaporation, nutrient outflow and pathogen infection [11,12,13]. The wound healing of cassava tuberous root is reportedly too late or inadequate to prevent the PPD response [14]. Therefore, deepening the research on cassava wound healing has become very necessary.

The phenylpropanoid metabolism pathway regulates wound healing by synthesizing phenolic compounds, lignin, and suberin [15]. For instance, activating phenylpropane metabolic pathway can promote wound healing in potato tubers [16]. Lignin is a phenolic heteropolymer in secondary cell walls, which plays a major role in mechanical support and defense against biotic and abiotic stresses, as well as affects the quality of fruits [17,18,19,20,21]. The biosynthesis of lignin involves many enzymes and corresponding genes, in which cinnamyl alcohol dehydrogenase (CAD) catalyzes the final step in monolignols biosynthesis [22]. Decreased lignin content has been reported in CAD down-regulation plants [23,24]. CADs are also involved in various stresses, such as drought and salinity stresses [25,26]. Until now, CAD gene families have been identified in many plants, including Arabidopsis [23], rice [27], poplar [28], mulberry [29], and pear [30]. However, the genome organization, gene structure, and expression profiling in PPD of MeCAD gene family in cassava is poorly understood.

In the present study, we identified the members of the MeCAD gene family in cassava, and then analyzed their phylogenetic relationships, gene structures, conserved motifs, synteny, and expression patterns with the PPD transcriptomic data. Additionally, the function of 2 MeCADs in lignin biosynthesis was verified by subjecting transgenic cassava to virus-induced gene silencing (VIGS). Our finding can provide a foundation for the in-depth functional analysis of MeCAD genes, as well as understanding the mechanisms underlying PPD regulation. Accordingly, effective postharvest management strategies can be developed to prolong the shelf life of tuberous roots and minimizing PPD losses.

2. Results

2.1. Identification and Phylogenetic Analysis of CAD in M. esculenta

A total of 36 CAD proteins were characterized from M. esculenta and denoted as MeCAD1 to MeCAD36, respectively. All MeCADs contained the CAD domain based on SMART tests and Pfam analysis (Figure S1). The MeCAD proteins’ lengths ranged from 353 (MeCAD28) to 435 (MeCAD17) amino acids, MW from 38.41 to 47.36 kDa, and pI from 5.50 to 8.22. Among the 36 MeCAD proteins, one was predicted to be peroxisomes protein (MeCAD10), five were located in the chloroplast (MeCAD17, 20, 22, 23 and 32), one was cytoskeleton protein (MeCAD18). The rest were localized in the cytoplasm. Detailed information including CAD name, gene accession, chromosome locus, protein length, MW, pI of all MeCAD proteins were shown in Table S1.

The CAD domain sequences of 36 MeCAD proteins, 12 OsCAD proteins and 9 AtCAD proteins were used to construct a phylogenetic tree through NJ method. Results indicated that the 36 MeCADs were classified into four groups (Figure 1). Group V contained the highest CAD members with 16 MeCADs, followed by group I with 14 MeCADs. Group IV and Group III had the least CAD members with only 3 MeCAD. However, Group II only including 8 OsCAD proteins. MeCAD13, MeCAD15 and MeCAD28 were clustered in group IV, including AtCAD4, AtCAD5 and OsCAD2, which were primary genes involved in lignin biosynthesis in Arabidopsis [31]. This finding suggested that CAD members in group IV of cassava con be the important genes involved in lignin biosynthesis in cassava.

2.2. Chromosomal Locations, Duplications, and Synteny Analysis of the MeCAD Genes

A total of 36 MeCAD were mapped on 12 chromosomes, and the distribution on each chromosome was uneven (Figure 2). For instance, chromosome 02 and 13 contained eight MeCAD genes, individually. however, chromosome 01, 05, 06, 07, 08, 10, 14 and 16 contained only one CAD genes on each chromosome. Interestingly, many MeCAD genes such as chromosome 02, 13 and 18 were clustered in a short distance. Six pairs (MeCAD2/MeCAD33, MeCAD13/MeCAD15, MeCAD15/MeCAD28, MeCAD16/MeCAD21, MeCAD18/MeCAD25 and MeCAD20/MeCAD23) of segmental duplications were identified in total (Figure 3 and Table S2). To detect the synteny of CAD genes, a collinearity analysis between M. esculenta, A. thaliana and O. sativa using MCScanX was performed. We found 7 paired collinearity relationships between 5 MeCAD and 4 AtCAD genes, as well as 1 paired collinearity relationships between 1 MeCAD and 1 OsCAD (Figure S2 and Table S3).

2.3. MeCAD Structure, Conserved Motifs and Cis-Acting Regulatory Element

The exon/intron organization and conserved motifs of all MeCAD genes were analyzed (Figure 4). The exon number ranged from 5 to 10, and 10 MeCAD genes contained 10 exons, the highest number of all genes. They included MeCAD1, MeCAD2, MeCAD9, MeCAD17, MeCAD31, MeCAD32, MeCAD33, MeCAD34, MeCAD35, MeCAD36. Meanwhile, MeCAD14 contained 7 exons, MeCAD22 contained 8 exons, MeCAD16 and MeCAD21 contained 9 exons, MeCAD10, MeCAD11, MeCAD12, MeCAD29, and MeCAD30 contained 6 exons, and the other 17 MeCAD genes contained 5 exons. Additionally, MeCAD genes in the same group had similar gene structure. Notably, members with high similarity in the same group shared a common motif composition. For instance, MeCAD13, MeCAD15 and MeCAD28 were found to contain 7 motifs. Thus, there three genes may have similar functions. MeCAD proteins contain 4-7 conserved motifs. However, only MeCAD12 contained 4 motif and 2 MeCAD proteins (MeCAD11, MeCAD22) contained five motifs. A total of 15 MeCAD proteins and 18 MeCAD proteins contain 6 motifs and 7 motifs, individually.

The cis-acting regulatory elements of each MeCAD promoter were characterized using PlantCARE and PLACE database. The most 10 cis-acting regulatory elements were founded in 36 MeCADs (Figure S3). The common regulatory elements such as CAAT-box and TATA-box were presented in all 36 MeCAD genes. We also identified cis-acting regulatory elements related to hormone responses, namely, MYB, ABRE and ARE, which were involved in abscisic acid (ABA) responsiveness. ERE, MeJA responsive elements. Interestingly, the promoter of 26 MeCAD genes contained elements related to wounding (WUN-motif), as well as other elements related to biotic and abiotic stress, such as drought (MBS), pathogen defense (W-box). In summary, MeCAD genes can be regulated by diverse hormone and different environment.

2.4. Expression Profile of MeCADs and Lignin Changes in Response to PPD

The expression levels of MeCADs in tuberous roots during PPD were analyzed by RNA-seq data. A total of 32 (88.89%) MeCAD genes exhibited differences in their expression levels in response to PPD (Figure 5A). Additionally, 27 MeCADs showed significant differences in compared with the control. Among them, 6 up-regulated genes (MeCAD9, MeCAD11, MeCAD13, MeCAD16, MeCAD26, MeCAD28) were verified by qRT-PCR in PPD samples. Thus, the expression of 6 MeCAD genes remarkably increased as PPD progressed (Figure 5B). However, the expression levels of MeCAD9, MeCAD28 at 5 d were less than those at 3 d, but still higher than those at 0 d. These results showed no significant difference with RNA-seq, where they were all up-regulated compared with those at 0 d.

Lignin biosynthesis was associated with PPD process, so we also detected the lignin content in tuberous roots during PPD. Results showed that lignin content significantly increased with PPD degree. The lignin content in tuberous roots was 67.45 μg/mg after harvest, whereas on the 1^st, 3^rd and 5^th after harvest, the values were 69.22, 70.16 and 72.22 μg/mg, respectively (Table 1).

2.5. MeCAD13 and MeCAD28 Were Necessary for Lignin Biosynthesis and Wound Response

Tissue-specific expression of MeCAD13 and MeCAD28 revealed that MeCAD13 was primarily expressed in tuberous roots and MeCAD28 was primarily expressed in fibrous roots (Figure S4). We also investigated the subcellular localization of MeCAD13 and MeCAD28 proteins in planta. Consistent with a previous study, GFP-fused MeCAD13 protein was localized in cytoplasm, and GFP-fused MeCAD28 protein was localized in cytoplasm even in tobacco leaves (Figure 6A), indicating their similar localization with prediction in Table S1. To study the in vivo roles of MeCADs in cassava lignin biosynthesis, two gene-silenced cassava lines were constructed (pCsCMV-CAD13 where MeCAD13 was silenced, and pCsCMV-CAD28 where MeCAD28 was silenced). The phenotypes between MeCADs-silenced lines and control plants are shown in Figure 6B. When the transcript levels of MeCAD13 and MeCAD28 were evaluated in MeCAD-silenced lines, the transcripts level of targeted MeCADs were drastically reduced compared with pCsCMV-NC (Figure 6C). We compared the lignin content of MeCADs-silenced lines in leaves with that of control plants, and found that lignin content significantly decreased in the leaves with MeCAD13- and MeCAD28-silenced lines (Table 1). In summary, the lignin content in pCsCMV-NC (control) leaves was 25.39 μg/mg, whereas the lignin content in the leaves of MeCAD28-silenced lines was 14.89 μg/mg, a decrease of 41.35% compared with the control. However, the lignin content in the leaves of MeCAD13-silenced lines decreased more significantly than in MeCAD28-silenced lines, reaching 60.06%. Therefore, MeCAD13 and MeCAD28 may be candidate genes involved in lignin biosynthesis, and the function of MeCAD13 in lignin synthesis may be superior to that of MeCAD28.

To understand the response of MeCAD13 and MeCAD28 to wound stress, the genes were transferred into yeast cells, and the growth of yeast cells was observed by simulating mechanical wound stress with upon systemin addition. As shown in Figure S5, the growth of pDR196-MeCAD13 and pDR196-MeCAD28 transgenic yeast strains was better than that of pDR196 empty vector after mechanical wound treatment. This indicated that the expression of MeCAD13 and MeCAD28 in yeast improved the tolerance of transgenic yeast cells to mechanical wound and enhanced their growth status.

3. Discussion

CAD is a major rate-limiting enzyme in lignin biosynthesis. It is related to stress responses, and is also closely associated with vegetative tissue and aging [32,33]. The CAD gene family has been identified and functionally analyzed genome-wide in many plant species, such as Arabidopsis [23], rice [27], wheat [34], sweet potato [17], sorghum [35], tea [36]. The present study aimed to identify and analyze CAD genes in cassava, as well as explore the candidate MeCADs genes for regulating lignin synthesis and PPD. A total of 36 MeCAD gene members were identified in cassava genome, and the number was higher than those in Arabidopsis and rice. The reason for the increase in MeCAD numbers may be the genome of “AM560” or the gene duplication events.

We subsequently analyzed and predicted gene functions through phylogenetic-tree analysis. Results revealed that 36 MeCADs can be classified into four subfamilies (Figure 1), MeCAD13/15/28 were clustered in cluster IV with AtCAD4/5 and OsCAD2. AtCAD4 and AtCAD5 in cluster IV played an important role in the biosynthesis of lignin in Arabidopsis thaliana [31], also OsCAD2 in Oryza sativa [37]. We speculated that the CADs of cluster IV in cassava may have a function similar to lignin synthesis. In our study, the lignin content significantly decreased in MeCAD13 silenced and MeCAD28 silenced leaves, especially in MeCAD13 silenced lines (Table 1). Therefore, MeCAD13 and MeCAD28 were the key genes involved in lignin synthesis in cassava. Furthermore, the homologous CAD genes may participate in lignin biosynthesis in one type of tissue during different developmental stage or different tissues [38]. For example, AtCAD4 was strongly expressed in flowers and roots, and AtCAD5 was expressed in lignified roots and strongly expressed in pathogen-infected tissues [39]. MeCAD13 was strongly expressed in tuberous roots, and MeCAD28 was strongly expressed in fibrous roots (Figure S4). MeCAD10/29/30 were clustered in cluster III with OsCAD1/4/6, consistent with a previous finding that OsCAD6 is not directly involved in lignin biosynthesis, but may be participated in lignans [27]. We speculated that MeCAD10/29/30 may be correlated with lignans synthesis. No MeCAD and AtCAD family members were in cluster II with 8 OsCADs, and cluster V included only 16 MeCADs, which was probably due to the methods and protein sequences. Notably, the function of MeCADs in these cluster requires further elucidation.

Gene duplication plays an essential role in species evolution and gene family evolution [40]. In this study, we identified six pairs of segmental duplications, while no tandem duplications in cassava (Figure 3 and Table S2). Tandem duplicated genes were relatively less, consistent with the results in tobacco [41]. Thus, segmental duplication was the main duplication model in cassava genome. A large number of hormone response elements (ABA, MeJA) and stress response elements (wound, drought) were founded in MeCAD gene promoter (Figure S3). It implied that MeCAD genes may be involved in various stress responses and can be induced by many hormones. For example, ABA has been shown to induce CADs gene expression in sweet potato after exposure to abiotic stress [17]. ABA-induced CmCADs expression is reportedly related to the upstream MYB response elements [32], whereas 32 MeCAD genes contained MYB elements in promoter regions. Nevertheless, this speculation should be demonstrated by further cloning and function analyses of MeCAD promoter. The hormone regulation on different MeCAD gene members, and the complex hormonal regulation of MeCAD genes under different stress conditions also warrant further investigation.

The phenylpropanoid pathway genes including CADs are closely related with biotic and abiotic stresses, including lignin deposition in secondary cell walls and defense-related compounds biosynthesis [42,43,44]. As mentioned above, the overexpression of MeCAD13 and MeCAD28 also enhanced the wound resistance of yeast (Figure S5). PPD is a kind of mechanical damage caused by wound, the transcriptional expression and enzyme activity of CAD continuously increased, and a net-like lignin layer formed in wound sites of SC9 tuberous roots [45]. A total of 32 MeCAD genes responded to PPD, whereas 27 MeCAD genes significantly changed, among which 6 were verified by qRT-PCR (Figure 5). These results were consistent with the changes in CAD enzyme activity in tuberous roots [6]. Lignin accumulation forms a physical barrier to limit pathogens invasion [46], thereby preventing the rapid postharvest deterioration of cassava roots [45]. Lignin content increased in cassava tuberous roots after PPD initiation, and significantly increased with prolonged storage time (Table 1), which may be due to the up-regulation expression of MeCAD genes. Wound healing involves cell wall reinforcement by lignin and other phenolic compounds, in which phenylpropane metabolism plays a vital role [15]. Phenylpropane metabolism was activated after cassava root injury indicated by an increase of phenylalanine ammonia-lyase (PAL) enzymatic activity, coupled with cinnamic acid 4-hydroxylase (C4H), 4-coumarate-CoA ligase (4CL) and CAD [45]. Moreover, MeCAD13 can interacted with MePOD12, delaying PPD occurrence through ROS elimination [6]. All these finding indicate that MeCADs are closely related with PPD regulation in cassava tuberous roots, whereas the specific functions of MeCADs in PPD remains to be demonstrated.

4. Materials and Methods

4.1. Plant Materials

The tuberous roots of South China 9 (SC9) were harvested after planting for 10 months. They were located in the National Cassava Germplasm Repository, Danzhou, Hainan, China. The tuberous roots then cultivated in a chamber at 26 °C and a 16/8 h photoperiod (day/night). After 1, 3 and 5 days, the tuberous roots were selected as materials, and PPD score was determined according to Rahmawati et al [47]. All samples were collected and frozen in liquid nitrogen immediately and stored at -80 °C for further analysis. For VIGS analysis, nine SC9 plants were inoculated with the recombinant VIGS vector after 25 days post inoculation (dpi). Leaves, petiole, stem, FEC, flowers, buds, fibrous roots, and tuberous roots of SC9 were used for tissue-specific analysis. Nicotiana benthamiana was used for subcellular localization.

4.2. Identification of MeCAD Genes in Cassava Genome

The whole genome sequences (v. 8.0) of cassava were downloaded from the phytozome database (https://phytozome-next.jgi.doe.gov/). The Hidden Markov Model (HMM) profile of CAD (PF00107 and PF08240) was retrieved from Pfam (http://pfam.sanger.ac.uk/). The HMMER program was used to search for CAD protein in cassava genome, and then all putative proteins were further filtered by e-Value < 0.001 and confirmed by Pfam and SMART database (http://smart.embl-heidelberg.de/). The molecular weight (MW) and theoretical isoelectric point (pI) of these proteins were predicted by ExPASy (https://www.expasy.org/). Finally, we named the proteins according to their locations on the chromosomes. Subcellular localization of MeCADs was predicted by WoLF PSORT (https://wolfpsort.hgc.jp/).

4.3. Chromosomal Mapping, Gene Structure, Conserved Motif, and Cis-Acting Regulatory Element Analysis

Chromosomal mapping of these MeCAD genes was constructed by MapChart program (v. 2.32). Gene structures were visualized using GSDS2.0 (http://gsds.gao-lab.org/). Additionally, the online software MEME (http://meme-suite.org/tools/meme) was used to identify the conserved motifs among all MeCAD genes. The number of motifs was set to 10. Then, 2000 bp upstream sequences of transcription start site ATG were extracted as the promoter sequence and screened using PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) before finally visualizing with GSDS2.0.

4.4. Phylogenetic Analysis, Gene Duplication, Multiple Alignments, and Synteny Analysis

Clustal W was performed for the multiple sequence-alignment analysis of CAD domain sequences. MEGA 11 was used to constructed the phylogenetic tree with the neighbor-joining (NJ) method. Clustal X (V.2.0) program was used to compare the coding sequences of repeated genes. Circos program was used to illustrate the relationships of duplicated genes. The synteny of CAD genes between Manihot esculenta, Arabidopsis thaliana and Oryza sativa was detected by MCScanX.

4.5. RNA-seq Analysis and qRT-PCR Verification

The tuberous roots of SC9 were storage for 0, 1, 3 and 5 days and used for transcriptome analysis, the platform for transcriptome sequencing was Illumina (HiSeq X-Ten), and three biological repeats were implemented for each sample. Heat maps were drawn using Mev software. Total RNA was extracted from 3 g of tuberous roots by using a RNAprep Pure Plant plus Kit (Tiangen, China). One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen, China) was used for first-strand cDNA synthesis. The qRT-PCR reactions were performed in 10 µL volume in a thermocycler (Thermo Fisher Scientific Inc., Göteborg, Sweden). The qRT-PCR primers were shown in Table S4. The reference gene of MeActin served as an internal reference. All experiments of each gene were performed in triplicate per sample. The formula 2^-△△Ct method was used to calculate the relative gene expressions.

4.6. Gene Cloning, Expression and Subcellular Localization of MeCAD13 and MeCAD28

Full-length MeCAD13 and MeCAD28 were amplified by PCR using the cDNA of SC9 tuberous roots as a template. They were then inserted into the transient expression vector pNC-Green-SubN to generate 35S::GFP-MeCAD13 and 35S::GFP-MeCAD28 recombinant vector, respectively [48]. Table S5 lists the primers used for subcellular localization and sequencing. Agrobacterium tumefaciens strain GV3101-psoup was co-transformed with the recombinant plasmid, and then transiently expressed in Nicotiana benthamiana leaves, as described by An et al. [40]. pNC-Green-SubN served as the positive control. The GFP signal (395-509 nm) were observed with a laser scanning confocal microscope (TCS SP8, Leica) after 3 days.

4.7. Virus-Induced Gene Silencing (VIGS) in Cassava and qRT-PCR Verification

For vector construction, 300 bp MeCAD13 and MeCAD28 DNA fragments were cloned into pCsCMV-NC as described by Tuo et al. [49]. The primers were listed in Table S5. The recombinant plasmid was co-transformed in A. tumefaciens strain GV3101-psoup, and cultured at 28 °C in chamber before resuspending in 10 mM MES, 10 mM MgCl₂ and 100 μM acetosyringone. The preparations were injected into the back of cassava leaves by using syringe and grown in greenhouse for 23 days. The expression level of the target gene in MeCAD13 and MeCAD28 silenced lines was detected, and the qRT-PCR primers are shown in Table S4.

4.8. Quantitative Analysis of Lignin

The Lignin content in tuberous roots and leaves of cassava were measured as described in Foster et al. [50], with minor modification. The main method involved using acetyl bromide to dissolve lignin in the sample, followed by determining the total lignin content through a colorimetric method. A multi-mode microplate reader (Multiskan GO, Thermo Fisher Scientific, USA) was used for the quantitative analysis of lignin content.

4.9. Stress of Yeast W303

The CDS of MeCAD13 and MeCAD28 was amplified (the specific primers are listed in Table S5), and then inserted into the pDR196 vector. The pDR196-MeCAD13 and pDR196-MeCAD28 construct were then transformed into Saccharomyces cerevisiae W303. After transformation, the yeast cells were streaked on SD-Ura medium and grown overnight in liquid SD-Ura medium. After the adjustment to optical density (OD) ₆₀₀=1.0, the cells were centrifuged briefly for 30 s, the supernatant was discarded, 20 nM systemin (wound stress) was added. Serial dilutions (OD₆₀₀ = 0.1, 0.01, 0.001) were plated on the SD-Ura medium. The plates were incubated at 30 °C for 3 days and photographed.

4.10. Statistical Analysis

Excel 2019 and Statistical Product and Service Solution program（version 20, SPSS Inc., Chicago, IL, USA）were used for all statistical analyses. One-way ANOVA (Tukey) was conducted for comparisons of the expression level and lignin content in cassava leaves.

5. Conclusions

A genome-wide analysis of MeCAD gene family was conducted with particular focus on their response to cassava PPD and lignin synthesis. A total of 36 MeCAD genes were identified and characterized, unevenly harbored in 12 chromosomes. Four subfamilies were clustered according to protein sequence. Lignin content in tuberous root increased during PPD, and 27 MeCAD genes may play a crucial role in PPD. Among them, MeCAD13 and MeCAD28 enhanced wound stress tolerance of transgenic yeast cells. Moreover, silencing MeCAD13 and MeCAD28 significantly decreased the lignin content in cassava leaves. Overall, the results provide new insight for cassava engineering programs, including improve PPD tolerance and increase lignin content.