Integrated Transcriptomic and Un-Targeted Metabolomics Analysis Reveals Mulberry Fruit (Morus atropurpurea) in Response to Sclerotiniose Pathogen Ciboria shiraiana Infection

Mulberry sclerotiniose caused by Ciboria shiraiana is a devastating disease of mulberry (Morus alba L.) fruit in Northwest China. At present, no disease-resistant varieties are used in production, as the molecular mechanisms of this disease are not well understood. In this study, to explore new prevention methods and provide direction for molecular breeding, transcriptomic sequencing and un-targeted metabolomics were performed on healthy (CK), early-stage diseased (HB1), and middle-stage diseased (HB2) mulberry fruits. Functional annotation, gene ontology, a Kyoto encyclopedia of genes and genomes (KEGG) analysis, and a Mapman analysis of the differentially expressed genes revealed differential regulation of genes related to plant hormone signal transduction, transcription factors, and phenylpropanoid biosynthesis. A correspondence between the transcript pattern and metabolite profile was observed in the phenylpropanoid biosynthesis pathway. It should be noted that the log2 ratio of eugenol (isoeugenol) in HB1 and HB2 are 85 times and 23 times higher than CK, respectively. Our study shows that phenylpropanoid biosynthesis may play an essential role in response to sclerotiniose pathogen infection and eugenol(isoeugenol) enrichment in mulberry fruit, which may provide a novel method for mulberry sclerotiniose control.


Introduction
Mulberry fruit is a valuable food because of its high nutritional value, pleasant taste, and large amount of biologically active components that might be associated with some potential pharmacological functions that are beneficial for human health [1,2]. In recent years, with both ecological and economic benefits, large-scale mulberry orchards have been set up and the mulberry fruit industry has developed rapidly in China, becoming a model for the diversified development of mulberry in Northwest China. However, incidences of disease in the fruit have also increased. Mulberry sclerotiniose is a devastating fungal disease of mulberry fruit; it not only causes mulberry fruit to lose their color and turn white, but also heavily affects mulberry fruit quality and yield.
Mulberry sclerotiniose can be classified into hypertrophy sorosis sclerotiniose, reduced sorosis sclerotiniose, and small particles sorosis sclerotiniose, according to its morphology and symptoms. Hypertrophy sorosis sclerotiniose is a common and major disease in Northwest China, and its pathogen is Ciboria shiraiana or Sclerotinia sclerotinia. The diseased fruit fall onto the ground and the pathogen overwinter in the soil as sclerotium. During the mulberry flowering time of the following spring,

Symptoms of Infected Mulberry Fruits and Fungal Detection
Healthy mulberry fruits were collected from a normal mulberry orchard and used as the control (CK). Diseased mulberry fruits were divided into three classes: HB1, HB2, and HB3. In the early stage of infection, only a few achenes turned white or gray at stage HB1. The ovary of the mulberry flower swelled due to hyphal proliferation, and the whole fruit was malformed at stage HB2. At stage HB3, the ovary swelled further, some of the diseased fruit turned pink, and the sclerotium was forming ( Figure 1A). This was the typical morphological characteristics of hypertrophy sorosis sclerotiniose. In addition, a lot of different levels of diseased mulberry fruits could be found on infected mulberry trees ( Figure S1A) and a large number of fungal fruiting bodies were discovered on the ground (Figure S1B,C). Plenty of matured ascospores existed in the apothecium of fruiting bodies ( Figure S1D). Moreover, the pathogen in diseased fruits was identified as C. shiraiana by fungal ITS primers and a phylogenetic tree analysis, while no pathogens were detected in mulberry fruit from a healthy mulberry orchard ( Figure S2A,B). The hyphal of C. Shiraiana that had proliferated could be observed by Calcofluor White fluorescent staining ( Figure 1B) and pectinase activity increased in diseased fruits ( Figure S2C). The plant tissue degraded rapidly, in keeping with hyphal proliferation, and mulberry fruit DNA was barely detected at stage HB3. Therefore, HB1 and HB2 were chosen for subsequent transcriptomic and metabolomic analyses.

