A Transcriptomic Analysis of Gene Expression in Chieh-Qua in Response to Fusaric Acid Stress

Fusarium wilt results in undesirable effects on the quality and production of chieh-qua (Benincasa hispida Cogn. var. Chieh-qua How). Fusaric acid (FA), a secondary metabolite of biotin produced by pathogens of genus Fusarium, induced resistant responses in chieh-qua; however, the physiological and molecular mechanism(s) of FA resistance remains largely unknown. In our study, ‘A39’ (FA-resistant cultivar) exhibited decreased malondialdehyde (MDA) content and increased superoxide dismutase (SOD) enzyme activity when exposed to FA compared with ‘H5’ (FA-susceptible cultivar). More apoptosis cells existed in ‘H5’ than ‘A39’ after 2 days of FA treatment. RNA-seq results revealed that a total of 2968 and 3931 differentially expressed genes (DEGs) were detected under normal conditions (1562 up-regulated and 1406 down-regulated) and FA treatment (2243 up-regulated and 1688 down-regulated), respectively. Interestingly, DEGs associated with pathogen-related protein and ethylene (ET) biosynthesis and signal pathways were most significantly changed during FA stress. Notably, several crucial genes encoding pathogenesis-related protein (CL4451.Contig2, CL2175.Contig4), peroxidase (Unigene49615 and CL11695.Contig2), and ET-responsive transcription factors (TFs) (CL9320.Contig1, CL9849.Contig3, CL6826.Contig2, CL919. Contig6, and CL518.Contig7) were specifically induced after FA treatment. Collectively, the study provides molecular data for isolating candidate genes involved in FA resistance, especially ET related genes in chieh-qua.


Introduction
Chieh-qua (Benincasa hispida Cogn. var. Chieh-qua How), a subspecies of wax gourd, is widely distributed in South China and Southeast Asian countries [1,2]. During growth, chieh-qua faces constantly changing environmental conditions including biotic stresses (herbivore attack, pathogen infection) and abiotic stresses (drought, high or low temperatures) [3]. Among unfavorable biotic stresses, fusarium wilt (FW), is a worldwide, disruptive soil borne disease, and is also the most threatening pathogen, causing severe decreases in production and quality in chieh-qua [2,4]. The FW pathogen first enters into the host organization by root wounds or micropores at the tip of root hairs, infecting the xylem vessels leading to a reddish-brown discoloration of the rhizome [5]. Fusaric acid (5-butylpyridine-2-carboxylic acid, FA), a non-host-specific phytotoxin of Fusarium species, is reported to be involved in pathogenecity [6,7]. After fusarium attack, much higher concentration was detected in plant tissues [8,9]. FA exerts essential roles in increasing FW symptoms in plants [10][11][12] and the content of FA is related to incidence rate in banana [9] and chieh-qua [13]. Therefore, FA has been considered as an identification index in screening for the FW response [1,14]. However, molecular studies of FA effects on chieh-qua have not been reported. 'A39' and 'H5' are both homozygous inbred lines with different FA resistance, identified in our previous study. We also found that 'A39' has a tolerance to high temperature and water-deficiency, while 'H5' showed sensitivity to both stresses [28,29]. To start, seeds were soaked in warm water for 6-8 h and germinated for 2-3 days in an incubator at 30 • C. Then, the germinated seeds were grown in feeding blocks containing soil under 14 h/10 h (5500 lux) for 30 • C/24 • C in day/night, respectively. After one month, entire seedlings (2 true leaves) was taken out of feeding blocks, rhizosphere soil was cut partially off, and roots were washed clearly. Subsequently, fibrous roots were partially dried, maintaining about 5 cm length of roots. Then, seedlings were transferred to water (normal condition) or FA solution (200 mg/L) without nutrient substances. Two days later, a total of 10 normal leaves were randomly sampled under normal and FA stress, respectively. Each biological replicate included 10 leaves from 10 randomly selected plants. The selected leaves were immediately frozen in liquid nitrogen and stored at −80 • C for further study.

