*Article* **Combined De Novo Transcriptome and Metabolome Analysis of Common Bean Response to** *Fusarium oxysporum* **f. sp.** *phaseoli* **Infection**

**Limin Chen <sup>1</sup> , Quancong Wu 1,\*, Weimin He <sup>1</sup> , Tianjun He <sup>1</sup> , Qianqian Wu <sup>2</sup> and Yeminzi Miao <sup>1</sup>**


Received: 24 November 2019; Accepted: 11 December 2019; Published: 12 December 2019

**Abstract:** Molecular changes elicited by common bean (*Phaseolus vulgaris* L.) in response to *Fusarium oxysproum* f. sp. *Phaseoli* (FOP) remain elusive. We studied the changes in root metabolism during common bean–FOP interactions using a combined de novo transcriptome and metabolome approach. Our results demonstrated alterations of transcript levels and metabolite concentrations in common bean roots 24 h post infection as compared to control. The transcriptome and metabolome responses in common bean roots revealed significant changes in structural defense i.e., cell-wall loosening and weakening characterized by hyper accumulation of cell-wall loosening and degradation related transcripts. The levels of pathogenesis related genes were significantly higher upon FOP inoculation. Interestingly, we found the involvement of glycosylphosphatidylinositol- anchored proteins (GPI-APs) in signal transduction in response to FOP infection. Our results confirmed that hormones have strong role in signaling pathways i.e., salicylic acid, jasmonate, and ethylene pathways. FOP induced energy metabolism and nitrogen mobilization in infected common bean roots as compared to control. Importantly, the flavonoid biosynthesis pathway was the most significantly enriched pathway in response to FOP infection as revealed by the combined transcriptome and metabolome analysis. Overall, the observed modulations in the transcriptome and metabolome flux as outcome of several orchestrated molecular events are determinant of host's role in common bean–FOP interactions.

**Keywords:** common bean; *Fusarium oxysproum*; plant–pathogen interaction; transcriptome; metabolome

#### **1. Introduction**

Fusarium wilt, caused by *Fusarium oxysporum* f. sp. *phaseoli* (FOP), is a destructive soil-borne common bean (*Phaseolus vulgaris* L.) disease. Since its first identification in USA in 1929, this disease has been detected in all bean growing areas such as Africa, East Asia, Europe, Latin America, United Sates, and China [1,2]. Sever fusarium wilt epidemics have been reported in Heilongjiang and other common bean growing areas of China where beans follow vegetables [3]. High moisture, excessive irrigation, or poorly drained fields and a lack of rotation encourage the disease and FOP can persist in the soil indeterminately because of the production of chlamydospores and to the colonization of plant residues including roots of non-susceptible crops cultivated in rotation [3,4]. The majority of studies showed that invasion begins with the hyphal network development around root hairs followed by penetration and colonization of the epidermis and subsequently into the vascular tissues of the root. The colonization in vascular tissues leads it to the stem or the whole plant causing phloem blockage, internal stem discoloration, and total plant wilt. Infected plants display stunting, complete wilting, extensive chlorosis, and necrosis on the leaves [4,5].

Host-*F. oxysporum* pathosystems have been characterized in limited crops i.e., banana (*Musa paradisiaca*) [6,7], melon (*Cucumis melo*) [8,9], chickpea (*Cicer arietinum*) [10–12], cotton (*Gossypium hirsutum*) [13],tomato (*Lycopersicon esculentum*) [14], *Arabidopsis*[15,16], *Medicago truncatula* [17]. For each plant species, the respective *Fusarium* pathogens and a variety of defense mechanisms have been observed, including wound responses and hypersensitive reactions as well as many gene expression and metabolic changes. After the infection, plants recognize *Fusarium oxysporum* (FO) attack by understanding endogenous signals originating from the cell-wall through surveillance of cellular intactness. In this regard, many genes such as subtilin-like proteases, leucine rich-repeat proteins, proline-rice glycoproteins, cellulose synthases, and syntaxins are regulated as revealed in melon and banana [6,9,18]. In addition, constitutive enzymatic responses to FO infection appear to be important with changes in glutathione S-transferases, peroxidases, and phenylalanine ammonia lyase enzyme levels and activities being significant upon pathogen attack [19]. Changes also occur in the types and levels of cell-wall proteins, proteinase inhibitors, hydrolytic enzymes, and pathogenesis-related (PR) proteins and phytoalexin biosynthetic enzymes also appear to play important roles in FO defense [19–21]. Upregulation of genes involved in shikimate phenylpropanoid-lignin and cellulose synthesis pathways is possibly the reason of resistance in many cultivars where reduced spores are attached and resistance to FO is enhanced [18]. Identification of microbial surface derived molecules i.e., pathogen/microbe-associated molecular patterns (PAMPs/MAMPs) via pattern recognition receptors (PRRs) followed by binding of PAMPs to specific PRRs activates them and sends downstream signals to trigger broad spectrum immunity [22]. In terms of chemical defense, several genes inducing chitinases, xylem proteinases, β-1,3-glucanases thaumatin-like proteins, and peroxidases have been reported to be induced in melon against FOM infection [9]. Hypersensitive response signal molecules such as salicylic acid (SA), jasmonic acid (JA), plant growth regulating hormones, antioxidants, defense metabolites (polyphenols, phenolic acids, and flavonoids), and certain organic acids have been reported to be induced as an active plant defense [9,23]. Certain changes at the metabolic level such as altered activity of genes involved in sugar metabolism (sucrose synthase, invertase, and β-amylase) are also considered a plant response to FO. Sugars offer a dual function in plants as a nutrient as well as a signal to onset of disease, hence, these reactions are important when considering plant defense response [10,24]. The redox status of the intracellular (symplastic) and extracellular (apoplastic) spaces also change with Fusarium wilt infection [25].

In common bean, only a limited number of studies have been conducted to understand the mechanism and pathways involved in response to FOP infection. A recent study using cDNA amplified fragment length polymorphism technique reported transcript-derived fragments functionally characterized as metabolism, signal transduction, protein synthesis and processing, development and cytoskeletal organization, redox reaction, defense and stress response, transport proteins, and gene expression and RNA metabolism related genes/proteins [3]. However proteomic and metabolomics scale responses of common bean against FOP infection is poorly understood.

Unbiased modern high-throughput technologies such as combination of metabolomics and transcriptomics are required to improve our understanding of the plant–fungus interactions in common bean. Such combined approaches have recently resulted in the elucidation of different pathways in plants such as reprogramming of metabolites in chickpea roots in response to FO [12], understanding system responses to brown planthopper and rice stem borer infestation in rice [26,27]. Considering the amount of work done in model plants and other plant species infected with FO, it is important to understand the metabolic changes, transcriptional regulation, or physiological responses of bioactive and signaling compounds during infection of common bean with FOP. We aimed at identifying the differentially expressed genes (DEGs) and metabolites in common bean in response to FOP. The results showed that the transcriptome and metabolome response included structural defense, pathogen recognition receptors, and other components of innate immune system. Hormones played a crucial role in signaling pathways in common bean-FOP pathosystem. The most significantly enriched pathway as per metabolome results and confirmed by transcriptome analysis was the flavonoid biosynthesis

pathway. We observed a highly orchestrated response with significant modulation in various metabolic processes. The results described here thus improve our fundamental knowledge of molecular responses to the common bean–FOP interaction and potentially useful in designing strategies against wilt disease in common bean. was the flavonoid biosynthesis pathway. We observed a highly orchestrated response with significant modulation in various metabolic processes. The results described here thus improve our fundamental knowledge of molecular responses to the common bean–FOP interaction and potentially useful in designing strategies against wilt disease in common bean.

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 3 of 19

#### **2. Results 2. Results**

Fusarium wilt progression in common bean (Liyun No. 2) infected with FOP (FO; inoculated) and the control (CK; non-inoculated) was monitored by phenotypic screening at 4, 8, 12, 18, and 24 h post inoculation (hpi). The fusarium infection symptoms i.e., disease incidence index (severity rate 1–5), root length, fresh weight, and root volume were recorded for each time point for nine plants. After 12 hpi, the symptoms started to appear, and disease could be confirmed to a scale of 3. The disease incidence continued to affect roots and shoots. At 24 hpi, the disease incidence reached to a scale of 5 for all nine plants confirming the establishment of the disease. At this time point we confirmed that all treated seedlings at 24 hpi showed significant changes in recorded characteristics (Table S1). The control seedlings in contrast showed normal growth as compared to infected seedlings, remained healthy, and showed no fusarium wilt confirming the successful inoculation of treated plants (Figure 1). Fusarium wilt progression in common bean (Liyun No. 2) infected with FOP (FO; inoculated) and the control (CK; non-inoculated) was monitored by phenotypic screening at 4, 8, 12, 18, and 24 h post inoculation (hpi). The fusarium infection symptoms i.e., disease incidence index (severity rate 1–5), root length, fresh weight, and root volume were recorded for each time point for nine plants. After 12 hpi, the symptoms started to appear, and disease could be confirmed to a scale of 3. The disease incidence continued to affect roots and shoots. At 24 hpi, the disease incidence reached to a scale of 5 for all nine plants confirming the establishment of the disease. At this time point we confirmed that all treated seedlings at 24 hpi showed significant changes in recorded characteristics (Table S1). The control seedlings in contrast showed normal growth as compared to infected seedlings, remained healthy, and showed no fusarium wilt confirming the successful inoculation of treated plants (Figure 1).

**Figure 1.** *Fusarium oxysproum* f. sp. *phaseoli* infected (FO) and non-infected (CK) common bean roots and seedlings at 4, 8, 12, 18, and 24 h post infection. **Figure 1.** *Fusarium oxysproum* f. sp. *phaseoli* infected (FO) and non-infected (CK) common bean roots and seedlings at 4, 8, 12, 18, and 24 h post infection.

#### *2.1. RNA Sequencing and Identification of Differentially Expressed Genes 2.1. RNA Sequencing and Identification of Di*ff*erentially Expressed Genes*

The transcriptome of six samples of infected and non-infected common bean seedlings were sequenced using the Illumina HiSeq High-throughput sequencing platform. Illumina reads ranging from 43.65 to 48.71 million/sample (on average 46.25 million reads) were obtained from the six samples (Table S2). After filtering low quality reads and adapter sequences, a total of 39.57 Gb clean data was obtained. The clean data of each sample reached 5 Gb, and the Q30 base percentage was 93% or more. The clean data of the six samples of infected and non-infected common bean seedlings were de novo assembled as the reference gene set using the Trinity software and 136,238 unigenes, comprising 195,895,876 bp, were obtained, with a mean length of 1438 bp, a N50 of 2150 bp, and a N90 of 682 bp. The transcriptome of six samples of infected and non-infected common bean seedlings were sequenced using the Illumina HiSeq High-throughput sequencing platform. Illumina reads ranging from 43.65 to 48.71 million/sample (on average 46.25 million reads) were obtained from the six samples (Table S2). After filtering low quality reads and adapter sequences, a total of 39.57 Gb clean data was obtained. The clean data of each sample reached 5 Gb, and the Q30 base percentage was 93% or more. The clean data of the six samples of infected and non-infected common bean seedlings were de novo assembled as the reference gene set using the Trinity software and 136,238 unigenes, comprising 195,895,876 bp, were obtained, with a mean length of 1438 bp, a N50 of 2150 bp, and a N90 of 682 bp.

Functional annotation of all the unigenes was conducted, and a total of 80,409 (59.02%), 105,731 (77.61%), 71,638 (52.58%), 7,102,724 (75.4), 61,960 (45.48), 88,302 (64.81%), and 82,599 (60.63%) unigenes were annotated to the Kyoto encyclopedia of genes and genomes (KEGG), non-redundant (NR), Swiss-Prot, Trembl, euKaryotic Ortholog Groups (KOG), Gene ontology (GO), and Pfam database, respectively. Out of all unigenes, 106,777 (78.38%) were annotated in at least one database (Figure S1A). The functional information of homologous sequences in related species showed that Functional annotation of all the unigenes was conducted, and a total of 80,409 (59.02%), 105,731 (77.61%), 71,638 (52.58%), 7,102,724 (75.4), 61,960 (45.48), 88,302 (64.81%), and 82,599 (60.63%) unigenes were annotated to the Kyoto encyclopedia of genes and genomes (KEGG), non-redundant (NR), Swiss-Prot, Trembl, euKaryotic Ortholog Groups (KOG), Gene ontology (GO), and Pfam database, respectively. Out of all unigenes, 106,777 (78.38%) were annotated in at least one database (Figure S1A). The functional information of homologous sequences in related species showed that the transcript sequences had 68.71%, 9.88%, 2.96%, 2.39%, 2.01%, 1.43%, 1.13%, 0.89%, and 10.59% similarity with *P. vulgaris*, *Quercus suber*, *Vigna angularis*, *Vigna radiate* var. *radiate*, *Cajanus cajan*, *Glycine max*, *Ricinus*