Identified Differentially Expressed Genes (DEGs), Their Gene Ontology (GO) Enrichments, and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichments by Transcriptome Analysis
In order to compare the difference in the transcriptome of mulberry fruit after fungal infection, we performed RNA-seq profiling with the samples CK, HB1, and HB2. A total of 5109 (with 1926 upand 3183 down-regulated) and 6760 (with 2655 up-and 4105 down-regulated) differentially expressed genes (DEGs) were identified in infection stages HB1 and HB2, respectively, as compared with CK ( Figure 2A). As shown in the Venn diagram in Figure 2B, 3315 genes of all DEGs were regulated by C. shiraiana at both HB1 and HB2. These results show the differences at the molecular level between healthy fruits and diseased fruits, and also between different disease levels at HB1 and HB2. A GO analysis was carried out to classify the functions of the DEGs of HB1 and HB2 compared with healthy mulberry fruit CK. Within the biological process categories, The DEGs from the comparisons of HB1 and HB2 were significantly enriched in the terms of cellular processes, metabolic processes, biological regulation, and response to stimuli. Within the cellular component category, the membrane, membrane parts, cells, and organelles were significantly enriched. Within the molecular function category, the DEGs reveal significant enrichments in categories for binding and catalytic activity.
The results indicate that the DEGs were functionally enriched in the same GO terms for the top 10 categories (Figure 3, Table S1). However, further study of GO enrichment showed that the DEGs of HB1 and HB2 were different ( Figure S3). KEGG pathway analysis of the DEGs identified the terms global and overview maps, carbohydrate metabolism, biosynthesis of other secondary metabolites, signal transduction, amino acid metabolism, and lipid metabolism as being significantly enriched as compared with the whole genomic background in both HB1 and HB2 ( Figure S4). KEGG function annotations were used to analyze the KEGG enrichment results as a scatter plot. In both HB1 and HB2, the top three highest significantly enriched KEGG pathways were flavonoid biosynthesis, plant hormone signal transduction, and phenylpropanoid biosynthesis ( Figure 4). To further validate the expression profiles of the RNA-Seq data, eight transcripts from top three of the highest significant terms were selected for analysis using quantification real-time PCR (qRT-PCR). The results from the qRT-PCR analysis are generally in agreement with the expression profiles obtained using the RNA-Seq data ( Figure S5, Table S2).

Transcriptomic Response of Mulberry Fruits Infected by C. shiraiana
According to the GO and KEGG enrichment results, plant hormone signal transduction, phenylpropanoid biosynthesis, and transcription factors were analyzed to help understand the response of mulberry fruits to fungal pathogen C. shiraiana infection.

Plant Hormone Signal Transduction Pathway
A group of genes changed involved in hormone biosynthesis, mobilization, and signal transduction in response to C. shiraiana infection. Transcripts related to plant hormone signal transduction pathway showed that the most were involved in auxin, brassinosteroid, cytokinine, jasmonic acid, gibberellin, abscisic acid, ethylene, salicylic acid, and so on. Log fold ratios for differentially expression genes in HB1 and HB2 compared with CK were shown ( Table 1). The transcripts' up-regulated or down-regulated trends in HB1 and HB2 were similar. Transcripts for abscisic acid receptor (PY1 and PY4), involved in the abscisic acid pathway, were more abundant in HB1 and HB2. Several transcripts related to jasmonic acid and brassinosteroid pathway were down-regulated in HB1 and HB2, including a large number of TIFY protein and protein kinase. Of auxin biosynthesis and signal transduction, transcripts for auxin-responsive protein IAA21 and SAUR36 were down-regulated, while IAA4, IAA29, SAUR24 and SAUR50 were up-regulated. Their detailed roles involved in the responses to C. shiraiana infection remain to be explored. Table 1. Differential expressed genes (DEGs) related to the plant hormone signal transduction pathway. Log-fold ratio is shown as a gradient between red (up-regulated) and green (down-regulated).