Measurement of Malondialdehyde (MDA) Content and Superoxide Dismuatse (SOD) Activity
A total of 0.5 g fresh leaves from each replicate plant were put into a tube including 0.5% (v/v) trichloroacetic acid (TCA), ground into a powder under ice using a grinding rod, and centrifuged 20 min at 3000× g. After that, about 2 mL supernatant was pipetted out and the same volume of 0.5% TCA was added. Then, the mixture was boiled at 100 • C for 30 min. Finally, the absorption of the supernatants was detected at 450, 532, and 600 nm.
The SOD activity was detected as described by Marklund et al. (1974) [30]. Briefly, the leaves (1 g) were homogenized in liquid nitrogen using a mortar and pestle in 10 mL of 0.2 M lysate (citrate phosphate buffer (pH 6.5) with 0.5% Triton X-100). The mixture was centrifuged at 10,000× g for 30 min at 4 • C. Then, 0.3 mL riboflavin (20 µmol/L) was added to the supernatant with enzyme OD values were measured at 560 nm. The SOD activity required to inhibit the autooxidation of pyrogallol to 50% was defined as one enzyme unit (1 U).

Terminal Deoxynucleotidyl Transferase-Mediated UTP Nick-End Labeling (TUNEL) Assay
For each replicate sample, about 1 g of fresh leaves were collected and fixed in FAA fixation solution for one day. The solution included 18:1:1 (v/v) mixture of formalin, 70% ethanol, and acetic acid. After that, these selected leaves were dehydrated through 90%, 70%, and 50% ethanol in order, and inserted into paraffin (Sigma-Aldrich). After that, tissue cross-sections were cut with a rotary microtome to about 10 µm in thickness, hydrated through an ethanol series (100%, 85%, 70%, and 50%), and treated with proteinase K in phosphate buffer (pH 7.4-8.0). The Green Fluorescence indicates the signal for apoptosis. The Merged represents the combination of PI staining and Green Fluorescence. Signals were observed under a laser scanning confocal microscope at 488 nm (excitation), 520 nm (detection), 488 nm (excitation), and 610 nm (detection), respectively (LSM 700; Carl Zeiss, Germany).

Transcriptome Sequencing
RNA was extracted from three replicate frozen leaf samples of each treatment group using TRIZOL reagent according to manufacturer's instructions (Waltham, MS, USA). Then, genomic DNA was treated with RNase-free DNase which was removed by phenolchloroform re-extraction. The quality of RNA was detected using an Agilent 2100 bioanalyzer (Agilent, China), mRNA was purified by poly T oligos attached to magnetic beads and cut into 200-300 bp fragments, then the fragments were used as templates to synthesize first-stand cDNA by using 6-base random primers, A paired-end library was constructed from the cDNA synthesized using a Genomic Sample Prep Kit (Illumina) followed by ligation of adaptors and purification with AMPureXP beads, then the library fragments were constructed with PCR amplification and produced a 300-400 bp library.
The DNA library quality was tested using the Agilent 2100 Bioanalyzer, and the samples were sequenced by Illumina paired-end (PE) on the Illumina HiSeq 2000 platform (Illumina, San Diego, CA, USA) at Gene Denovo (Guangzhou, China). Three biological replicates were performed separately.

Screening and Significant Test for DEGs
By quantifying reads according to the FPKM (per kilo-base per million reads) method, gene expression levels were calculated. For determining gene transcriptional activity, the cutoff value was over 1 [31]. Then, NOISeq was used to detect differently expressed genes (DEGs) of both normal and FA stress transcriptome libraries based on the following criteria: fold change ≥2 or ≤0.5, divergence probability ≥ 0.8 and p value < 0.05. Gene Ontology (GO) enrichment for these DEGs was performed using WEGO software [32] by performing Fisher's exact test with a robust false discovery rate (FDR) correction to obtain an adjusted p-value between certain test gene groups and whole genome annotation. A pathway enrichment analysis was carried out on the basis of KEGG database in order to understand DEG biological functions [33].