*communis*, *Vigna angularis* var. *angularis*, and other genomes, respectively (Figure S1B). The GO annotation indicated 88,302 unigenes were categorized into 60 functional terms in three categories. Among them, genes associated with metabolic and cellular process in the category 'biological process'; cell, organelle, and cell part in the category 'cellular components'; and catalytic activity and transport activity in the category 'molecular function', were the most abundant (Figure S1C). The KEGG pathway database was used to analyze intracellular metabolic processes, and 80,409 unigenes were assigned to 144 KEGG pathways (Figure S1D). *Glycine max*, *Ricinus communis*, *Vigna angularis* var. *angularis*, and other genomes, respectively (Figure S1B). The GO annotation indicated 88,302 unigenes were categorized into 60 functional terms in three categories. Among them, genes associated with metabolic and cellular process in the category 'biological process'; cell, organelle, and cell part in the category 'cellular components'; and catalytic activity and transport activity in the category 'molecular function', were the most abundant (Figure S1C). The KEGG pathway database was used to analyze intracellular metabolic processes, and 80,409 unigenes were assigned to 144 KEGG pathways (Figure S1D).

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 4 of 19

the transcript sequences had 68.71%, 9.88%, 2.96%, 2.39%, 2.01%, 1.43%, 1.13%, 0.89%, and 10.59% similarity with *P. vulgaris*, *Quercus suber*, *Vigna angularis*, *Vigna radiate* var*. radiate*, *Cajanus cajan*,

For the FO-24 samples 84.51–87.26% and for the CK-24 samples 89.17–89.82% reads could be mapped to reference gene set (Table S3). Overall the Fragments Per Kilobase of Transcript per Million fragments mapped (FPKM) for FO-24 seedlings was higher for all three samples as compared to control (Figure 2a). Pearson correlations between FOP inoculated replicates ranged from 0.78 to 0.94 (Figure 2b). The CK replicates clustered together while the FOP inoculated replicates showed variability suggesting that inoculation within the treated plants differed (Figure 2c). For the FO-24 samples 84.51–87.26% and for the CK-24 samples 89.17–89.82% reads could be mapped to reference gene set (Table S3). Overall the Fragments Per Kilobase of Transcript per Million fragments mapped (FPKM) for FO-24 seedlings was higher for all three samples as compared to control (Figure 2A). Pearson correlations between FOP inoculated replicates ranged from 0.78 to 0.94 (Figure 2B). The CK replicates clustered together while the FOP inoculated replicates showed variability suggesting that inoculation within the treated plants differed (Figure 2C).

**Figure 2.** (**a**) Overall distribution of sample gene expression, (**b**) principle component analysis of expressed genes, and (**c**) Pearson correlations between CK-24 and FO-24 replicates. **Figure 2.** (**a**) Overall distribution of sample gene expression, (**b**) principle component analysis of expressed genes, and (**c**) Pearson correlations between CK-24 and FO-24 replicates.

#### *2.2. Differential Gene Expression Analysis Related to FOP Infection 2.2. Di*ff*erential Gene Expression Analysis Related to FOP Infection*

The screening conditions for the differentially expressed genes (DEG) were |log 2 Fold Change| ≥ 1, and False discovery rate (FDR) < 0.05. A total of 22,040 unigenes were expressed differentially with 8269 downregulated and 13,771 upregulated (Figure 3A). The screening conditions for the differentially expressed genes (DEG) were |log 2 Fold Change| ≥ 1, and False discovery rate (FDR) < 0.05. A total of 22,040 unigenes were expressed differentially with 8269 downregulated and 13,771 upregulated (Figure 3a).

Of the DEGs, 2780 downregulated genes were exclusively expressed in CK-24 while, 8231 upregulated DEGs were exclusively expressed in FO-24. These results show that the transcriptional changes are intense in FO infected common bean seedlings at 24 hpi. We further performed KEGG analysis to look at the key biological pathways involved in response to infection of FOP at 24 hpi. Spliceosome was the significantly enriched pathway with highest ratio of the number of differential genes annotated to this pathway to the number of annotated differential genes i.e., 314 out of 7204 genes in response to FOP infection 24 hpi. Other pathways such as endocytosis, RNA transport, mitogen-activated protein kinase (MAPK) signaling pathway-plant, mRNA surveillance pathway, amino sugar and nucleotide sugar metabolism pathway, and peroxisomes were significantly Of the DEGs, 2780 downregulated genes were exclusively expressed in CK-24 while, 8231 upregulated DEGs were exclusively expressed in FO-24. These results show that the transcriptional changes are intense in FO infected common bean seedlings at 24 hpi. We further performed KEGG analysis to look at the key biological pathways involved in response to infection of FOP at 24 hpi. Spliceosome was the significantly enriched pathway with highest ratio of the number of differential genes annotated to this pathway to the number of annotated differential genes i.e., 314 out of 7204 genes in response to FOP infection 24 hpi. Other pathways such as endocytosis, RNA transport, mitogen-activated protein kinase (MAPK) signaling pathway-plant, mRNA surveillance pathway, amino sugar and nucleotide sugar metabolism pathway, and peroxisomes were significantly enriched (Figure 3b).