Phenylpropanoid Biosynthesis Pathway
The phenylpropanoid biosynthesis pathway plays an important role in the regulation of secondary metabolites in plants. In this study, 40 DEGs were detected in both HB1 and HB2 ( Table 2). Most of the DEGs were up-regulated in HB1 and HB2 compared with CK, such as phenylalanine ammonia-lyase, eugenol synthase, cinnamyl alcohol dehydrogenase, and cinnamoyl-CoA reductase. However, β-glucosidase was down-regulated. It should be noted that both up-regulated and down-regulated genes existed in peroxidase, caffeic acid 3-O-methyltransferase, and 4-coumarate-CoA ligase. Table 2. Classification of the differentially expressed genes identified from HB1 and HB2 compared with CK in the phenylpropanoid biosynthesis pathway. Log-fold ratio is shown as a gradient between red (up-regulated) and green (down-regulated).

Transcription Factors
In the present study, 1345 transcription factors (TFs) genes were detected, 306 and 443 TFs were differentially expressed in HB1 and HB2, respectively. There were 122 up-regulated and 184 down-regulated TFs in HB1, and 164 up-regulated and 279 down-regulated TFs in HB2 (Table 3). Of all the differentially expressed TFs, the MYB cluster showed the largest percentage, followed by AP2-EREBP, bHLH, NAC, WRKY, C2H2, MADS, etc. Most of the TFs in HB2, both up-and downregulated, were increased compared with HB1. Specially, some TFs were down-regulated only in HB2, such as SBP, GRF, C3H, TCP and Tify.

Overview of the Metabolites in the Phenylpropanoid Pathway by Metabolimics Data
In order to obtain an overview of metabolic changes in response to C. shiraiana of mulberry fruits, comparative metabolite analyses were performed at CK, HB1, and HB2 using an untargeted metabolomics approach. A principal component analysis (PCA), and a partial least squares discriminant analysis (PLS-DA) were used to analyze these data (Figures S6 and S7). Evaluation parameters of PLS-DA are shown in Table S3. The results show that the three samples were easy to distinguish from the data and the model was stable and reliable.