Quantitative Real-Time PCR Analysis
After RNA extraction, a total of 20 µL cDNA was synthesized from 1 g RNA using the Quanti Tect Reverse Transcription Kit (Qiagen, China). Then, qRT-PCR was performed using a 20 µL reaction volume with 0.5 µL of cDNA, 0.2 µM of primer mix, and the TB Green Premix Ex Taq Kit (TaKaRa, Japan). The internal reference gene of chieh-qua was CqActin (F: TCAACCCAAAGGCTAACAG; R: CTTGTCCATC AGGCAGTTC) [28,29]. The assay was performed with three technical and biological replicates of each treatment group, respectively. All qRT-PCR primers are listed in Table S1.

Statistical Analyses
A cluster analysis was performed using MEGA6 software (Mega Limited, Auckland, New Zealand). The significant differences were detected by IBM SPSS Statistics 20 (Armonk, NY, USA). The GraphPad Prism 5 (GraphPad Software Inc., San Diego, CA, USA) was used for chart preparation. Gene relative expressions were calculated using the 2 −∆∆ Ct method [34].

Phenotype of 'A39' and 'H5' upon FA Stress
'A39' is a cultivar with comprehensive resistance to heat, drought, and FA stress, while 'H5' is sensitive to the three stresses [28,29]. When plants grew to two-true leaves, 'A39' and 'H5' seedlings showed vigorous growth and development under normal condition ( Figure 1A). Then, the roots of their seedlings were exposed to FA stress for 2 days ( Figure 1B). Results showed that 'H5' leaves, especially the cotyledons, exhibited more severe wilting and chlorosis than 'A39' ( Figure 1B,C). Although there was no difference in MDA content or SOD activity between 'H5' and A59 at normal conditions, significant increases of MDA ( Figure 1D) and less increase in SOD activity occurred in 'H5' in comparison to 'A39' after 2 days of FA stress ( Figure 1E), indicating that 'A39' had a better response to FA stress than 'H5'. formed using a 20 μL reaction volume with 0.5 μL of cDNA, 0.2 μM of primer mix, the TB Green Premix Ex Taq Kit (TaKaRa, Japan). The internal reference gen chieh-qua was CqActin (F: TCAACCCAAAGGCTAACAG; R: CTTGTCCATC GCAGTTC) [28,29]. The assay was performed with three technical and biological r cates of each treatment group, respectively. All qRT-PCR primers are listed in Table   2. 7