enriched (Figure 3B). KEGG database (http://www.genome.jp/kegg/) was used to perform pathway mapping of the DEGs involved in common bean–FO interactions to facilitate the inspection of the plant gene networks. KEGG analysis revealed that unigenes were significantly enriched in various components involved in pathogen resistance mechanisms or signaling pathways (Figure 4).

To identify the most potential candidate genes related to resistance mechanism in common bean against FOP, we focused subsequent analysis on the 11,950 DEGs with fold change > 2 (Table S4).

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 5 of 19

**Figure 3.** (**a**) Differential gene MA map. The ordinate represents the log2 fold change value; the abscissa represents the average value of gene expression in the two samples; the red dot represents the upregulation of the gene expression, and the green dot represents the downregulation of the expression. Blue indicates no significant difference in gene expression. (**b**) Kyoto encyclopedia of genes and genomes (KEGG) enrichment scatter plot. The ordinate represents the KEGG pathway. The abscissa represents the Rich factor. The larger the Rich factor, the greater the enrichment. The larger the point, the greater the number of differential genes enriched in the pathway. The redder the color of the dots, the more significant the enrichment. KEGG database (http://www.genome.jp/kegg/) was used to perform pathway mapping of the **Figure 3.** (**a**) Differential gene MA map. The ordinate represents the log2 fold change value; the abscissa represents the average value of gene expression in the two samples; the red dot represents the upregulation of the gene expression, and the green dot represents the downregulation of the expression. Blue indicates no significant difference in gene expression. (**b**) Kyoto encyclopedia of genes and genomes (KEGG) enrichment scatter plot. The ordinate represents the KEGG pathway. The abscissa represents the Rich factor. The larger the Rich factor, the greater the enrichment. The larger the point, the greater the number of differential genes enriched in the pathway. The redder the color of the dots, the more significant the enrichment. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 6 of 19

**Figure 4.** KEGG orthology map (ko04626, plant–pathogen interaction) of common bean-*Fusarium oxysproum* f. sp. *Phaseoli* (FOP) pathosystem. For the treatment group, the red box labeled enzyme is associated with the upregulated gene, and the green box labeled enzyme is associated with the **Figure 4.** KEGG orthology map (ko04626, plant–pathogen interaction) of common bean-*Fusarium oxysproum*

downregulated gene. The blue labeled enzyme is related to both upregulation and downregulation. This pathway map is associated with DEGs. The enzymes are all marked with different colors.

2.2.1. Structural Defense

To identify the most potential candidate genes related to resistance mechanism in common bean

Plants recognize FOP attack for effective defense. Plants perceive an arsenal of endogenous signals originating from their own cell-wall by surveillance of cellular interactions. The plant cellwall is the first line of defense for plant cells and defines the primary strength of plants to restrict entry of pathogen to cell. We found 21 subtilisin-like protease and 46 DEGs encoding for leucine richrepeat (LRR) proteins were of the highly expressed proteins in common bean seedlings infected by FOP 24 hpi. Several genes involved in shikimate phenylpropanoid-lignin and cellulose biosynthesis pathways are reported to strengthen the cell-wall in resistant plants in response to FO infection [18]. f. sp. *Phaseoli* (FOP) pathosystem. For the treatment group, the red box labeled enzyme is associated with the upregulated gene, and the green box labeled enzyme is associated with the downregulated gene. The blue labeled enzyme is related to both upregulation and downregulation. This pathway map is associated with DEGs. The enzymes are all marked with different colors.

#### 2.2.1. Structural Defense

Plants recognize FOP attack for effective defense. Plants perceive an arsenal of endogenous signals originating from their own cell-wall by surveillance of cellular interactions. The plant cell-wall is the first line of defense for plant cells and defines the primary strength of plants to restrict entry of pathogen to cell. We found 21 subtilisin-like protease and 46 DEGs encoding for leucine rich-repeat (LRR) proteins were of the highly expressed proteins in common bean seedlings infected by FOP 24 hpi. Several genes involved in shikimate phenylpropanoid-lignin and cellulose biosynthesis pathways are reported to strengthen the cell-wall in resistant plants in response to FO infection [18]. We observed higher expression of 3-deoxy-d-arabino-heptulosonate-7-phosphate synthase, a lower expression of coumarate-CoA ligase gene which is reported to strengthen the cell-wall, indicating cell-wall weakening in response to FOP infection 24 hpi. A polyphenol oxidase gene had a very high expression in FOP infected plants as compared to the CK. Similarly, other genes such as UDP-glucuronic acid decarboxylase and cellulose synthases were downregulated in FOP infected plants. Other genes for cell-wall reinforcement have been reported to express during FO infection in tomato [9], we also noticed a shift in expression of four proline-rich glycoprotein, four hydroxyproline-rich glycoprotein, and 18 syntaxin genes. In addition, the bean response to FOP infection 24 hpi was characterized by hyper accumulation of transcripts coding for cell-wall degradation i.e., pectate lyases, pectin methylesterase inhibitors (PMEI), pectin methylesterases (PME), and Polygalacturonases (PG) (Table S4).

#### 2.2.2. Signaling

Plants recognize pathogen surface-derived molecules i.e., (PAMP/MAMP) via PRRs. This binding of PAMP to specific PRR activates these receptors and relays the signal downstream to convergent signaling pathways triggering broad-spectrum immunity. A total of 15 pathogenesis-related (PR) genes were differentially expressed. The fragments released in response to disruption of the first line of defense i.e., cell-wall (galacturonic acid-containing fragments) act as signals and mediate defense response by strengthening defensive barriers i.e., Chitin elicitor-binding protein (CEBiP) and chitin elicitor receptor kinase (CERK). In beans infected with FO-24, seven fungal elicitor immediate early-responsive genes showed higher expression as compared to CK-24. Further, 46 receptor kinases belonging to different gene families were expressed in FO-24 with a limited or zero expression in CK-24.

Involvement of glycosylphosphatidylinositol-anchored proteins (GPI-Aps) with extracellular ligands such as pathogen molecules as well as other ligands i.e., phytohormones, signaling polypeptides, leads to the phosphorylation of the intracellular kinase domain, which consequently activate cytoplasmic signaling components and switch on the response mechanisms. We also observed that glycosylphosphatidylinositol (GPI)-anchor biosynthesis pathway (K000563) was within significantly enriched pathways (Figure 3b). Apart from this, calmodulin (CaM) related DEGs were also induced suggesting the involvement of the CaM dependent signaling pathway.

The role of hormones in signaling pathways is well established which involves systematic acquired resistance. It has been suggested that resistance to FOP is mediated by SA, jasmonate (JA), and ethylene (ET) pathways [9,18]. In this regard 24 DEGs encoding for Ankyrin repeat containing protein genes were highly expressed in FO-24 beans. Genes responsive to ET are activated during early infection. We observed that four ET-insensitive protein 2 genes were highly responsive to FOP infection in FO-24 bean seedlings. This suggests that ET responsive genes might also be involved in latter infection stages of FOP and should be investigated further. Levels of core JA-signaling related genes i.e., four ZIM domain containing proteins, 12 TIFY, 26 ethylene-responsive transcription factor genes were highly expressed. Among these ethylene-responsive transcription factor-RAP2-7 and RAP2,

AP2-like ethylene-responsive transcription factor, ethylene-responsive transcription factor ERF113, ethylene-responsive transcription factor 1, ERF15 were uniquely expressed in FO-24 bean seedlings. There were nine superoxide dismutases expressed in FO-24 in higher fold change values as compared to CK-24 supporting the notion that JA levels are quite high in infected common bean seedlings. Further, the elevated levels of the transcripts of lipoxygenases, linoleate 13/9 S-lipoxygenases, and allene oxide cyclases suggested higher JA levels. However, the role of JA pathway in susceptibility and tolerance FO is still controversial. Almost a five-fold increase in transcript abundance of a defensin-like protein (*PHAVU\_005G071400g*) was observed. Homologs of this gene have been reported in *Arabidopsis* and are induced in response to FO infection. It is known that FO infection activates the transcription of auxin-related genes leading to a higher auxin biosynthesis. Three auxin influx career genes, nine auxin induced proteins, and one auxin responsive protein had higher expression in FO-24 seedlings as compared to CK-24. We observed the activation of the MAPK cascade. This pathway was also found to be significantly enriched (Figure 3b). *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 8 of 19

#### *2.3. Validation of DEGs by qRT-PCR 2.3. Validation of DEGs by qRT-PCR*

We validated the expression profiles of eight common bean genes of particular interest (Figure 5). The *Actin* gene was used as an internal control to standardize the data, and the amount of target gene transcript was normalized compared to the constitutive abundance of common bean *Actin* gene [3,28]. Among the common bean DEGs analyzed, five genes were upregulated in FO-24 as compared to CK-24 i.e., *PHAVU\_007G070400g*, *PHAVU\_004G134300g*, *PHAVU\_011G042100g*, *PHAVU\_008G232600g*, and *PHAVU\_007G185300g*. All these upregulated genes were characterized by similar trend in transcript accumulation at tested hpi supporting the RNA-Seq data. The other three genes i.e., *PHAVU\_003G141800g*, *PHAVU\_007G0495001g*, and *PHAVU\_007G236300g* were downregulated in response to FOP infection in FO-24 seedlings confirming the reliability of our RNA-Seq data. The three downregulated genes are cell wall related, a serine/threonine kinase activity related gene, and an Ankyrin repeat containing gene, respectively. We validated the expression profiles of eight common bean genes of particular interest (Figure 5). The *Actin* gene was used as an internal control to standardize the data, and the amount of target gene transcript was normalized compared to the constitutive abundance of common bean *Actin* gene [3,28]. Among the common bean DEGs analyzed, five genes were upregulated in FO-24 as compared to CK-24 i.e., *PHAVU\_007G070400g*, *PHAVU\_004G134300g*, *PHAVU\_011G042100g*, *PHAVU\_008G232600g*, and *PHAVU\_007G185300g*. All these upregulated genes were characterized by similar trend in transcript accumulation at tested hpi supporting the RNA-Seq data. The other three genes i.e., *PHAVU\_003G141800g*, *PHAVU\_007G0495001g*, and *PHAVU\_007G236300g* were downregulated in response to FOP infection in FO-24 seedlings confirming the reliability of our RNA-Seq data. The three downregulated genes are cell wall related, a serine/threonine kinase activity related gene, and an Ankyrin repeat containing gene, respectively.

**Figure 5.** qRT-PCR validation of the selected common bean differentially expressed genes (DEGs) in control (CK-24) and FOP infected plants (FO-24) 24 h post infection. **Figure 5.** qRT-PCR validation of the selected common bean differentially expressed genes (DEGs) in control (CK-24) and FOP infected plants (FO-24) 24 h post infection.

#### *2.4. Metabolite Profile 2.4. Metabolite Profile*

prospects for new bioactive compound discovery.

A combination of UPLC-MS/MS detection platform, self-built database, and multivariate statistical analysis was used to study the differences in metabolome between CK-24 and FO-24. It offers a platform to detect a great diversity of metabolites as previously reported in tomato [29], *Prunus mira* [30], and hulless barley [31–33]. In total, 754 metabolites were successfully detected in A combination of UPLC-MS/MS detection platform, self-built database, and multivariate statistical analysis was used to study the differences in metabolome between CK-24 and FO-24. It offers a platform to detect a great diversity of metabolites as previously reported in tomato [29], *Prunus mira* [30], and hulless barley [31–33]. In total, 754 metabolites were successfully detected in both sample types (Table S5).

both sample types (Table S5). The diverse set of detected molecules could be roughly grouped into 23 major classes, predominantly organic acids and derivatives, amino acids and their derivatives,

study have not yet been reported in the *Phaseolus vulgaris* metabolic network, our work offers

We compared the quantitative metabolic profiles between CK-24 and FO-24 roots in order to identify the compounds that differentially accumulated after infection. A series of pairwise OPLS-DA were applied to maximize the discrimination between experimental samples and to focus on metabolic variations significantly contributing to the resulting classifications. The differences between the control and infected groups in the OPLS-DA suggested that significant biochemical perturbation occurred in these samples (Figure S2). All significantly differentially expressed The diverse set of detected molecules could be roughly grouped into 23 major classes, predominantly organic acids and derivatives, amino acids and their derivatives, nucleotides and their derivatives, lipids, phenylpropanoids, and flavones. Collectively, phenolics were the major components of metabolome (flavanone, flavone, flavonoid, flavonol, isoflavone, polyphenol, anthocyanins, proanthocyanidins, phenolamides, phenylpropanoids) accounting for 1/3 of the total metabolites detected (Table S5). Moreover, since most of the identified metabolites in this study have not yet been reported in the *Phaseolus vulgaris* metabolic network, our work offers prospects for new bioactive compound discovery. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 9 of 19 metabolites (fold change ≥ 2 or ≤ 0.5), with the variable importance in the projection (VIP ≥ 1.0) between FO-24 and CK-24 roots are listed in Table S6 and Figure 6A. In total, 158 metabolites were differentially expressed, with 110 upregulated and 48 downregulated. We found that most of the significantly altered metabolites between FO-24 and CK-24 are

We compared the quantitative metabolic profiles between CK-24 and FO-24 roots in order to identify the compounds that differentially accumulated after infection. A series of pairwise OPLS-DA were applied to maximize the discrimination between experimental samples and to focus on metabolic variations significantly contributing to the resulting classifications. The differences between the control and infected groups in the OPLS-DA suggested that significant biochemical perturbation occurred in these samples (Figure S2). All significantly differentially expressed metabolites (fold change ≥ 2 or ≤ 0.5), with the variable importance in the projection (VIP ≥ 1.0) between FO-24 and CK-24 roots are listed in Table S6 and Figure 6a. In total, 158 metabolites were differentially expressed, with 110 upregulated and 48 downregulated. phenolic compounds (Table S6). The top 10 most differentially expressed metabolites are listed in Figure 6B. Among them, the upregulated metabolites include All-trans-13,14-dihydroretinol, Phillyroside, Isoeugenol, Quinone, N-Acetyl-L-tyrosine, D-Mannitol, E-3,4,5'-Trihydroxy-3' glucopyranosylstilbene, L-Carnitine, Prunetin, and L-Cystathionine. Similarly, the top 10 most downregulated metabolites include Luteolin 3',7-di-O-glucoside, 6-Gingerol, 5-Hydroxytryptophol, Peonidin O-malonylhexoside, 3,4,5-Trimethoxycinnamic acid, (3,4-Dimethoxyphenyl) acetic acid, Hinokitiol, 4-Hydroxycoumarin, N-Isovaleroylglycine, and Guanosine monophosphate. Differentially accumulated metabolites were functionally annotated using the KEGG database. It was observed that flavonoid biosynthesis, Glycerolipid metabolism, and Glycerophospholipid metabolism pathways were the most enriched (Figure 7).

**Figure 6.** (**a**) Heatmap hierarchical clustering of differentially expressed metabolites. Hierarchical trees were drawn based on differentially accumulated metabolites in CK-24 and FO-24. (**b**) Top 10 differentially accumulated metabolites in CK-24 and FO-24. **Figure 6.** (**a**) Heatmap hierarchical clustering of differentially expressed metabolites. Hierarchical trees were drawn based on differentially accumulated metabolites in CK-24 and FO-24. (**b**) Top 10 differentially accumulated metabolites in CK-24 and FO-24.

We found that most of the significantly altered metabolites between FO-24 and CK-24 are phenolic compounds (Table S6). The top 10 most differentially expressed metabolites are listed in Figure 6b. Among them, the upregulated metabolites include All-trans-13,14-dihydroretinol, Phillyroside, Isoeugenol, Quinone, N-Acetyl-L-tyrosine, D-Mannitol, E-3,4,5'-Trihydroxy-3'-glucopyranosylstilbene, L-Carnitine, Prunetin, and L-Cystathionine. Similarly, the top 10 most downregulated metabolites include Luteolin 3',7-di-O-glucoside, 6-Gingerol, 5-Hydroxytryptophol, Peonidin O-malonylhexoside, 3,4,5-Trimethoxycinnamic acid, (3,4-Dimethoxyphenyl) acetic acid, Hinokitiol, 4-Hydroxycoumarin, N-Isovaleroylglycine, and Guanosine monophosphate.

Differentially accumulated metabolites were functionally annotated using the KEGG database. It was observed that flavonoid biosynthesis, Glycerolipid metabolism, and Glycerophospholipid metabolism pathways were the most enriched (Figure 7).

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 10 of 19

**Figure 7.** Differential metabolite KEGG enrichment. The larger value indicates that the degree of enrichment is greater. The closer the *p-*value is to 0, the more significant the enrichment is. The size of the points in the graph represents the number of distinct significant metabolites enriched into the corresponding pathway. **Figure 7.** Differential metabolite KEGG enrichment. The larger value indicates that the degree of enrichment is greater. The closer the *p-*value is to 0, the more significant the enrichment is. The size of the points in the graph represents the number of distinct significant metabolites enriched into the corresponding pathway.

#### **3. Discussion 3. Discussion**

#### *3.1. Structural Defense in Response to FOP Infection in Common Bean 3.1. Structural Defense in Response to FOP Infection in Common Bean*

The present study reports the first combined de novo metabolome and RNA-seq analysis designed to describe the response of common bean infected with FOP. Previously, plant–FOP interaction has been demonstrated in many experiments in different crop plants [6–16]. In common bean it is known that colonization of FOP induces the defense responses in both a constitutive and inducible way. Pathogens are perceived by two different recognition systems that initiate patterntriggered immunity and effector triggered immunity in order to repel pathogens via induced defense responses [22]. Common bean, like other plants, employs cell-wall as the first barrier and defines primary or basic strength to encounter FOP infection. In this regard, the lower expression of UDPglucuronic acid decarboxylase and cellulose synthases, coumarate-CoA ligase and hyper accumulation of pectate lyases, PMEIs, PMEs, and PGs indicate that FOP has established itself at 24 hpi and cell-wall weakening in response to FOP infection has started. The higher expression of 3 deoxy-d-arabino-heptulosonate-7-phosphate synthase and polyphenol oxidase highlighted the timely recognition of FOP invasion and induction of the defense system [9,18]. These results confirmed that upon infection and establishment of FOP in common bean root tissues, the FOP secreted enzymes loosen and degrade the cell-wall i.e., pectin, cellulose, and hemicellulose [34]. The present study reports the first combined de novo metabolome and RNA-seq analysis designed to describe the response of common bean infected with FOP. Previously, plant–FOP interaction has been demonstrated in many experiments in different crop plants [6–16]. In common bean it is known that colonization of FOP induces the defense responses in both a constitutive and inducible way. Pathogens are perceived by two different recognition systems that initiate pattern-triggered immunity and effector triggered immunity in order to repel pathogens via induced defense responses [22]. Common bean, like other plants, employs cell-wall as the first barrier and defines primary or basic strength to encounter FOP infection. In this regard, the lower expression of UDP-glucuronic acid decarboxylase and cellulose synthases, coumarate-CoA ligase and hyper accumulation of pectate lyases, PMEIs, PMEs, and PGs indicate that FOP has established itself at 24 hpi and cell-wall weakening in response to FOP infection has started. The higher expression of 3-deoxy-d-arabino-heptulosonate-7-phosphate synthase and polyphenol oxidase highlighted the timely recognition of FOP invasion and induction of the defense system [9,18]. These results confirmed that upon infection and establishment of FOP in common bean root tissues, the FOP secreted enzymes loosen and degrade the cell-wall i.e., pectin, cellulose, and hemicellulose [34].

#### *3.2. Modulation of Defense Related Proteins in Common Bean 3.2. Modulation of Defense Related Proteins in Common Bean*

In order to trigger immunity, plants recognize the pathogen surface derived molecules PAMP/MAMP by employing PRRs which in turn activates the receptors and transcends the signals to different pathways and triggers broad spectrum immunity in plants. The PR proteins, among all defense related proteins, are induced and accumulated in response to FOP infection at 24 hpi and are considered an indispensable component of the innate immune system [35]. The 15 differentially expressed PRs in FO-24 suggest that in common bean, the innate immunity system is activated at this In order to trigger immunity, plants recognize the pathogen surface derived molecules PAMP/MAMP by employing PRRs which in turn activates the receptors and transcends the signals to different pathways and triggers broad spectrum immunity in plants. The PR proteins, among all defense related proteins, are induced and accumulated in response to FOP infection at 24 hpi and are considered an indispensable component of the innate immune system [35]. The 15 differentially expressed PRs in FO-24 suggest that in common bean, the innate immunity system is activated at this stage. Once, the FOP has established itself in bean tissues, the fragments released in response to disruption of the first line of defense i.e., cell-wall (galacturonic acid-containing fragments) act as signals and mediate defense response by strengthening defensive barriers i.e., CEBiP and CERK. In this regard, the high activation of fungal elicitor immediate early-responsive genes in FO-24 as compared to CK-24 confirms that such signals are received by common bean root tissues [18,22,36]. Involvement of GPI-APs with extracellular ligands such as pathogen molecules as well as other ligands i.e., phytohormones, signaling polypeptides, leads to the phosphorylation of the intracellular kinase domain, which consequently activate cytoplasmic signaling components and switch on the response mechanisms. In FO-24 bean roots, the GPI-anchor biosynthesis pathway was one of the significantly enriched pathways both in transcriptome as well metabolome analysis. This confirms that in common bean, GPI-APs is involved in signal transduction in response to FOP infection [37] (Table S3; Figure 3b).

#### *3.3. Crucial Role of Hormones in Signaling Pathways in Common Bean-FOP Pathosystem*

Role of hormones in signaling pathways is well established which involves systematic acquired resistance. It has been suggested that resistance to FOP is mediated by SA, JA, and ET pathways [9,18]. Previous reports on functional characterization of rice Ankyrin repeat containing proteins confirmed their role in defense against pathogen attack, particularly against *Magnaporthe oryzae* [38]. The contrasting expression of 24 Ankyrin repeat containing genes in FO-24 and CK-24 observed in our transcriptome clearly indicates that bean roots, under FOP infection, employ them as a defense response. Activation of the ET responsive genes at FO-24 is an important observation as previously it was known that some ET responsive proteins such as ET-insensitive protein-2 genes are involved in early infection in banana [18]. Hence their activation/expression at 24 hpi suggests that these genes might be involved in latter FOP infection stages and should be investigated further. The unique expression of JA-signaling related genes in FO-24 seedlings clearly indicates that hormone signaling pathways are involved defense responses in common bean against FOP infection (Table S4). However, the role of JA pathway in susceptibility and tolerance FO is still controversial [39,40]. The emerged response to FOP infection in our study confirms the involvement of ET/JA-dependent pathways together with the activation of TIFY, ET-responsive TFs against FOP infection. Similar response has been observed by Sebastiani et al. [9] in melon against FO infection. Our data suggests that FOP infection activates the transcription of auxin related genes exclusively in FO-24 roots which in turn increases the auxin biosynthesis and indicates direct involvement of auxin in common bean-FOP pathosystem. Together with signaling and structural defense responses, our results envisage that common bean employs a cascade of defense mechanisms including structural and signaling responses and that the auxin, ET, and JA are the main hormones involved in common bean-FOP pathosystem similar to what has been reported in tomato and banana [9,18,22].

#### *3.4. FOP Induced Energy Metabolism and Nitrogen Mobilization*

Obligate biotrophs depend on host metabolism for intake of nutrients, which is a measure of pathogenicity of the fungus within the host. Many plant defensive compounds are derived from amino acid precursors such as glucosinolates and secondary metabolites [41,42]. Our results confirm that in common bean roots infected with FOP, amino acids and their derivatives such as N-acetlyl-L-tyrosin, L-cystathionin, glutathione oxidized, glutathione reduced form, 5-aminovaleric acid, and Nα-Acetyl-L-arginine were upregulated play a crucial role in defense against FOP. In relation to this, it is important to relate the significantly enriched pathway observed in KEGG enrichment scatter plot i.e., amino sugar and nucleotide sugar metabolism (Figure 3b). This reveals that primary metabolites i.e., amino acids and sugars are playing a critical role in innate defense pathways. The activation and rapid accumulation of amino acids and sugars affect FOP susceptibility as observed in chickpea [12]. The upregulation of glutamate synthase, glutamate dehydrogenase, and aspargine synthase indicate the role of nitrogen mobilization. Significant upregulation of these genes correlated well with metabolite outcome (Tables S4 and S5).

### *3.5. FOP Resistance in Common Bean is Mediated by Flavonoid Biosynthesis Pathway*

It has been established that plants respond to pathogens by increased activation of the phenylpropanoid pathway leading towards biosynthesis of flavonoids, isoflavonoids, and phenolics. In our metabolome, the most significantly enriched pathway was flavonoid biosynthesis pathway (Figure 7). These compounds are regulated by chalcone synthase (CHS), chalcone isomerase (CHI), isoflavone synthase (IFS), isoflavone reductase (IFR), flavanone 3-hydroxylase (F3H), dihydroflavonol 4-reductase (DFR), anthocyanidin synthase (ANS), and leucoanthocyanidin reductase (LAR) genes and manipulated by MYB transcription factors [43]. As detailed in Table S5, these genes (except LAR gene(s)) were upregulated in FO-24 infected common bean roots. Similar results were revealed in chickpea infected with FO [12]. Previously, it was also known that fungal extract can induce these genes (CHI, IFS, and IFR) in Medicago cell suspension culture [44]. Consistently the accumulation of polyphenols, anthocyanins, flavanones, flavones, flavonoids, isoflavones, phenolamides, quinones, and terpenes identified by the metabolome profiles of CK-24 vs. FO-24 infected roots in the present study correlated well with transcriptomic data. Hence, the accumulation of flavonoid biosynthesis related genes and metabolites in FO-24 infected roots suggested their potential involvement in FOP resistance.

#### **4. Materials and Methods**

#### *4.1. Plant Growth and In Vivo Inoculations*

Seeds of common bean (Liyun No. 2) were obtained from the Lishui Institute of Agricultural Sciences, China. Seeds were sown in 15 cm diameter pots. The growing material filled in pots was sterile vermiculite and clay mixed in a 3:1 volume/volume ratio. The seedlings were allowed to grow under normal conditions i.e., 25 ◦C/18 ◦C day/night temperatures with a 16-h light/8-h dark photoperiod, and 60% humidity for 5 days. Plants were separated into two groups: control (CK) and *Fusarium oxysporum* f. sp. *phaseoli* (FOP) treated plants. Five individual seedlings were monitored at 4, 8, 12, and 24 h post infection. Then, 5-day-old seedlings were used for the infection at the fully expanded trifoliate leaves. The inoculum was prepared as reported earlier [3]. The control plants were supplemented with sterile ultrapure water. All the treated and control plants were evaluated for root quantitative traits such as root length, root volume, and fresh weight. As 24 h time point provided the best contrasting phenotype between CK and FO treated plants (Table S1, Figure 1), whole root samples from three individuals (biological replicates) in CK-24 and FO-24 were used for transcriptome and metabolome analysis. All treatments were grown in the same greenhouse with a 16 h light and 8 h dark cycle.

#### *4.2. RNA Extraction, cDNA Library Construction, and Transcriptome Sequencing*

Total RNAs were extracted using Spin Column Plant total RNA Purification Kit following the manufacturer0 s protocol (Sangon Biotech, Shanghai, China) [45]. Purity of the extracted RNAs was assessed on 1% agarose gels as well as by NanoPhotometer spectrophotometer (IMPLEN, Los Angeles, CA, USA). For RNA quantification we used a Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). Further, RNA integrity was assessed by the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Sequencing libraries was created using NEB Next Ultra RNA Library Prep Kit following manufacturer0 s instructions. The index codes were added to each sample. Briefly, the mRNA was purified from 3 µg total RNA of each of three replicate using poly-T oligo-attached magnetic beads. Subsequently, the fragmentation buffer was used to break the RNA into short fragments, and the short-fragment RNA was used as a template to synthesize the first strand cDNA with random hexamers, followed by buffer, dNTPs (dUTP, dATP, dGTP, and dCTP). The double-stranded cDNA was synthesized with DNA polymerase I, and then the double-stranded cDNA was purified using AMPure XP beads. The purified double-stranded cDNA was subjected to terminal repair, A tail was added, and the sequencing linker was ligated, and then AMPure XP beads were used for fragment size

selection, and finally PCR enrichment was performed to obtain a final cDNA library. Library quality was initially quantified using Qubit 2.0 using the 2100 to test the insert size of the library followed by accurately quantifying the effective concentration of the library (>2 nM) by Q-PCR. Finally, six paired-end cDNA libraries with an insert size of 300 bp were constructed for transcriptome sequencing and sequenced on Illumina HiSeq platform (Illumina Inc., San Diego, CA, USA) by Wuhan MetWare Biotechnology Co., Ltd. (www.metware.cn).

#### *4.3. De Novo Assembly, Functional Annotation, Classification, and Metabolic Pathway Analysis*

The clean reads were retrieved after trimming adapter sequences, removal of low quality (containing>50% bases with a Phred quality score<20) and reads with unknown nucleotides (more than 1% ambiguous residues N) using the FastQC tool (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/). GC content distribution check was performed. To stitch clean reads, Trinity was used (Version r20140717, [46]). For hierarchical clustering, Corset was used (https://code.google.com/p/corset-project/). The longest cluster sequence was obtained by clustering with Corset hierarchy as Unigene for subsequent analysis. The assembled unigenes were then aligned with various databases such as KEGG [47], GO [48], Clusters of Orthologous Groups (COG) [49], PfAM, Swissprot [50], egNOG [51], NR [52], KOG [53] using BLAST [54] with a threshold of E-value <sup>&</sup>lt; 1.0 <sup>×</sup> <sup>10</sup>−<sup>5</sup> .

The software KOBAS2.0 [55] was employed to get the unigene KEGG orthology. The analogs of the unigene amino acid sequences were searched against the Pfam database [56] using HMMER tool [57] with a threshold of E-value <sup>&</sup>lt; 1.0 <sup>×</sup> <sup>10</sup>−10. The sequenced reads were compared with the unigene library using Bowtie [58], and the level of expression was estimated in combination with RSEM [59]. The gene expression level was determined according to the FPKM.

#### *4.4. Di*ff*erential Expression and Enrichment Analysis*

The read count was normalized and EdgeR Bioconductor package [60] was used to determine the differential expression genes (DEGs) between CK-24 and FO-24 with the fold change > 2 [61] and FDR correction set at *p* < 0.01. GO enrichment analysis was performed using the topGO method based on the wallenius noncentral hypergeometric distribution with *p* < 0.05 [62]. KEGG pathway enrichment analysis of the DEGs was done using KOBAS2.0 [55]. The FDR correction was employed (*p* < 0.05) to reduce false positive prediction of enriched GO terms and KEGG pathways.

#### *4.5. Quantitative RT-PCR Analysis*

Eight DEGs, characterized by interesting expression profiles in response to FOP infection in FO-24 common bean plants were selected for qRT-PCR. First strand cDNAs were synthesized from 100 ng of total RNA using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystem, Carlsbad, CA, USA). Primers were designed using Primer3 Software (http://frodo.wi.mit.edu/primer3/; Table S7) and the specificity was checked by blasting their sequences in the NCBI database. The *Actin* constitutively expressed gene was used as the reference gene [3]. All qRT-PCR reactions were carried out on a Rotor-Gene 6000 machine (Qiagen, Shanghai, China) with the following thermal cycling profile: 50 ◦C for 2 min and 95 ◦C for 2 min, followed by 40 cycles at 95 ◦C for 3 s and 60 ◦C for 30 s. Melting curve analysis was performed to verify single product amplification with temperature ranging from 55 to 95 ◦C by increasing of 1 ◦C every step. All reactions were performed in a total volume of <sup>10</sup> <sup>µ</sup>L containing 30 ng of cDNA, 5 <sup>µ</sup>L 1 <sup>×</sup> SYBR® Select Master Mix (Applied Biosystem, Carlsbad, CA, USA), and 0.2 µL (20 µM) of each primer. For each sample, two biological replicates were analyzed in independent runs and a no-template control was included for each gene. Intra-assay variation was evaluated by performing all reactions in triplicate. The quantification cycle (Cq) was automatically determined using Rotor-Gene 6000 Series Software, version 1.7 as reported earlier [9].

#### *4.6. Widely Targeted Metabolomics*

The sample preparation, extract analysis, metabolite identification, and quantification were performed at Wuhan MetWare Biotechnology Co., Ltd. (www.metware.cn) following their standard procedures [30].

#### *4.7. Sample Preparation*

The vacuum freeze-dried root samples were crushed using a grinder (MM 400, Retsch, Haan, Germany) to a powder. A total of 100 mg powder was weighed and aliquots were extracted at 4 ◦C with 0.6 mL 70% aqueous methanol and vortexed six times to increase the extraction rate. After centrifuging at 10,000× *g* for 10 min, the supernatant was aspirated, and the sample was filtered through a microporous membrane (0.22 µm pore size) and stored in a sample bottle for UPLC-MS/MS analysis.

#### *4.8. Chromatographic Mass Spectrometry Acquisition Conditions*

The data acquisition instrument system included Ultra Performance Liquid Chromatography (UPLC) (Shim-pack UFLC SHIMADZU CBM30A, https://www.shimadzu.com.cn/) and tandem mass spectrometry (SHIMADZU Corp., Kyoto, Japan) (MS/MS) (Applied Biosystems 4500 QTRAP, http: //www.appliedbiosystems.com.cn/). The liquid phase conditions included (1) column: waters ACQUITY UPLC HSS T3 C18 1.8 µm, 2.1 mm × 100 mm; (2) mobile phase: phase A = ultrapure water (0.04% acetic acid was added), phase B = acetonitrile (0.04% acetic acid was added); (3) elution gradient: 0.00 min B = 5% in comparison, B was linearly increased to 95% in 10.00 min, and maintained at 95% 1 min, 11.00–11.10 min, B was reduced to 5%, and was 5% balanced to 14 min; (4) flow rate 0.35 mL/min; column temperature 40 ◦C; injection volume 4 µL. Whereas the mass spectrometry conditions were as following: the electrospray ionization (ESI) temperature was 550 ◦C, the mass spectrometry voltage was 5500 V, the curtain gas (CUR) was 30 psi, and the collision-induced dissociation (CAD) parameter was set high. In the triple quadrupole (QQQ), each ion pair was scanned for detection based on optimized decolusting potential (DP) and collision energy (CE) [63].

Based on the self-built database MWDB (metware database) at Wuhan MetWare Biotechnology Co., Ltd. (www.metware.cn), the material was characterized according to the secondary spectrum information. The isotope signal was removed during the analysis, and the repeated signals including K+ ions, Na+ ions, NH4+ ions, and fragment ions which are themselves other larger molecular weight substances were removed.

Metabolite quantification was performed using multiple reaction monitoring (MRM, as shown below) in triple quadrupole mass spectrometry. In the MRM mode, the fourth-stage rod first screens the precursor ions (parent ions) of the target substance, and excludes the ions corresponding to other molecular weight substances to initially eliminate the interference; the precursor ions break through the collision chamber to induce ionization and break to form a lot of fragment ions. The triple quadrupole filter is then used to select a desired feature fragment ion to eliminate non target ion interference, which makes the quantification more accurate and repeatable. After obtaining metabolite mass spectrometry data for different samples, peak area integration was performed on the mass spectral peaks of all the substances, and the mass spectral peaks of the same metabolite in different samples were integrated [64].

#### *4.9. Metabolomics Data Analysis*

Data matrices with the intensity of metabolite features under FOP and control conditions were uploaded to the Analyst 1.6.1 software (AB SCIEX, Ontario, Canada). For statistical analysis, missing values were assumed to be below the limits of detection, and these values were imputed with a minimum compound value [63]. The relative abundance of each metabolite was log transformed before analysis to meet normality. A Dunnett's test was used to compare the abundance of each metabolite between control and FOP. False discovery rate was used for controlling multiple testing. The supervised multivariate method, partial least squares-discriminant analysis (PLS-DA), was used to maximize the metabolome difference between the control and FOP treated samples. The relative

importance of each metabolite to the PLS-DA model was checked using a parameter called the variable importance in projection (VIP). Metabolites with VIP > 1.0 were considered as differential metabolites for group discrimination. Principal Component Analysis (PCA), Hierarchical Cluster Analysis (HCA), and KEGG pathway analysis were performed in R software (www.r-project.org).

### **5. Conclusions**

In the present study the whole transcriptome and metabolome of common bean infected by FOP 24 hpi were characterized. The differences in terms of DEGs between the inoculated and non-inoculated common bean showed that nitrogen metabolism and energy metabolism participated in defense response to FOP infection. Flavonoid pathway was the main defense response in common bean. Transcriptome analysis showed that the spliceosome, RNA transport, ribosome biogenesis in eukaryotes, proteasome, and phenylalanine metabolism were the top five significantly enriched pathways. Cell-wall related genes proved to be the first response to FOP attack and started a cascade of signaling leading to accumulation of cell wall degradation related transcripts. PAMP/MAMP, PRRs, and PRs were being regulated in response to FOP infection suggesting triggering of immunity in common bean. Activation of systematic acquired resistance was also observed in our study where the role of hormones in the signaling pathway was observed. These results demonstrate the common bean in response to FOP utilizes different and effective defense pathways comprising of a complex resistance network of structural, signaling, and chemical responses. Further investigations will be focused on functional validation and mapping of the DEGs, which could represent a helpful tool for developing common bean resistant varieties toward FOP.

**Supplementary Materials:** Supplementary files can be found at http://www.mdpi.com/1422-0067/20/24/6278/s1. Supplementary Table S1. Comparison of means for root quantitative data. Supplementary Table S2. Summary of sequencing output statistics. Supplementary Table S3. Mapping results of RNA-Seq data from inoculated and non-inoculated common bean seedlings. Supplementary Table S4. Complete list of differentially expressed genes regulated by infection of FOP at 24 hpi in control and infected roots. Supplementary Table S5. List of all detected metabolites. Supplementary Table S6. List of differentially expressed metabolites. Supplementary Table S7. List of primers used for qRT-PCR analysis. Supplementary Figure S1. Differential expression gene function annotation and enrichment analysis. Supplementary Figure S2. Orthogonal partial least squares discriminant analysis.

**Author Contributions:** Conceptualization, L.C., and Q.W.; methodology, L.C.; validation, W.H., T.H., and L.C.; formal analysis, Y.M.; data curation, L.C.; writing—original draft preparation, L.C., Q.W.; writing—review and editing, Q.W., Y.M.; visualization, T.H.; supervision, Q.W.; project administration, W.H.; funding acquisition, Q.W.

**Funding:** This research was funded by Zhejiang Natural Science Foundation, grant number No. LGN18C140004.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Availability of Data:** The data for RNA-sequencing is available at National Center for Biotechnology Information (NCBI) with accession number PRJNA589754.

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Molecular Analysis of UV-C Induced Resveratrol Accumulation in** *Polygonum cuspidatum* **Leaves**

## **Zhongyu Liu, Junxiong Xu, Xiang Wu, Yanyan Wang, Yanli Lin, Duanyang Wu, Hongjie Zhang and Jianbing Qin \***

College of Life Science, Yangtze University, Jingzhou 434025, China; yzrs91@163.com (Z.L.); 2017714456@yangtzeu.edu.cn (J.X.); xufyz20016@126.com (X.W.); llqww17@126.com (Y.W.); linyanli1998@foxmail.com (Y.L.); kdysuperme@163.com (D.W.); zhk18672558067@163.com (H.Z.) **\*** Correspondence: yzpc19@126.com; Tel.: +186-7256-3466

Received: 11 November 2019; Accepted: 4 December 2019; Published: 7 December 2019

**Abstract:** Resveratrol is one of the most studied plant secondary metabolites owing to its numerous health benefits. It is accumulated in some plants following biotic and abiotic stress pressures, including UV-C irradiation. *Polygonum cuspidatum* represents the major natural source of concentrated resveratrol but the underlying mechanisms as well as the effects of UV-C irradiation on resveratrol content have not yet been documented. Herein, we found that UV-C irradiation significantly increased by 2.6-fold and 1.6-fold the resveratrol content in irradiated leaf samples followed by a dark incubation for 6 h and 12 h, respectively, compared to the untreated samples. De novo transcriptome sequencing and assembly resulted into 165,013 unigenes with 98 unigenes mapped to the resveratrol biosynthetic pathway. Differential expression analysis showed that *P. cuspidatum* strongly induced the genes directly involved in the resveratrol synthesis, including phenylalanine ammonia-lyase, cinnamic acid 4-hydroxylase, 4-coumarate-CoA ligase and stilbene synthase (*STS*) genes, while strongly decreased the chalcone synthase (*CHS*) genes after exposure to UV-C. Since CHS and STS share the same substrate, *P. cuspidatum* tends to preferentially divert the substrate to the resveratrol synthesis pathway under UV-C treatment. We identified several members of the MYB, bHLH and ERF families as potential regulators of the resveratrol biosynthesis genes.

**Keywords:** regulation; RNA-seq; abiotic stress; biosynthesis pathway; chalcones; stilbenes

#### **1. Introduction**

*Polygonum cuspidatum* Sieb. et Zucc. is a member of the buckwheat family (Polygonaceae). It is a tall and resilient herbaceous perennial with woody rhizomes [1], native to East Asia in countries such as Korea, Japan and China. Although this plant has been recognized as an invasive species in Europe and North America [2,3], *P. cuspidatum* has an extraordinary value in phytotherapy. In China, the roots of *P. cuspidatum* have been long employed in traditional medicine to combat cough, hepatitis, infection, arthralgia, tumors, bronchitis, jaundice, bleeding, amenorrhea and hypertension [4–6]. Analytic investigations of the major health-promoting molecules of *P. cuspidatum* roots have revealed the presence of several lead compounds belonging to the group of polyphenols [7,8]. Distinctly, one of the most important compounds in *P*. *cuspidatum* roots, which has drawn growing interest on this species, is resveratrol, a molecule with proven anti-inflammatory and antioxidant activity [9,10].

The "French paradox", a curious observation referring to the low level of coronary heart disease in France despite high intake of dietary cholesterol and saturated fat [11], has been linked to the high consumption of red wine containing resveratrol. Resveratrol (3,5,40-trihydroxy-trans-stilbene) is a naturally occurring stilbene metabolite found in grapes, berries, nuts and other plants such as *P. cuspidatum*. Over the past decades, extensive researches have been carried out on the physiological functions of resveratrol in human and it has been suggested that resveratrol is a potential remedy for a range of diseases, including heart disease, diabetes, cancer and Alzheimer's disease [12]. The compound was first discovered in *P. cuspidatum* [13], which is the most important concentrated source of free or glycosylated resveratrol. Therefore, it is expected that this species should be a model plant to study resveratrol biosynthesis and engineering in plant. Surprisingly, very limited advances in these fields have been made so far in *P. cuspidatum* [14], contrasting with the extensive knowledge generated in grape (reviewed by Hasan and Bae, [15]). In fact, the lack of genomic resources for *P*. *cuspidatum* hinders genetic discoveries. Particularly, the key genes and molecular mechanisms underlying the striking accumulation of resveratrol in this species are still unknown. On the opposite, early genome sequencing of grape [16] has triggered much research on the biosynthesis of stilbenes. The biosynthetic pathway of resveratrol has been well characterized and involves four key enzymes, including phenylalanine ammonia lyase (PAL), cinnamic acid 4-hydroxylase (C4H), 4-coumarate: CoA ligase (4CL) and stilbene synthase (STS) [17]. p-coumaroyl-CoA is a product of PAL, which is abundant in plants and used as a precursor for the synthesis of both resveratrol and chalcone. Therefore, in stilbene-synthesizing plants, STS competes with chalcone synthase (CHS) for the synthesis of resveratrol [18].

Resveratrol is basically a defense compound (phytoalexin) in plants but it is produced in very small amounts [19]. Therefore, how to induce a strong accumulation of resveratrol in plant organs in order to satisfy the increasing global demand of resveratrol has become one of the hot topics in secondary metabolite research. Manipulation of different synthetic enzymes and the identification of their regulator genes such as *MYB14*, *MYB15* and *WRKY8* [20,21] provide currently prospects to increase resveratrol production in planta. Moreover, it has been demonstrated that biotic and abiotic factors, including fungal attack, UV-C irradiation, jasmonic acid, salicylic acid, H2O<sup>2</sup> and AlCl3, can induce the accumulation of resveratrol in grape [22–28].

Hao et al. [14] previously identified 54 unigenes predicted to participate in the resveratrol biosynthesis pathway in *P. cuspidatum* roots. However, the mechanism leading to the high resveratrol accumulation in this species and the responses of these genes to abiotic factors such as UV-C irradiation remain open to study. Here, we investigated the effect of UV-C irradiation on leaf resveratrol content in *P. cuspidatum* and further de novo sequenced the transcriptome, offering a novel insight into the UV-C induced biosynthesis of resveratrol in plants.

#### **2. Results**

#### *2.1. E*ff*ect of UV-C Irradiation on Leaf Resveratrol Content in Polygonum cuspidatum*

*Polygonum cuspidatum* leaves were treated with UV-C light for 10 min and then incubated in the dark for 6 h (PC6H) and 12 h (PC12H). We quantified the leaf resveratrol accumulation in UV-C treated leaves and untreated leaves (PC). As shown in Figure 1, we observed a significant difference in the resveratrol contents of UV-C treated leaves as compared to untreated leaves (*p* < 0.001) and between the different incubation times of treated leaves (*p* < 0.001). UV-C significantly increased the leaf resveratrol contents and 6 h incubation time after UV-C treatment showed the highest accumulation of resveratrol in the leaves of *P. cuspidatum*.

**2. Results** 

synthesis of resveratrol [18].

induce the accumulation of resveratrol in grape [22–28].

novel insight into the UV-C induced biosynthesis of resveratrol in plants.

*2.1. Effect of UV-C Irradiation on Leaf Resveratrol Content in Polygonum cuspidatum* 

a range of diseases, including heart disease, diabetes, cancer and Alzheimer's disease [12]. The compound was first discovered in *P. cuspidatum* [13], which is the most important concentrated source of free or glycosylated resveratrol. Therefore, it is expected that this species should be a model plant to study resveratrol biosynthesis and engineering in plant. Surprisingly, very limited advances in these fields have been made so far in *P. cuspidatum* [14], contrasting with the extensive knowledge generated in grape (reviewed by Hasan and Bae, [15]). In fact, the lack of genomic resources for *P*. *cuspidatum* hinders genetic discoveries. Particularly, the key genes and molecular mechanisms underlying the striking accumulation of resveratrol in this species are still unknown. On the opposite, early genome sequencing of grape [16] has triggered much research on the biosynthesis of stilbenes. The biosynthetic pathway of resveratrol has been well characterized and involves four key enzymes, including phenylalanine ammonia lyase (PAL), cinnamic acid 4-hydroxylase (C4H), 4-coumarate: CoA ligase (4CL) and stilbene synthase (STS) [17]. p-coumaroyl-CoA is a product of PAL, which is abundant in plants and used as a precursor for the synthesis of both resveratrol and chalcone. Therefore, in stilbene-synthesizing plants, STS competes with chalcone synthase (CHS) for the

Resveratrol is basically a defense compound (phytoalexin) in plants but it is produced in very small amounts [19]. Therefore, how to induce a strong accumulation of resveratrol in plant organs in order to satisfy the increasing global demand of resveratrol has become one of the hot topics in secondary metabolite research. Manipulation of different synthetic enzymes and the identification of their regulator genes such as *MYB14*, *MYB15* and *WRKY8* [20,21] provide currently prospects to increase resveratrol production in planta. Moreover, it has been demonstrated that biotic and abiotic factors, including fungal attack, UV-C irradiation, jasmonic acid, salicylic acid, H2O2 and AlCl3, can

Hao et al. [14] previously identified 54 unigenes predicted to participate in the resveratrol biosynthesis pathway in *P. cuspidatum* roots. However, the mechanism leading to the high resveratrol accumulation in this species and the responses of these genes to abiotic factors such as UV-C irradiation remain open to study. Here, we investigated the effect of UV-C irradiation on leaf resveratrol content in *P. cuspidatum* and further de novo sequenced the transcriptome, offering a

*Polygonum cuspidatum* leaves were treated with UV-C light for 10 min and then incubated in the dark for 6 h (PC6H) and 12 h (PC12H). We quantified the leaf resveratrol accumulation in UV-C treated leaves and untreated leaves (PC). As shown in Figure 1, we observed a significant difference in the resveratrol contents of UV-C treated leaves as compared to untreated leaves (*p* < 0.001) and

**Figure 1.** Effect of UV-C irradiation and incubation times on the resveratrol contents in *Polygonum*.

#### *2.2. De Novo Transcriptome Sequencing and Unigene Assembly*

accumulation of resveratrol in the leaves of *P. cuspidatum*.

In order to elucidate the molecular mechanism underlying the UV-C induced resveratrol accumulation in *Polygonum cuspidatum* leaves, we synthesized nine cDNA libraries from leaf samples of PC12H, PC6H, PC and generated de novo RNA-sequencing data.

A total of 73.63 Gb raw read data was generated. After cleaning, 98.7% of the reads were kept as clean reads with 93.5% of bases scoring Q30 and above (Table 1). The assembly was performed using the Trinity software and a total of 165,013 unigenes were obtained with N50 length about 1744 bp (Table 2). The unigene lengths ranged from 200 to 17,269 bp with a significant number of genes having length of 200–300 bp (Figure 2).

mature leaves. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively. The letters above the bar represent significant difference between samples. **Total Clean Clean Reads Clean Reads Detected**

**Table 1.** Overview of the transcriptome sequencing dataset and quality check. *Polygonum cuspidatum*


**Table 2.** Statistics of the unigene assembly results.


**Figure 1.** Effect of UV-C irradiation and incubation times on the resveratrol contents in *Polygonum*.

In order to elucidate the molecular mechanism underlying the UV-C induced resveratrol accumulation in *Polygonum cuspidatum* leaves, we synthesized nine cDNA libraries from leaf samples

A total of 73.63 Gb raw read data was generated. After cleaning, 98.7% of the reads were kept as clean reads with 93.5% of bases scoring Q30 and above (Table 1). The assembly was performed using the Trinity software and a total of 165,013 unigenes were obtained with N50 length about 1744 bp (Table 2). The unigene lengths ranged from 200 to 17,269 bp with a significant number of genes having

**Table 1.** Overview of the transcriptome sequencing dataset and quality check. *Polygonum cuspidatum* mature leaves. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively. The letters above the bar represent significant difference between samples.

**Type Read Length Total Clean Reads Clean Reads Q20 (%) Clean Reads Q30 (%) Detected Genes**  P12H1 150 55,184,328 97.81 93.76 122,534 P12H2 150 57,685,708 97.47 93.21 123,131 P12H3 150 60,727,642 97.83 93.85 124,005 P6H1 150 51,083,608 97.66 93.41 122,200 P6H2 150 53,076,416 97.64 93.36 120,570 P6H3 150 48,587,712 97.74 93.56 119,464 PC1 150 53,966,128 97.52 93.14 123,390 PC2 150 51,267,598 97.5 93.19 117,489 PC3 150 52,633,030 97.77 93.69 125,163

**Table 2.** Statistics of the unigene assembly results.

All 165,013 200 17,269 1,102.73 1,744 489 0.4216

*2.2. De Novo Transcriptome Sequencing and Unigene Assembly* 

length of 200–300 bp (Figure 2).

of PC12H, PC6H, PC and generated de novo RNA-sequencing data.

 **Figure 2.** Unigene length distribution. The numbers above the bar represent the number of unigenes with a sequence length comprised in the classes of size displayed in the x-axis. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 4 of 16

We performed the functional annotation of the unigenes in various databases, including NR, NT, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG) and Gene Ontology (GO), which resulted in 104,903 annotated unigenes (Figure 3A) with 32,448 unigenes successfully annotated based on all the five databases (Figure 3B). Overall, these unigenes contributed to various GO terms of which the most represented ones were binding, catalytic activity and transporter activity (molecular function). They were mainly located in the cell, cell part and membrane (cell component) and involved in the cellular process, metabolic process, etc. (biological process; Figure S1). We further analyzed the genes encoding for transcription factors within the annotated unigenes. As shown in Table 3, we identified 5526 TFs classified into 58 different TF families with bHLH, MYB, bZIP and ERF as the most abundant TFs. In addition, we searched for all genes belonging to the resveratrol biosynthetic pathway and identified 98 unigenes including, 26 phenylalanine ammonia-lyase (PAL), 15 cinnamic acid 4-hydroxylase (C4H), 20 4-coumarate-CoA ligase (4CL), 4 stilbene synthase (STS) and 33 chalcone synthase (CHS; Table S1). **Figure 2.** Unigene length distribution. The numbers above the bar represent the number of unigenes with a sequence length comprised in the classes of size displayed in the x-axis. We performed the functional annotation of the unigenes in various databases, including NR, NT, Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG) and Gene Ontology (GO), which resulted in 104,903 annotated unigenes (Figure 3A) with 32,448 unigenes successfully annotated based on all the five databases (Figure 3B). Overall, these unigenes contributed to various GO terms of which the most represented ones were binding, catalytic activity and transporter activity (molecular function). They were mainly located in the cell, cell part and membrane (cell component) and involved in the cellular process, metabolic process, etc. (biological process; Figure S1). We further analyzed the genes encoding for transcription factors within the annotated unigenes. As shown in Table 3, we identified 5526 TFs classified into 58 different TF families with bHLH, MYB, bZIP and ERF as the most abundant TFs. In addition, we searched for all genes belonging to the resveratrol biosynthetic pathway and identified 98 unigenes including, 26 phenylalanine ammonia-lyase (PAL), 15 cinnamic acid 4-hydroxylase (C4H), 20 4-coumarate-CoA ligase (4CL), 4 stilbene synthase (STS) and 33 chalcone synthase (CHS; Table S1).

**Figure 3.** Unigene functional annotation in *P. cuspidatum*. (**A**) The annotated gene number in different databases and (**B**) Venn diagram showing the number of shared and unique annotated genes in the different databases. **Figure 3.** Unigene functional annotation in *P. cuspidatum*. (**A**) The annotated gene number in different databases and (**B**) Venn diagram showing the number of shared and unique annotated genes in the different databases.

**Table 3.** Statistics of genes encoding for transcription factors. **TF Count TF Count TF Count**  LFY 3 M-type 5 HB-PHD 17

HRT-like 3 MYB 402 SRS 15 HB-other 45 LSD 14 NF-X1 8 FAR1 67 WOX 15 Nin-like 43 ERF 341 C2H2 204 MYB\_related 272 Whirly 11 CAMTA 23 VOZ 18 ARR-B 57 GRAS 109 G2-like 157 HSF 90 GRF 32 NF-YA 56 NZZ/SPL 8 EIL 36 Dof 119


**Table 3.** Statistics of genes encoding for transcription factors.

Next, the clean reads data were mapped to the assembled unigene libraries and the statistics of mapping results are presented in Table S2. A total of 157,665 genes were expressed with the number of fragments per kilobase of exon per million fragments mapped (FPKM) values raging from 0.02 to 18,562.92 (Figure 4A). Principal component analysis (PCA) of the samples based on FPKM showed that all the biological replicates clustered together, suggesting a high reliability of our RNA-sequencing data (Figure 4B). Furthermore, the PCA clearly distinguished the control and the UV-C treated leaf samples, indicating that a large number of genes were differentially expressed after exposure to UV-C irradiation. Next, the clean reads data were mapped to the assembled unigene libraries and the statistics of mapping results are presented in Table S2. A total of 157,665 genes were expressed with the number of fragments per kilobase of exon per million fragments mapped (FPKM) values raging from 0.02 to 18,562.92 (Figure 4A). Principal component analysis (PCA) of the samples based on FPKM showed that all the biological replicates clustered together, suggesting a high reliability of our RNAsequencing data (Figure 4B). Furthermore, the PCA clearly distinguished the control and the UV-C treated leaf samples, indicating that a large number of genes were differentially expressed after exposure to UV-C irradiation.

**Figure 4.** Overview of the transcriptome sequencing in *P. cuspidatum* leaves. (**A**) Gene expression profiles in the nine libraries. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively and (**B**) 3D principal component analysis showing clustering pattern among leaf samples based on global gene expression profiles. **Figure 4.** Overview of the transcriptome sequencing in *P. cuspidatum* leaves. (**A**) Gene expression profiles in the nine libraries. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively and (**B**) 3D principal component analysis showing clustering pattern among leaf samples based on global gene expression profiles.

#### *2.3. Differential Expressed Genes (DEG) between Control and UV-C Treated Leaf Samples 2.3. Di*ff*erential Expressed Genes (DEG) between Control and UV-C Treated Leaf Samples*

contributed by these DEGs (Figure 5B).

family members (Table S3).

We compared the gene expression levels between control samples (PC) to UV-C treated samples PC6H and PC12H with the aim to detect the genes affected by the UV-C treatment. Concerning PC vs. PC6H, we identified 38,985 differentially expressed genes (DEGs) with 17,859 and 21,126 up- and down-regulated, respectively. A slightly lower number of DEGs (32,312) was found for PC vs. PC12H, including 14,416 and 17,896 up and down-regulated genes, respectively. Cross-comparison of the two sets of DEGs showed that in total 45,222 genes were affected by the UV-C treatment and 26,075 DEGs were constantly conserved after 6 h and 12 h post UV-C irradiation (Figure 5A). These genes represent the key genes involved in the metabolic changes in response to UV-C treatment. To get insight into the biological pathways contributed by these DEGs, we performed KEGG enrichment analysis. The results indicated that the DEGs play various roles but were mainly involved in We compared the gene expression levels between control samples (PC) to UV-C treated samples PC6H and PC12H with the aim to detect the genes affected by the UV-C treatment. Concerning PC vs. PC6H, we identified 38,985 differentially expressed genes (DEGs) with 17,859 and 21,126 up- and down-regulated, respectively. A slightly lower number of DEGs (32,312) was found for PC vs. PC12H, including 14,416 and 17,896 up and down-regulated genes, respectively. Cross-comparison of the two sets of DEGs showed that in total 45,222 genes were affected by the UV-C treatment and 26,075 DEGs were constantly conserved after 6 h and 12 h post UV-C irradiation (Figure 5A). These genes represent the key genes involved in the metabolic changes in response to UV-C treatment. To get insight into the biological pathways contributed by these DEGs, we performed KEGG enrichment analysis. The results indicated that the DEGs play various roles but were mainly involved in metabolic pathways

metabolic pathways and biosynthesis of secondary metabolites. In addition, ribosome, plant hormone transduction and phenylpropanoid biosynthesis were the other enriched pathways

C exposure in the *Polygonum cuspidatum* leaf. Surprisingly, nearly half (2461) of the total annotated TF genes (5526) in *P. cuspidatum* was found differentially expressed and included 2132 and 1895 TFs DEGs detected at 6 h and 12 h, respectively. The majority of these genes were MYB, bHLH and ERF

Gene expression levels are modulated by regulators such as transcription factors [29]. We

and biosynthesis of secondary metabolites. In addition, ribosome, plant hormone transduction and phenylpropanoid biosynthesis were the other enriched pathways contributed by these DEGs (Figure 5B). *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 6 of 16

**Figure 5.** Differentially expressed genes (DEG) between the UV-C treated samples and the untreated control. (**A**) Venn diagram showing the unique and conserved DEGs between PC vs. PC6H and PC vs. PC12H; (**B**) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of all the identified DEGs. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively. **Figure 5.** Differentially expressed genes (DEG) between the UV-C treated samples and the untreated control. (**A**) Venn diagram showing the unique and conserved DEGs between PC vs. PC6H and PC vs. PC12H; (**B**) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of all the identified DEGs. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively.

*2.4. Mapping of DEGs Related to the Resveratrol Biosynthetic Pathway*  Deamination of phenylalanine ammonia by PAL is the first step in the resveratrol biosynthesis pathway. Then, the conversions of cinnamic acid into p-coumaric acid and subsequently into 4 coumarate-CoA are catalyzed by C4H and 4CL, respectively. The last step in the pathway consists of the conversion of one 4-coumarate-CoA and three malonyl-CoA units into resveratrol or naringenin chalcone by STS or CHS, respectively. Later, resveratrol is converted into pterostilbene by resveratrol Gene expression levels are modulated by regulators such as transcription factors [29]. We extended the analysis to detect the major transcription factor families involved in the response to UV-C exposure in the *Polygonum cuspidatum* leaf. Surprisingly, nearly half (2461) of the total annotated TF genes (5526) in *P. cuspidatum* was found differentially expressed and included 2132 and 1895 TFs DEGs detected at 6 h and 12 h, respectively. The majority of these genes were MYB, bHLH and ERF family members (Table S3).

O-methyltransferase (ROMT). We searched within the DEGs all genes encoding the key enzymes

#### involved in the resveratrol biosynthesis pathway and successfully identified 37 related DEGs, *2.4. Mapping of DEGs Related to the Resveratrol Biosynthetic Pathway*

of the leaves in the dark post UV-C irradiation as compared to 12 h.

including 36 and 34 at 6 h and 12 h post UV-C exposure, respectively. These genes included 10 PAL, eight *C4H*, nine *4CL*, two *STS*, six *CHS* and two *ROMT*. We mapped the 37 DEGs related to the resveratrol biosynthesis along with their gene expression fold change in order to understand the effect of UV-C radiation on resveratrol biosynthesis (Figure 6). Interestingly, we observed that all the 29 genes directly involved in the resveratrol biosynthesis (*PAL, C4H, 4CL* and *STS*) were significantly up-regulated after exposure to UV-C. Meanwhile, the six *CHS* genes were found strongly down-regulated. Concerning the DEGs encoding ROMT, all of the two genes were found sharply induced by UV-C treatment in *Polygonum cuspidatum* leaf (Figure 6). Globally, it appeared that the regulation of these 37 key genes was more intense after 6 h incubation Deamination of phenylalanine ammonia by PAL is the first step in the resveratrol biosynthesis pathway. Then, the conversions of cinnamic acid into p-coumaric acid and subsequently into 4-coumarate-CoA are catalyzed by C4H and 4CL, respectively. The last step in the pathway consists of the conversion of one 4-coumarate-CoA and three malonyl-CoA units into resveratrol or naringenin chalcone by STS or CHS, respectively. Later, resveratrol is converted into pterostilbene by resveratrol O-methyltransferase (ROMT). We searched within the DEGs all genes encoding the key enzymes involved in the resveratrol biosynthesis pathway and successfully identified 37 related DEGs, including 36 and 34 at 6 h and 12 h post UV-C exposure, respectively. These genes included 10 PAL, eight *C4H*, nine *4CL*, two *STS*, six *CHS* and two *ROMT*.

We mapped the 37 DEGs related to the resveratrol biosynthesis along with their gene expression fold change in order to understand the effect of UV-C radiation on resveratrol biosynthesis (Figure 6). Interestingly, we observed that all the 29 genes directly involved in the resveratrol biosynthesis (*PAL, C4H, 4CL* and *STS*) were significantly up-regulated after exposure to UV-C. Meanwhile, the six *CHS* genes were found strongly down-regulated. Concerning the DEGs encoding ROMT, all of the two genes were found sharply induced by UV-C treatment in *Polygonum cuspidatum* leaf (Figure 6). Globally, it appeared that the regulation of these 37 key genes was more intense after 6 h incubation of the leaves in the dark post UV-C irradiation as compared to 12 h. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 7 of 16

**Figure 6.** Proposed model of the molecular mechanism leading to the high accumulation of resveratrol after treatment with UV-C in *Polygonum cuspidatum* leaf. Deamination of phenylalanine ammonia by phenylalanine ammonia-lyase (PAL) is the first step in the resveratrol biosynthesis pathway. Then, the conversions of cinnamic acid into *p*-coumaric acid and subsequently into 4-coumarate-CoA are catalyzed by cinnamic acid 4-hydroxylase (C4H) and 4-coumarate-CoA ligase (4CL), respectively. The last step in the pathway consists of the conversion of one 4-coumarate-CoA and three malonyl-CoA units into resveratrol or naringenin chalcone by stilbene synthase (STS) or chalcone synthase (CHS), respectively. Later, resveratrol is converted into pterostilbene by resveratrol O-methyltransferase (ROMT). PAL, C4H, 4CL and STS were found strongly up-regulated in PC6H and PC12H as compared to PC, while CHS displayed the opposite trend. STS and CHS share the same substrate. *P. cuspidatum* tends to prioritize resveratrol accumulation by diverting the substrate **"**one 4-coumarate-CoA and three malonyl-CoA units" to the resveratrol synthesis pathway over the naringenin chalcone synthesis pathway through up-regulation of STS genes and down-regulation of CHS genes in response to UV-C exposure. High induction of ROMT also suggests that pterostilbenes may also be accumulated. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively. **Figure 6.** Proposed model of the molecular mechanism leading to the high accumulation of resveratrol after treatment with UV-C in *Polygonum cuspidatum* leaf. Deamination of phenylalanine ammonia by phenylalanine ammonia-lyase (PAL) is the first step in the resveratrol biosynthesis pathway. Then, the conversions of cinnamic acid into *p*-coumaric acid and subsequently into 4-coumarate-CoA are catalyzed by cinnamic acid 4-hydroxylase (C4H) and 4-coumarate-CoA ligase (4CL), respectively. The last step in the pathway consists of the conversion of one 4-coumarate-CoA and three malonyl-CoA units into resveratrol or naringenin chalcone by stilbene synthase (STS) or chalcone synthase (CHS), respectively. Later, resveratrol is converted into pterostilbene by resveratrol O-methyltransferase (ROMT). PAL, C4H, 4CL and STS were found strongly up-regulated in PC6H and PC12H as compared to PC, while CHS displayed the opposite trend. STS and CHS share the same substrate. *P. cuspidatum* tends to prioritize resveratrol accumulation by diverting the substrate "one 4-coumarate-CoA and three malonyl-CoA units" to the resveratrol synthesis pathway over the naringenin chalcone synthesis pathway through up-regulation of STS genes and down-regulation of CHS genes in response to UV-C exposure. High induction of ROMT also suggests that pterostilbenes may also be accumulated. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively.

In order to validate the results obtained through the RNA-seq, we selected 12 candidate DEGs from the enzymes involved in the resveratrol biosynthesis and some highly regulated transcription factors genes and conducted a qRT-PCR analysis (Table S4). As shown in the Figure 7, the expression levels of all the tested genes were significantly altered after treatment with UV-C and the qRT-PCR

*2.5. qRT-PCR Validation of Candidate Genes* 

#### *2.5. qRT-PCR Validation of Candidate Genes*

In order to validate the results obtained through the RNA-seq, we selected 12 candidate DEGs from the enzymes involved in the resveratrol biosynthesis and some highly regulated transcription factors genes and conducted a qRT-PCR analysis (Table S4). As shown in the Figure 7, the expression levels of all the tested genes were significantly altered after treatment with UV-C and the qRT-PCR results were strongly correlated with the RNA-seq reports (R<sup>2</sup> = 0.81, Figure S2). These results demonstrate the reliability of the RNA-seq data obtained in the present study. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 8 of 16

**Figure 7.** Quantitative real time PCR validation of selected candidate genes. The error bar represents the SD of three biological replicates. The *Actin* gene was used as the internal reference gene for normalization; PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C **Figure 7.** Quantitative real time PCR validation of selected candidate genes. The error bar represents the SD of three biological replicates. The *Actin* gene was used as the internal reference gene for normalization; PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively.

#### treatment, respectively. **3. Discussion**

**3. Discussion**  Resveratrol exhibits diverse beneficial properties in humans, including anti-inflammatory effects, anti-tumor activities and anti-aging effects, which have drawn extensive studies on this precious molecule [12]. As a phytoalexin, resveratrol plays important functions as an antimicrobial and antioxidant compound in plant responses to environmental stresses, such as UV irradiation and fungal infection [15]. Plants have developed various alleviation mechanisms to mitigate short wavelength UV-C radiation damaging effects, including a strong accumulation of UV-absorbing phenolic and flavonoid molecules in epidermal cells to block light penetration and anti-oxidative molecules to limit photo-oxidative damage [30–32]. In grape, UV-C irradiation induces a strong accumulation of resveratrol [28,33,34]. It was reported that UV-C could significantly stimulate the synthesis of resveratrol in *Vitis vinifera* × *V. amurensis*, mainly in the form of trans-resveratrol [35]. Similarly, UV-C treatments increased resveratrol synthesis in *Gnetum parvifolium* [36]. These reports suggest that the UV-C induction of resveratrol is quite conserved in stilbenoid-synthesizing plants. Although, *P. cuspidatum* represents the highest natural source of resveratrol known to date, no information is available regarding the effect of abiotic stress such as UV-C treatment on the level of resveratrol. Herein, we studied the effect of UV-C irradiation on the resveratrol metabolism in *P.*  Resveratrol exhibits diverse beneficial properties in humans, including anti-inflammatory effects, anti-tumor activities and anti-aging effects, which have drawn extensive studies on this precious molecule [12]. As a phytoalexin, resveratrol plays important functions as an antimicrobial and antioxidant compound in plant responses to environmental stresses, such as UV irradiation and fungal infection [15]. Plants have developed various alleviation mechanisms to mitigate short wavelength UV-C radiation damaging effects, including a strong accumulation of UV-absorbing phenolic and flavonoid molecules in epidermal cells to block light penetration and anti-oxidative molecules to limit photo-oxidative damage [30–32]. In grape, UV-C irradiation induces a strong accumulation of resveratrol [28,33,34]. It was reported that UV-C could significantly stimulate the synthesis of resveratrol in *Vitis vinifera* × *V. amurensis*, mainly in the form of trans-resveratrol [35]. Similarly, UV-C treatments increased resveratrol synthesis in *Gnetum parvifolium* [36]. These reports suggest that the UV-C induction of resveratrol is quite conserved in stilbenoid-synthesizing plants. Although, *P. cuspidatum* represents the highest natural source of resveratrol known to date, no information is available regarding the effect of abiotic stress such as UV-C treatment on the level of resveratrol. Herein, we studied the effect of UV-C irradiation on the resveratrol metabolism in *P. cuspidatum* mature leaves. As expected, we found that *P. cuspidatum* leaf contained a very high level of resveratrol (1000 µg/g

In grape, irradiation strength ranging from 30 to 510 W up to 1 min, followed by incubation in 2–4 days, resulted in 10.8-fold higher accumulation of resveratrol than that observed in the untreated control [37,38]. Cantos et al. [39] also showed that a distance of 40 cm, irradiation time of 30 s, source power of 500 W and storage time of 3 days generated the highest resveratrol accumulation. In this

*cuspidatum* mature leaves. As expected, we found that *P. cuspidatum* leaf contained a very high level of resveratrol (1000 µg/g FW) and UV-C could significantly induce the resveratrol accumulation FW) and UV-C could significantly induce the resveratrol accumulation (Figure 1), indicating that UV-C irradiation represented a good prospect for increasing resveratrol content in *P. cuspidatum*. Different UV irradiation treatments resulted in various levels of resveratrol, therefore an optimized UV treatment protocol is important for a significant induction of resveratrol. In grape, irradiation strength ranging from 30 to 510 W up to 1 min, followed by incubation in 2–4 days, resulted in 10.8-fold higher accumulation of resveratrol than that observed in the untreated control [37,38]. Cantos et al. [39] also showed that a distance of 40 cm, irradiation time of 30 s, source power of 500 W and storage time of 3 days generated the highest resveratrol accumulation. In this study, we incubated the UV-C irradiated leaves for 6 and 12 h. We observed an increase of 2.6-fold and 1.6-fold of resveratrol after 6 h and 12 h incubation times, respectively, suggesting that the effect of the UV-C treatment diminished after 6 h. Although the induction levels were not high as compared to studies in grape [37,38], based on our results and pending more optimization of the UV-C treatment protocol, we might recommend 10 min UV-C irradiation and 6 h incubation in the dark for obtaining a relatively high resveratrol in *P. cuspidatum* leaves.

To understand the molecular mechanism underlying the strong accumulation of resveratrol in *P. cuspidatum* after exposure to UV-C, we de novo sequenced the transcriptome and assembled the unigenes (Tables 1 and 2; Figure 2). In total, 165,013 unigenes were identified in this study, a number that is quite the double of the unigenes number (86,418) detected by Hao et al. [14]. The high discrepancy between the detected unigene numbers in both studies could result from differences in the tissue types, the employed software for assembly and more importantly the sequencing platform. Differential gene expression (DEG) analysis showed that UV-C treatments affected nearly <sup>1</sup> 4 of the expressed genes and various cellular processes including hormones and secondary metabolism were altered (Figure 5). Our results were in perfect concordance with works by Yin et al. [40], who demonstrated that more than 100 functional subcategories were contributed by the DEGs between UV-treated grape berries and untreated samples and particularly, "response to stress" and "metabolic processes" were the most represented terms.

Previous studies have demonstrated that UV-C irradiation alters the activity levels of several structural genes participating in the resveratrol biosynthesis pathway [35,36,40–42]. Within the DEGs detected in this study, 37 were mapped to the resveratrol biosynthesis pathway (Figure 6). Interestingly, all the genes directly involved in the resveratrol biosynthesis (10 phenylalanine ammonia-lyase (PAL), 8 cinnamic acid 4-hydroxylase (C4H), 9 4-coumarate-CoA ligase (4CL) and 2 stilbene synthase (STS)) were significantly up-regulated in UV-C treated samples as compared to untreated samples (Figure 6). Conversely, the six chalcone synthase (CHS) DEGs were all down-regulated in samples exposed to UV-C (Figure 6). In fact, CHS and STS enzymes use the same substrate (one 4-coumarate-CoA and three malonyl-CoA units) for the production of naringenin chalcone and resveratrol, respectively. This result indicates that after exposure to UV-C, *P. cuspidatum* diverts the substrate "one 4-coumarate-CoA and three malonyl-CoA units" to the resveratrol synthesis pathway over the naringenin chalcone synthesis pathway by up-regulating STS and down-regulating CHS. The conclusions of several authors, including Xi et al. [41], Suzuki et al. [42] and Yin et al. [40] support well our findings as they demonstrated that since STS and CHS share the same substrate, there may be a competitive and/or inhibitory relationship between them in response to UV-C exposure, which may in turn play a vital role in resveratrol accumulation in grape berries.

The tyrosine ammonia-lyase (TAL) enzyme can also utilize L-tyrosine as a substrate to produce p-coumaric acid, which subsequently is converted into resveratrol by STS [43]. In this study, we did not find any *TAL* gene among the DEGs, indicating that the high accumulation of resveratrol after UV-C treatments was essentially due to the strong conversion of phenylalanine and cinnamic acid through PAL and C4H, respectively (Figure 6).

Trans-resveratrol is highly instable upon exposure to light and oxygen or under harsh PH, leading to the reduction its bioavailability and bioactivity [44]. Other natural stilbenes derived from resveratrol such as pterostilbene and pinostilbene, display higher oral bioavailability and bioactivity, therefore, efforts are ongoing to develop strategies to obtain these bioactive resveratrol derivatives in abundance. An efficient technique to achieve this goal is the manipulation of resveratrol O-methyltransferase (ROMT), which is the enzyme that converts resveratrol into pterostilbenes [27,45,46]. Here, we observed that 2 ROMT genes were significantly up-regulated in UV-C treated *P. cuspidatum* leaf samples (Figure 6), suggesting that UV-C treatments not only increase the trans-resveratrol levels but may also increase the methylated derivatives of resveratrol. Comparing 6 h and 12 h incubation times with respect to the expression levels of resveratrol biosynthesis related DEGs, we found that gene induction/repression was more vigorous at 6 h than 12 h, which correlates with the higher resveratrol content observed at this time point through resveratrol quantification.

Transcription factors (TF) such as MYB and bHLH were reported to regulate the expression levels of several genes involved in the phenylpropanoid pathways [47,48]. Following UV-C treatment, the gene Myb14 has been shown to activate *STS* genes in grape [20,35,40,49]. In this study, we found 2461 TFs among the DEGs with members of MYB, bHLH and ERF being the most active after UV-C irradiation (Table S3) and might play crucial roles in UV-C induced responses in *P. cuspidatum*. We did not observe a clear trend in the differential expression of genes belonging to these families as many family members were either down-regulated or up-regulated after exposure to UV-C. This highlights the complex mechanisms of UV-C transcriptional regulation and renders the identification of potential master regulators of resveratrol biosynthetic genes very difficult. Gene co-expression analysis is a widely used technique to break down large transcriptome datasets into co-expressed modules in order to pinpoint key TFs modulating important structural genes [50–56]. For example, this approach has been used to discover WRKY TFs (WRKY24/28/29/37/41) that were co-expressed simultaneously with eight *STS* genes (*STS12*/*13*/*16*/*17*/*18*/*24*/*27*/*29*) in roots and leaves of *Vitis vinifera* [57]. Later on, TFs belonging to different families such as MYB, WRKY and ERF were shown to putatively contribute to STS regulation [58,59]. Therefore, we proposed further RNA-sequencing in *P. cuspidatum* using various genotypes and tissues to comprehensively undertake an integrated gene co-expression analysis in order to find out the key regulators of the resveratrol biosynthetic genes in this important medicinal plant species.

#### **4. Materials and Methods**

#### *4.1. Plant Materials*

*Polygonum cuspidatum* Sieb. et Zucc. was used as plant material in this study. Plants were maintained in the medicinal plant garden at Yangtze University, Jingzhou, China. Healthy, mature (30-days old) leaves of similar size were detached from the shoot and were immediately immersed in water and subsequently transferred into ddH2O [60]. All leaves were incubated in the dark at 25 ◦C for 30 min. Then, the leaf abaxial surfaces were exposed for 10 min to 6 W/m<sup>2</sup> UV-C irradiation provided by a UV-C lamp (Model CJ-1450, Sujie Purification, China). Samples were collected at 0 h (before UV-C treatment, PC), 6 h (PC6H) and 12 h (PC12H) after the initiation of treatment. Each treatment was performed three times and each replication consisted of six leaves. After sampling, the leaves were ground into powder in liquid nitrogen and stored at −80 ◦C until analysis.

#### *4.2. Quantification of Resveratrol in Polygonum cuspidatum Leaves*

The sample preparation, extract analysis, resveratrol identification and quantification were performed at Wuhan MetWare Biotechnology Co., Ltd, Wuhan, China. (www.metware.cn) following their standard procedures and previously described by Zhang et al. [61]. The experimental measurements were carried out in triplicate and results were presented as average of three analyses ± standard deviation. Statistical analysis of the data was done using the (www.r-project.org) using the one-way analysis of variance for testing statistical significance between different samples. Mean comparisons were done using the Tukey HSD test.

#### *4.3. RNA Extraction, cDNA Library Construction and Transcriptome Sequencing*

The leaf samples from PC, PC6H and PC12H were used for total RNAs extraction employing the Spin Column Plant total RNA Purification Kit according to the manufacturer's protocol (Sangon Biotech, Shanghai, China). Purity of the extracted RNAs was assessed on 1% agarose gel followed by a NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA quantification was performed using the Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). Next, RNA integrity was checked by the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Sequencing libraries was created using NEB Next Ultra RNA Library Prep Kit following the manufacturer's instructions. The index codes were added to each sample. Briefly, the mRNA was purified from 3 µg total RNA of each of three replicate using poly-T oligo-attached magnetic beads and then broken into short fragments to synthesize first strand cDNA. The second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. PCR was carried out with Phusion High Fidelity DNA polymerase using universal PCR primers and index (X) primer. Finally, six paired-end cDNA libraries with an insert size of 150 bp were constructed for transcriptome sequencing and sequenced on HiSeqTM 2000 platform (Illumina Inc., San Diego, CA, USA).

#### *4.4. Quality Check, Cleaning and de novo Assembly*

The clean reads were retrieved after trimming adapter sequences, removal of low quality (containing >50% bases with a Phred quality score <10) and reads with unknown nucleotides (more than 5% ambiguous residues N) using the FastQC tool (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/). The high quality reads from all the nine libraries were de novo assembled into transcripts using the Trinity software (version r20140717, [62]) with the following parameters: -min-contig-length, 100-min-glue and 3–path-sing-distance-85-min-kmer-cov3. Next, the transcripts were realigned to construct unigenes and the software TGICL, version v2.1 [63] was employed to eliminate the redundant unigenes and to get the final unigene list based on the following parameters: -l, 40, -c, 10, -v, 25, -O, '-repeat\_stringency, 0.95, -minmatch, 35 and -minscore 35'.

#### *4.5. Functional Annotation of the Unigenes*

The assembled unigenes were annotated by searching against various databases such as Kyoto Encyclopedia of Genes and Genomes (KEGG) [64], Gene Ontology (GO) [65], Clusters of Orthologous Groups (COG) [66], Swissprot [67], NR [68], euKaryotic Orthologous Groups (KOG) [69] and NT using BLAST version:v2.2.26 [70] with a threshold of *<sup>e</sup>*-value <sup>&</sup>lt; <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> . Based on the functional annotation results, according to the database priorities of NR, SwissProt, KEGG and COG, we selected the best comparison fragment of unigenes as the CDS for the unigenes. For unigenes, which failed to return a hit, we then used ESTScan v3.0.2 [71] to make predictions using default parameters.

#### *4.6. Analysis of Transcription Factors*

To identify the gene encoding transcription factors (TF), we first collected the HMM profiles of known transcription factors from various databases (PlantTFDB [72], AnimalTFDB [73], FTFD [74] and DBD [75]) and the literature, and then used the hmmsearch of the HMMER package version 3.1b2 to compare the protein sequence to the HMM of known TFs, with the *<sup>e</sup>*-value <sup>&</sup>lt; <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> .

#### *4.7. Gene Expression and Di*ff*erential Expression Analysis*

The clean reads were compared to the assembled unigenes using the tool Bowtie2 version 2.1.0 [76] and then the number of reads on each unigene was calculated using the RSEM version 1.2.21 [77]. The gene expression level was determined according to the fragments per kilobase of exon per million fragments mapped (FPKM).

The read count was normalized and DESeq2 tool version 1.4.5 [78] was used to detect the differential expressed genes (DEGs) between the control samples (PC) and the UV-C treated samples (PC6H and PC12H) with the fold change of >2 [79] and false discovery rate (FDR) correction set at *p* < 0.01. GO enrichment analysis was performed using the GOseq tool version 1.16.2 [80] with *p* < 0.05. KEGG pathway enrichment analysis of the DEGs was done using KOBAS2.0 [81] with FDR correction (*p* < 0.05).

#### *4.8. Gene Expression Using Quantitative Real Time-PCR*

The qRT-PCR was performed on RNA extracted from control and stressed leaf samples as described by Dossa et al. [82] using the *Actin* gene as the internal control. Specific primer pairs of 12 selected genes were designed using the Primer Premier 5.0 [83] (Table S4). The qRT-PCR was conducted on a Roche Lightcyler®480 instrument using the SYBR Green Master Mix (Vazyme), according to the manufacturer's protocol. Each reaction was performed using a 20 µL mixture containing 10 µL of 2× ChamQ SYBR qPCR Master Mix, 6 µL of nuclease-free water, 1 µL of each primer (10 mM) and 2 µL of four-fold diluted cDNA. All of the reactions were run in 96-well plates and each cDNA was analyzed in triplicate. The following cycling profile was used: 95 ◦C for 30 s, followed by 40 cycles of 95 ◦C/10 s, 60 ◦C/30 s. Data are presented as relative transcript level based on the 2−∆∆Ct method [84].

#### **5. Conclusions**

In summary, we showed that UV-C treatments induced the accumulation of trans-resveratrol in *P. cuspidatum* leaves. Various UV-C treatments generated different levels of accumulation suggesting that there was room for optimization of the UV-C treatment protocol in order to yield maximum content of resveratrol. We further provided the putative genes participating in the resveratrol biosynthetic pathway and highlighted the key genes differentially expressed upon exposure to UV-C. It was evident that *P. cuspidatum* prioritized the resveratrol synthesis by strongly up-regulating the genes directly involved in this pathway and shutting down the expression of chalcone synthase genes. The results from this study provided insights into the mechanism of UV-C induced accumulation of resveratrol and probably its methylated forms in a species other than grape. It lays the foundation for further enhancement of stilbenes in *P. cuspidatum*.

**Supplementary Materials:** Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/24/ 6185/s1. Figure S1 Gene ontology classification of all unigenes detected in this study; Figure S2 Correlation of the qRT-PCR expression profiles and the RNA-seq based on the selected genes. PC, PC6H and PC12H represent samples collected before, 6 h and 12 h after UV-C treatment, respectively; Table S1 List of the genes encoding enzymes involved in the resveratrol biosynthetic pathway; Table S2 Statistics of the read mapping results; Table S3 List of the differentially expressed transcription factors detected in this study; Table S4 The primer sequences of genes used for real time quantitative PCR.

**Author Contributions:** Conceptualization, H.Z., J.Q.; methodology, Z.L., J.X., X.W., Y.W.; validation, Y.L., D.W.; formal analysis, Z.L.; investigation, Z.L., J.X., X.W., Y.W.; resources, Y.L., D.W.; writing—original draft preparation, Z.L.; writing—review and editing, J.Q.; supervision, J.Q.; project administration, H.Z., J.Q.; funding acquisition, H.Z., J.Q.

**Funding:** This work was financially supported by the National Natural Science Foundation of China (Grant No. 81803670).

**Conflicts of Interest:** The authors declare no conflict of interest

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