Discussion
In order to maintain predominance, plants are engaged in a co-evolutionary fight against their pathogens [17]. The widespread changes in plants during pathogen infection make understanding the interactions a huge task [18]. To have a global view of all DEGs in mulberry fruits at HB1 and HB2, a large number of genes involved in metabolism were assessed by a Mapman analysis. Mapman is a common software for displaying large data sets onto pictorial diagrams that symbolically depict areas of biological function [19]. In this study, mainly DEGs of plant hormone signal transduction, phenylpropanoid pathway and transcription factors were discussed. Figure 5 summarizes the global view of DEGs in response to C. shiraiana infection in both HB1 and HB2 using Mapman analysis. Plant hormones play key roles in the regulation of immune responses to microbial pathogens [20][21][22][23]. Jasmonic acid and salicylic acid coordinate growth and defense response upon fungal infection in poplar [24]. Jasmonic acid and terpenoids play roles in resistance against Phytophthora cinnamomic in Zea mays roots [25]. Ethylene plays a positive role in adapting mulberry to cope with the drought and salinity condition [26]. Transcription factors are also an important factor of gene expression change [27]. PbMYB10b and PbMYB9 regulate anthocyanins and flavonols of the flavonoid biosynthesis in pear fruit [28]. The bHLH factor VvMYC1 regulates anthocyanin and/or proanthocyanidin (PA) synthesis of the flavonoid pathway in grapevine [29]. In this study, the down-regulation of most MYB and bHLH clusters in HB1 and HB2, which indicates that they may affect the mulberry flavonoid biosynthesis. As shown in Figure 6, in the case of pathogen infection, primary metabolism, such as raffinose synthesis, lipid synthesis, hemicellulose synthesis, cellulose synthesis, and cell wall proteins were decreased, while cell wall modification and degradation, and lipid degradation increased in plants. Light reactions, the Calvin cycle, and glycolysis were also indicated to be suppressed in plants. Changes in plant secondary metabolites were more obvious. Terpenes, phenols, and waxy compounds were suppressed, and phenylpropanoids and flavonoids were increased (Figures 5 and 6). The Mapman analysis results show that most of genes related to secondary metabolism had decreased expression levels, which is similar to the results of the GO and KEGG enrichment analyses. With the infection of pathogens, the infection hyphae rapidly expanded, which took away a lot of plant nutrients, resulting in the inhibition of plant metabolism and the destruction of plant tissue. These also appear during infection by other pathogens. In rice blast disease, the photoassimilates were shifted to mannitol and glycerol for mycelial growth [30]. The bean metabolites that varied in leaves included amines/amino acids, organic acids, phytoalexins, and ureides during Sclerotinia sclerotiorum infection [31]. Terpene down-regulation triggers pathogen defense responses in transgenic orange [32]. Figure 6. Metabolism overview in Mapman depicting differential regulated genes in HB1 and HB2. Log-fold ratio is shown as a gradient between red (up-regulated) and green (down-regulated).
Today, numerous questions about the mechanism of pathogen infection still need to be addressed. Metabolomics is an important platform for studying stress in plant-fungal pathogen interactions [25,33]. We can qualify and quantify the metabolites of plants with biotic stress using modern analytical techniques. When pathogens infect, a large number of plant secondary metabolites are found to be up-regulated. In barley, higher levels of flavonoids, phenylpropanoids, and metabolites of fatty acids and terpenoid pathways were found in the resistant barley lines compared with the susceptible lines during infection with Fusarium [34]. Phenylpropanoid pathway intermediates, like 4-hydroxybenzoate, cinnamic acid, ferulic acid, and caffeic acid, were highly accumulated in the resistant soybean line against Sclerotinia sclerotiorum infection [35]. Moreover, in trees displaying symptoms of wood decay, large amounts of lignans were observed [36]. Mulberry trees contain large quantities of flavonoids and other secondary metabolites, such as phenylpropanoids and alkaloids, which have wide roles in defenses against biotic and abiotic stresses [26] and have long been used in Chinese traditional medicine for human disease treatment [37]. A large number of flavonoid metabolites from mulberry fruit have been evaluated both in vitro and in vivo for their antioxidant activity [38], such as morcin N [39], anthocyanins [40], and so on. Combined with transcriptome and metabolite profiling data, the significant genes and metabolites were found by co-expression analysis of the phenylpropanoid biosynthesis pathway of mulberry fruit in response to C. shiraiana (Figure 7, Table S5). We found that a large number of phenylpropanoid pathway intermediates, cinnamic acid, coumaric acid, caffeic acid, caffeyl aldehyde, coniferyl aldehyde were decreased. However, the final metabolites, like chavicol, isochavicol, methychavicol, anethole, eugenol, isoeugenol, methyleugenol, isomethyleugenol, and sinapic acid were increased. In particular, the large increase in eugenol/isoeugenol levels indicates that this compound may play an important role in response to pathogen infection. Eugenol is widely present in plants and is mainly used for antibacterial purposes, and can also be used in various cosmetic flavors and soap flavor formulations. There is no clear evidence that eugenol is carcinogenic. Eugenol exhibits an antimicrobial effect against carbapenem-resistant Klebsiella pneumoniae (CRKP) strains and could be potentially used to control CRKP-related infections [41]. The antimicrobial sachet containing the synergistic inhibitory effect of eugenol and citral (SEC) was developed for bread preservation [42]. Eugenol-casein nanoparticles (EC-NPs) exhibited a great antifungal activity against spore germination of fungus and have shown potential as an environmentally friendly preservative in the food industry [43]. Our data reveal important roles of eugenol/isoeugenol enrichment in mulberry fruit in response to C. shiraiana, which may provide a novel method for mulberry sclerotiniose control.