. Statistical Analyses
A cluster analysis was performed using MEGA6 software (Mega Limited, Auck New Zealand). The significant differences were detected by IBM SPSS Statistics 20 monk, NY, USA). The GraphPad Prism 5 (GraphPad Software Inc., San Diego, CA, U was used for chart preparation. Gene relative expressions were calculated using 2 −△ △ Ct method [34].

Phenotype of 'A39′ and 'H5′ upon FA Stress
'A39′ is a cultivar with comprehensive resistance to heat, drought, and FA st while 'H5′ is sensitive to the three stresses [28,29]. When plants grew to two-true lea 'A39′ and 'H5′ seedlings showed vigorous growth and development under normal dition ( Figure 1A). Then, the roots of their seedlings were exposed to FA stress for 2 ( Figure 1B). Results showed that 'H5′ leaves, especially the cotyledons, exhibited m severe wilting and chlorosis than 'A39′ ( Figure 1B, C). Although there was no differ in MDA content or SOD activity between 'H5′ and A59 at normal conditions, signif increases of MDA ( Figure 1D) and less increase in SOD activity occurred in 'H comparison to 'A39′ after 2 days of FA stress ( Figure 1E), indicating that 'A39′ had a ter response to FA stress than 'H5′. (D,E) MDA content and SOD activity before treatment and after 2 days of FA stress. Data is presented as the mean ± standard deviation (n = 9). ** Student's t-test at p < 0.01.

Cell Death under FA Stress
In order to examine possible cellular changes under FA stress, we performed a TUNEL assay to detect nuclear DNA fragmentation. Results showed that no positive TUNEL signals were detected in the leaf cells under normal conditions ( Figure 2A). However, more intense TUNEL signals were found in 'H5' than 'A39' after 1 day ( Figure 2B,C) and 2 days ( Figure 2D,E) of FA stress, suggesting that there was more cell degradation due to FA stress.
MDA content and SOD activity before treatment and after 2 days of FA stress. Data is presented as the mean ± standard deviation (n = 9). ** Student's t-test at p < 0.01.

Cell Death under FA Stress
In order to examine possible cellular changes under FA stress, we performed a TUNEL assay to detect nuclear DNA fragmentation. Results showed that no positive TUNEL signals were detected in the leaf cells under normal conditions ( Figure 2A). However, more intense TUNEL signals were found in 'H5′ than 'A39′ after 1 day ( Figure  2B,C) and 2 days ( Figure 2D,E) of FA stress, suggesting that there was more cell degradation due to FA stress.

Transcripts Assembly and SNP Analysis
From RNA-seq assays, we obtained a total range of 48~104 million clean reads for each sample after removing low-quality and adaptor-containing raw reads (Table S2). All clean reads of Q20 (%) were over 96% and the average numbers of detected genes was 72458 in 'A39' (normal), 71437 in 'H5' (normal), 81203 in 'A39' (FA), and 89939 in 'H5' (FA), respectively. In order to know the relationship among the 12 samples, we conducted PCA analysis with clean data of 12 samples. Results showed that samples clustered together from the two-and three-dimensional diagrams ( Figure 3). All the original data has been uploaded in the NCBI SRA database (SRP180414). After removing undifferentiated genes, we identified 2968 (Table S3) and 3931 (Table S4) differentially expressed genes (DEGs) under normal conditions and FA stress, respectively. Among them, 1562 genes were up-regulated and 1406 down-regulated (gene expression in 'A39' compared with that in 'H5') under normal conditions ( Figure 4A). Additionally, 2243 up-regulated and 1688 down-regulated genes were identified at 2 days after FA stress ( Figure 4B). pressed genes (DEGs) under normal conditions and FA stress, respectively. Among them, 1562 genes were up-regulated and 1406 down-regulated (gene expression in 'A39′ compared with that in 'H5′) under normal conditions ( Figure 4A). Additionally, 2243 up-regulated and 1688 down-regulated genes were identified at 2 days after FA stress ( Figure 4B).
In addition, to further verify the accuracy of transcriptome data found by RNA-Seq, a total of 17 DEGs (Table S5) were selected and performed using the qRT-PCR method. A strong positive correlation (R 2 = 0.938) between the RNA-seq data and qRT-PCR results was detected ( Figure S1), suggesting the accuracy of RNA-seq results. In addition, we applied the GATK (Genome Analysis Tool Kit) tool to analyze SNP variation. SNP detection analysis revealed 139, 819 SNPs in the twelve samples (Tables S6 and S7). Of these SNPs, the transition type occupied 72.01%, while the transversion type was 27.99% (Figure S2). Among transition SNPs, C/T and G/A were the larger part with 19.41% and 19.85%, respectively. The transversions included 8 types and each ratio of SNP type was similar ( Figure S2). These SNPs will be useful in improving the associated polymorphism for mapping potential candidate genes.

Functional Classification of FA Stress Response Genes
In order to determine the functional categories of the DEGs for the two cult Gene Ontology (GO) analysis was implemented. The Cellular Components (CC) cat occupied the largest number of DEGs and genes related to the integral compone In addition, to further verify the accuracy of transcriptome data found by RNA-Seq, a total of 17 DEGs (Table S5) were selected and performed using the qRT-PCR method. A strong positive correlation (R 2 = 0.938) between the RNA-seq data and qRT-PCR results was detected ( Figure S1), suggesting the accuracy of RNA-seq results. In addition, we applied the GATK (Genome Analysis Tool Kit) tool to analyze SNP variation. SNP detection analysis revealed 139, 819 SNPs in the twelve samples (Tables S6 and S7). Of these SNPs, the transition type occupied 72.01%, while the transversion type was 27.99% ( Figure S2). Among transition SNPs, C/T and G/A were the larger part with 19.41% and 19.85%, respectively. The transversions included 8 types and each ratio of SNP type was similar ( Figure S2). These SNPs will be useful in improving the associated polymorphism for mapping potential candidate genes.

Functional Classification of FA Stress Response Genes
In order to determine the functional categories of the DEGs for the two cultivars, Gene Ontology (GO) analysis was implemented. The Cellular Components (CC) category occupied the largest number of DEGs and genes related to the integral component of membrane were the maximum quantity under normal conditions ( Figure S3A). However, under FA stress, the largest number was also mainly in the category of CC (the integral component of membrane), followed by the nucleus ( Figure S3; Table S8).
To examine the DEG-associated pathways, the KEGG was analyzed and exhibited top 20 pathway enrichment of DEGs under normal conditions. The main pathways were starch and sucrose metabolism, splicesome, and plant pathogen interaction ( Figure 5A, Table S9). When exposed to FA stress, genes related to plant hormones signal transduction, splicesome, purine metabolism, and plant pathogen interaction were the most enriched pathways ( Figure 5B, Table S10), indicating these progresses might be involved in FA resistance.

Functional Classification of FA Stress Response Genes
In order to determine the functional categories of the DEGs for the two culti Gene Ontology (GO) analysis was implemented. The Cellular Components (CC) cate occupied the largest number of DEGs and genes related to the integral compone membrane were the maximum quantity under normal conditions ( Figure S3A). Ho er, under FA stress, the largest number was also mainly in the category of CC (the gral component of membrane), followed by the nucleus (Figure S3; Table S8).
To examine the DEG-associated pathways, the KEGG was analyzed and exhi top 20 pathway enrichment of DEGs under normal conditions. The main pathways starch and sucrose metabolism, splicesome, and plant pathogen interaction (Figur Table S9). When exposed to FA stress, genes related to plant hormones signal trans tion, splicesome, purine metabolism, and plant pathogen interaction were the mos riched pathways ( Figure 5B, Table S10), indicating these progresses might be involv FA resistance.

Expression of Genes Encoding Pathogen-Related Proteins
When plants are faced with attacks by potentially pathogenic organisms, they c acquire resistance by activating physiological defense mechanisms, c plant-pathogen interactions [35,36]. This induction course is always accompanied b activation of a set of genes and the accumulation of corresponding gene products, cially the pathogenesis-related (PR) proteins [37]. Interestingly, plant-pathogen int

Expression of Genes Encoding Pathogen-Related Proteins
When plants are faced with attacks by potentially pathogenic organisms, they could acquire resistance by activating physiological defense mechanisms, called plant-pathogen interactions [35,36]. This induction course is always accompanied by the activation of a set of genes and the accumulation of corresponding gene products, especially the pathogenesisrelated (PR) proteins [37]. Interestingly, plant-pathogen interaction GO terms enriched among the DEGs were differentially expressed between 'A39' and 'H5' under FA treatment. In particular, four genes encoding pathogen-related proteins (CL4451.Contig2, Unigene25112, CL1685.Contig27, and CL1130.Contig5) showed significantly higher expression levels in 'A39' seedlings than 'H5', respectively (Table 1; Figure 6). In response to various stresses, plants could accumulate reactive oxygen species (ROS) and some scavenger enzymes. For example, peroxidase played important roles in the control of ROS, which alleviated damage from stresses [38]. In our study, the relative expression of two peroxidase genes (Unigene49615 and CL11695.Contig2) were prominently induced in 'A39' after FA stress. These observations indicate that genes of PR proteins and peroxidase might contribute to plant resistance in 'A39'.  Figure 6). In response to various stresses, plants could accumulate reactive oxygen species (ROS) and some scavenger enzymes. For example, peroxidase played important roles in the control of ROS, which alleviated damage from stresses [38]. In our study, the relative expression of two peroxidase genes (Unigene49615 and CL11695.Contig2) were prominently induced in 'A39′ after FA stress. These observations indicate that genes of PR proteins and peroxidase might contribute to plant resistance in 'A39′.