Plant Growth Conditions and Fungus Identification
Mulberry cultivar Hongguo2 (Morus atropurpurea), which is a widely cultivated mulberry fruit tree in Northwest China, was bred by the Institute of Sericulture and Silk, Northwest A&F University.
A four-year high-yielding fruit mulberry orchard with area of about 3000 m 2 , and a space of 3 m between the rows and 1.5 m between the plants, was located in Guangji Town, Zhouzhi County, Shannxi Province, China. The climate and soil conditions of the orchard were consistent. A total of 20 mulberry trees suffered a serious outbreak of mulberry sclerotiniose on one side of the mulberry orchard near the low slope. A large number of diseased fruits dropped into the ground and formed sclerotia later. About 30 sclerotia per tree were discovered in winter, and PCR amplification was done using ITS sequences as primers, C. shiraiana was identified by morphological characters and phylogenetic tree. The diseased fruit could be divided into three stages according to the degree of the disease. Only a few achenes turning light white from green were marked as HB1, and a few appearing white and swelling slightly marked as HB2, while most achenes of mulberry fruit turned white and the whole fruit showed malformation, being marked as HB3. On the other side of mulberry orchard, about 200 m from the diseased area, mulberry fruits grew healthy and no pathogens in the fruit were detected by ITS sequencing, and marked as CK.
HB1, HB2, HB3, CK, fungal fruiting bodies, and sclerotia were collected at the same time. The samples were frozen in liquid nitrogen for storage at −80 • C. Three biological replicates were analyzed using transcriptome sequencing, qRT-PCR validation, and the enzyme activity test. Six biological replicates were used for metabolite profiling. Each replicate consisted of three single fruit, and a single fruit collected from an independent tree.

Pathogen Detection and Enzyme Activity Analysis
Genome DNA of CK, HB1, HB2, and HB3 were extracted from 0.1 g tissue using a rapid plant genomic DNA isolation kit (Sangon, Shanghai, China). Fungal genome DNA were extracted with a rapid fungi genomic DNA isolation kit (Sangon, Shanghai, China), and PCR amplification using fungal universal ITS1 and ITS4 primers.
Pectinase activity was measured by a pectinase assay kit (Jiancheng Biotechnology, Nanjing, China). The activity of the pectinase enzyme was expressed as units per milligrams of protein.
Statistical differences were based on one-way ANOVA analysis followed by Tukey's HSD t-test using SPSS 19.

Tissue Paraffin Section and Fluorescent Staning
Several samples were fixed in 4% buffered formalin, embedded in paraffin, and sectioned. For toluidine blue staining, the sections were cut and stained with toluidine blue solution (0.5% toluidine blue, 1% borax). For lactophenol-cotton blue staining, the staining included 0.05 g cotton blue, 20 g phenol crystals, 40 mL glycerol, 20 mL lactic acid, and 20 mL distilled water. For Calcofluor White Staining, one drop of Calcofluor White Stain (18909, Sigma-Aldrich, St. Louis, MO, USA) and one drop of 10% potassium hydroxide were added to the sections, a coverslip was placed over the specimen, let stand for 1 min, and finally examined under UV light (355 nm) [44].