Expression of Genes Related to Ethylene (ET)
Plant hormones play crucial roles when plants are faced with abiotic and biotic stresses [39]. DEGs related to plant hormone signal transduction were differentially expressed between 'A39' and 'H5' under FA treatment. Among these hormones, ethylene (ET) is an important regulator of the plant immune system [40]. In our study, six genes encoding ET-responsive transcription factors (CL9320.Contig1, CL6826.Contig2, CL919.Contig6, CL518.Contig7), ACC oxidase (CL787.Contig1) ET insensitive-like 1 protein isoform X1(CL9849. Contig3) were selected and their expression detected. Combining the RNA-Seq data, six genes were highly induced in 'A39' compared with 'H5' under FA stress, which was consistent with our qRT-PCR data (Figure 7, Table 2). These results indicated that ET-related genes were involved in the response to FA stress.

Discussion
In this study, 'A39′ showed increased tolerance to FA stress compared with 'H5 our previous studies, 'A39′ also exerted high temperature and drought resistanc contrast to 'H5′ [28,29]. These findings seem to indicate that the chieh-qua cultivar 'A exhibited simultaneous tolerance to multiple stresses, including heat and drought, wh are two negative environmental factors that drastically affect crop growth. Fusarium w as a serious biotic stress, significantly influences chieh-qua productivity and qua Here, our study mainly focused on exploring genes involved in FA resistance. Howe the relationship among different stresses would be necessary for studying the molec mechanisms of stresses on chieh-qua.
At present, the genome sequence information on wax gourd [2,41] and the RNA data of chieh-qua [28,29] provide molecular data for performing genetic research chieh-qua. Nevertheless, the publicly available data is not sufficient to explain the lecular mechanisms underlying resistance to FW/FA in chieh-qua. Therefore, much m extensive molecular data is needed to discover genes related to FW resistance chieh-qua. This could serve as a good source for deeply isolating useful genes involve disease resistance or plant-pathogen interaction. Using transcriptome analysis, mult DEGs were screened under normal conditions and after FA stress. Among them, sev DEGs related to pathogen-related proteins, peroxidase, and ET-responsive transcrip factors were possibly involved in FA resistance with significant differences in expres between 'A39′ and 'H5′. In view of 'A39′ having comprehensive resistance to diffe stresses, we speculated that common regulatory networks contributing to its greater sistance might exist. For instance, under drought stress [29] and FA treatment, we fo that there was a common pathway of plant pathogen interaction by KEGG analysis.