RNA Extraction, RNA Sequencing, and Functional Annotations
Total RNA in each sample was extracted using the MiniBEST Plant RNA Extraction Kit (Takara Bio Inc., Dalian, China) by following the manufacturer's procedure. An ethanol precipitation protocol and CTAB-PBIOZOL reagent was used for the purification of total RNA. The total RNA quantity and purity were analyzed by a Nano Drop and Agilent 2100 bioanalyzer (Thermo Fisher Scientific, Massachusetts, USA). Oligo (dT)-attached magnetic beads were used to purify mRNA. Sequencing was carried out using the BGISEQ-500 platform (BGI, Shenzhen, China).
The raw reads were filtered by removing adaptor sequences and low-quality sequences which contain more than 5% unknown nucleotides (N), and low-quality reads (those with >20% of bases with a quality value of <15) with SOAPnuke v1.5.2 [45]. The clean reads were aligned to the reference index by HISAT2 v2.0.4 [46]. String Tie (http://ccb.jhu.edu/software/stringtie) v1.0.4 was used to assemble the transcripts depending on the mulberry (Morus notabilis) annotation. The gene expression level was calculated with fragments per kilobase of transcript per million (FPKM) mapped reads. FPKM > 1, Fold Change > 2, and Adjusted p-value < 0.001 were used as the thresholds for the differential gene expression (DEGs) for the comparisons of HB1 vs. CK and HB2 vs. CK, and enrichment analyses of the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of all the DEGs were analyzed.

Validation of Gene Expression by Quantitative Real-Time PCR (qRT-PCR)
To confirm the RNA-Seq results, 10 genes were selected to examine the consistency of their expression using qRT-PCR assays. All of the gene-specific primers were designed by Primer-BLAST (https://blast.ncbi.nlm.nih.gov/blast.cgi), and the primers are listed in Table S2. PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Dalian, China) was used for cDNA synthesis. The ABI 7500 Fast Real-Time PCR Systems (Applied Biosystems, Waltham, Massachusetts, USA) and the abm ® EvaGreen qPCR MasterMix-No dye kit (Applied Biological Materials Inc., Vancouver, BC, Canada) were used for all of the qRT-PCR assays. The 2 −∆∆Ct method was used to calculate the relative gene expression based on the qRT-PCR data [47].

Metabolite Profiling Assay by Ultra-Performance Liquid Chromatography (UPLC) and Mass Spectrometry (MS)
All samples were prepared according to the LC-MS system machine orders. All chromatographic contents were performed using an ultra-performance liquid chromatography (UPLC) system (Waters, Manchester, UK), and ACQUITY UPLC HSS T3 column (100 mm × 2.1 mm, 1.8 µm, Waters, UK) was used for the reversed phase separation. The mobile phase composed of solvent A (water, 0.1% formic acid) and solvent B (acetonitrile, 0.1% formic acid), and the flow rate of the mobile phases was 0.4 mL/min. Gradient elution conditions were set as follows: 0~2 min, 100% A; 2~11 min, 0% to 100% B; 11~13 min, 100% B; 13~15 min, 0% to 100% A. The column oven was maintained at 50 • C and the injection volume for each sample was 10 µL.
A high-resolution tandem mass spectrometer Xevo G2 XS QTOF (Waters, Manchester, UK) was used to test metabolites eluted from the column and operated in both positive and negative ion modes. The capillary and sampling cone voltages were set at 3.0 kV and 40.0 V, respectively, in positive ion mode, while 2.0 kV and 40.0 V were set in negative ion mode. The TOF mass range was between 50 and 1200 Da, and the scan time was set as 0.2 s. For the MS/MS detection, all precursors were fragmented using 20-40 eV, and the scan time was 0.2 s. In the data acquisition process, real-time quality correction was performed on the LE signal every 3 s. At the same time, quality control samples were collected every 10 samples to evaluate the stability of the instrument state during sample collection.

Metabolite Data Processing and Analysis
The raw data, both ESI negative and positive, were pretreated by Progenesis QI v 2.2, (Waters, UK) and metaX [48] (BGI, Shenzhen, China). Principal component analysis (PCA) and partial least square discriminant analysis (PLS-DA) were used for comparing the groups among CK, HB1, and HB2. The metabolites were initially identified based on their retention time, molecular weight, fragmentation patterns, and ultraviolet absorption, compared to standards and published references. The variable importance in the projection (VIP) value and ANOVA methods, under the condition of VIP > 1.0, p < 0.05, and fold change < 0.5 or > 1.5, were used for preliminary assessment of significantly different metabolites identified from UPLC-QTOF-MS analysis.

Mapman Analysis
Protein functions were categorized using MapMan software (available online: http://mapman. gabipd.org/). Annotations were transferred to the Arabidopsis genome with consideration for orthologous genes to predict functions of identified proteins from M. notabilis.