Discussion
In this study, 'A39' showed increased tolerance to FA stress compared with 'H5'. In our previous studies, 'A39' also exerted high temperature and drought resistance in contrast to 'H5' [28,29]. These findings seem to indicate that the chieh-qua cultivar 'A39' exhibited simultaneous tolerance to multiple stresses, including heat and drought, which are two negative environmental factors that drastically affect crop growth. Fusarium wilt, as a serious biotic stress, significantly influences chieh-qua productivity and quality. Here, our study mainly focused on exploring genes involved in FA resistance. However, the relationship among different stresses would be necessary for studying the molecular mechanisms of stresses on chieh-qua.
At present, the genome sequence information on wax gourd [2,41] and the RNAseq data of chieh-qua [28,29] provide molecular data for performing genetic research in chieh-qua. Nevertheless, the publicly available data is not sufficient to explain the molecular mechanisms underlying resistance to FW/FA in chieh-qua. Therefore, much more extensive molecular data is needed to discover genes related to FW resistance in chieh-qua. This could serve as a good source for deeply isolating useful genes involved in disease resistance or plant-pathogen interaction. Using transcriptome analysis, multiple DEGs were screened under normal conditions and after FA stress. Among them, several DEGs related to pathogen-related proteins, peroxidase, and ET-responsive transcription factors were possibly involved in FA resistance with significant differences in expression between 'A39' and 'H5'. In view of 'A39' having comprehensive resistance to different stresses, we speculated that common regulatory networks contributing to its greater resistance might exist. For instance, under drought stress [29] and FA treatment, we found that there was a common pathway of plant pathogen interaction by KEGG analysis.

More Apoptosis were Detected in 'H5'
When plants encounter different adversities, several physiological indicators change [3]. For example, MDA is the end-product of membrane lipid peroxidation caused by ROS [42], which could reflect the degree of damage derived from stress. In our study, we found that much more MDA was accumulated in 'H5' than 'A39' after FA stress, indicating the lipid peroxidation of membranes in 'H5' was severe compared with 'A39'. Meanwhile, by TUNEL assay, we found that more cells in 'H5' were apoptotic with stronger signals than 'A39', suggesting 'H5' was more affected by FA stress.

Analysis of Pathogen-Related (PR) Protein under FA Stress
For defending against pathogens, plants have evolved a multiple of internal strategies including pathogenesis-related (PR) proteins [43], WRKY TFs [44], and plant hormones [15]. Thus, getting an understanding of the complexity and regulation of FW resistance could promote the development and breeding for disease resistance of chieh-qua.
PR proteins accumulated upon pathological conditions in at least two or more plantpathogen combinations [15,45]. Several external and internal factors could induce PR expression. PR10 was up-regulated when plants were exposed to salicylic acid (SA), abscisic acid (ABA), gibberellin (GA 3 ), and H 2 O 2 in Solanum capsicoides [10], which is also suitable to pepper (Capsicum annuum L.) [46], and rice (Oryza sativa L.) [47]. Another PR protein, PRB1-2 was highly induced in leaf tissues in peanut (Arachis hypogaea Linn.) [48]. Our previous study showed that the gene PRB1-2-like (pathogenesis-related protein) was significantly up-regulated in 'H5' in contrast to 'A39' under drought stress [29]. Here, we detected five PR proteins including PRB1-2 and PR-4A with up-regulated expression under FA stress, especially in 'A39', indicating the transcription level of the PR protein could be induced by FA and increased PR transcripts might partly contribute to plant disease resistance in chieh-qua.

Analysis of Genes Related to Enzymes under FA Stress
The generation of ROS when plants are exposed to oxidative stresses and could be eliminated by enzymes (peroxidase, SOD, and catalase) and non-enzymatic pathways [3,49]. Previous studies identified that improving ROS scavenging capacity could significantly enhance plant stress tolerance [50,51]. Peroxidase has a central role in the control of ROS, and the H 2 O 2 could be highly induced after pathogen infection [52]. When Arabidopsis was exposed to a pathogen, a peroxidase-dependent apoplastic oxidative burst occurred to improve resistance to the disease [53]. A class III peroxidase was specifically expressed in pathogen-attacked barley epidermis, leading to basal resistance [54]. In our study, 'A39' exhibited higher SOD enzyme activity than 'H5' and two peroxidase genes were up-regulated in 'A39' in comparison with 'H5', indicating 'A39' exhibited stronger ROS scavenging ability compared with 'H5'. This result suggested that the increased antioxidant capacity might contribute to the enhanced FA stress in chieh-qua.

Analysis of Genes Related to Hormones under FA Stress
Phytohormones play important roles in regulating plant responses to a wide range of biotic and abiotic stresses [55][56][57]. Several host reactions to pathogen infections are influenced by phytohormones such as ET and SA [24,58,59]. Enhanced ET production could induce defense-related effector genes through a cascade of events [60]. It is also a critical component of responses to herbivory, pathogen attack, and mechanical damage [60]. Overexpression of OsACS2 (a key enzyme of ET biosynthesis) significantly increased the level of endogenous ET and related gene expression, especially in response to pathogen (Magnaporthe oryzae and Rhizoctonia solani) infection [61]. A wheat ET response factor mediated host responses to a necrotrophic pathogen [62]. The AtPDF1.2 encoding the plant defense gene is widely considerd a marker of ET-induced signaling in Arabidopsis defense responses [63,64]. As the gene contains GCC box promoter elements, it is inducible by both ET and JA through activation of AtERF1 [63,65]. In this study, the KEGG analysis of DEGs revealed the enriched pathways for plant hormones signal transduction, which included nine DEGs related to SA, ET, and GA signaling pathways. Among them, genes related to the ET signaling pathway were only up-regulated in 'A39' compared with 'H5'. This indicated that the activation of the ET signal pathway might perform as an indirect mechanism for enhanced FA resistance in chieh-qua.
Taken together, this study revealed that 'A39' exhibited enhanced FA resistance compared with 'H5'. The compared global transcriptome analysis allowed us to put forward a potential molecular pathway to explain the increased FA stress in 'A39' (Figure 8). When exposed to FA stress, the up-regulated genes of PR protein might cause PR response and the activation of ET-related genes may lead to efficient ET signal transduction and synthesis. In addition, the increased expression of peroxidase and SOD activity in 'A39' may lead to enhanced ROS scavenging ability compared with 'H5'. In conclusion, the study revealed the DEGs expression changes underlying FA stress resistance of 'A39' and will help illuminate the mechanisms and isolate genes in response to abiotic stresses.

PEER REVIEW 11 of 14
bidopsis defense responses [63,64]. As the gene contains GCC box promoter elements, it is inducible by both ET and JA through activation of AtERF1 [63,65]. In this study, the KEGG analysis of DEGs revealed the enriched pathways for plant hormones signal transduction, which included nine DEGs related to SA, ET, and GA signaling pathways. Among them, genes related to the ET signaling pathway were only up-regulated in 'A39′ compared with 'H5′. This indicated that the activation of the ET signal pathway might perform as an indirect mechanism for enhanced FA resistance in chieh-qua. Taken together, this study revealed that 'A39′ exhibited enhanced FA resistance compared with 'H5′. The compared global transcriptome analysis allowed us to put forward a potential molecular pathway to explain the increased FA stress in 'A39′ ( Figure  8). When exposed to FA stress, the up-regulated genes of PR protein might cause PR response and the activation of ET-related genes may lead to efficient ET signal transduction and synthesis. In addition, the increased expression of peroxidase and SOD activity in 'A39′ may lead to enhanced ROS scavenging ability compared with 'H5′. In conclusion, the study revealed the DEGs expression changes underlying FA stress resistance of 'A39′ and will help illuminate the mechanisms and isolate genes in response to abiotic stresses. Figure 8. A model for the mechanism of FA resistance in chieh-qua. Under FA stress, 'A39′ exhibited enhanced resistance compared with 'H5′ due to induction of genes encoding pathogen-related protein and ET biosynthesis and signal transduction. Also, increased expression of peroxidase and SOD activity would improve ROS scavenging ability, reducing stress effects.

Conclusions
FA produced by pathogens of the genus Fusarium could lead to negative effects on plant development. The mechanism of plant resistance to FA stress remains unclear. In our study, we compared the transcriptome results between 'A39′ and 'H5′ under normal and FA stress in chieh-qua. Results showed that genes related to pathogen-related protein and ethylene were most related to chieh-qua FA resistance. Our data provides a basis for FA resistance gene isolation and be helpful in the molecular mechanism of FA resistance and in molecular breeding of Fusarium wilt resistance in chieh-qua.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: qRT-PCR validation of DEGs under FA stress. Correlation between the fold change analyzed by Figure 8. A model for the mechanism of FA resistance in chieh-qua. Under FA stress, 'A39' exhibited enhanced resistance compared with 'H5' due to induction of genes encoding pathogen-related protein and ET biosynthesis and signal transduction. Also, increased expression of peroxidase and SOD activity would improve ROS scavenging ability, reducing stress effects.

Conclusions
FA produced by pathogens of the genus Fusarium could lead to negative effects on plant development. The mechanism of plant resistance to FA stress remains unclear. In our study, we compared the transcriptome results between 'A39' and 'H5' under normal and FA stress in chieh-qua. Results showed that genes related to pathogen-related protein and ethylene were most related to chieh-qua FA resistance. Our data provides a basis for FA resistance gene isolation and be helpful in the molecular mechanism of FA resistance and in molecular breeding of Fusarium wilt resistance in chieh-qua.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/horticulturae7040088/s1, Figure S1: qRT-PCR validation of DEGs under FA stress. Correlation between the fold change analyzed by RNA-seq (y-axis) and data obtained using qRT-PCR (x-axis). Figure S2: SNP analysis of A39 and H5 under normal and FA stress. Figure S3: GO analysis of A39 and H5 under normal condition and FA stress. Table S1: Primers used in qRT-PCR, Table S2: Mapping results of RNA-seq reads of A39 and H5 under Normal condition (N) and FA stress (F). Table S3: DEGs between A39 and H5 before FA stress treatment. Table S4: DEGs between A39 and H5 under FA stress treatment. Table S5: Genes of qRT-PCR identification. Tables S6 and S7: SNP variation and detection analysis. Table S8: Gene Ontology (GO) analysis.   Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: All the sequencing data were submitted to NCBI Sequence Read Archive database under accession number SRP180414.

Conflicts of Interest:
All authors declare that the research was conducted in the absence of any commercial or financial relationships which could be construed as potential conflicts of interest.