Combined Transcriptome and Metabolome Analysis Reveals Adaptive Defense Responses to DON Induction in Potato

Phytophthora infestans poses a serious threat to potato production, storage, and processing. Understanding plant immunity triggered by fungal elicitors is important for the effective control of plant diseases. However, the role of the potato stress response to Fusarium toxin deoxynivalenol (DON)-induced stress is still not fully understood. In this study, the metabolites of DON-treated potato tubers were studied for four time intervals using UPLC-MS/MS. We identified 676 metabolites, and differential accumulation metabolite analysis showed that alkaloids, phenolic acids, and flavonoids were the major differential metabolites that directly determined defense response. Transcriptome data showed that differentially expressed genes (DEGs) were significantly enriched in phenylpropane and flavonoid metabolic pathways. Weighted gene co-expression network analysis (WGCNA) identified many hub genes, some of which modulate plant immune responses. This study is important for understanding the metabolic changes, transcriptional regulation, and physiological responses of active and signaling substances during DON induction, and it will help to design defense strategies against Phytophthora infestans in potato.


Introduction
To date, more than 10,000 varieties of potato (Solanum tuberosum L.) have been cultivated in more than 160 countries and regions, making it the world's most important non-cereal food crop [1,2]. Its tubers are rich in starch, dietary fiber, protein, minerals, vitamins, and other nutrients required by the human body [3]. As a high-yielding crop, potatoes are widely used as food ingredients, feed processing, bioenergy reserves, and in environmental material development, playing an important role in ensuring food security and increasing income. Global losses due to Fusarium infection account for 17.2% of total potato production, causing huge economic losses [4]. At moderate temperatures (16)(17)(18)(19)(20)(21)(22) • C) and high humidity (over 97%), the pathogen usually infects tubers by being washed down through soil [5]. Potatoes are susceptible to infection during planting and storage, resulting in red, brown or granular rot of the tubers, which destroys the entire crop within days of infection, severely affecting potato quality and yield. Storage rot can also occur if the pathogen is present, and infections can lead to secondary infections of fungi and bacteria. Although chemical control can play a better role in a short period of time, there are problems with resistance and pesticide residues. Biological control has the advantages of being environmentally friendly and less prone to resistance during disease prevention and control, and it is an important complementary means of chemical control. For example, tuber biocontrol of potato late blight was achieved using phenazine-1-carboxylic acid produced by Pseudomonas spp. [6]. Studies have shown that P. infestans also produces virulence factors (cell wall-degrading enzymes, virulence-related effector proteins, and In recent years, combined transcriptional and metabolic analyses have become an important tool for discovering new signaling pathways and responses to environmental stress induction, and they are an effective method for revealing plant stress responses, plantpathogen interactions, mechanisms of metabolite synthesis, and identification of resistance genes [39]. For example, transcriptome analysis of P. infestans infected potatoes revealed SA-JA-ABA(abscisic acid)-mediated signaling pathways [40]; Li et al. [41] elucidated the molecular mechanisms of metabolite accumulation and flavonoid biosynthesis in jujube leaves by metabolomic and transcriptomic analyses; and Yang et al. [42] elucidated the biosynthetic regulatory network of flavonoid metabolites in Salvia miltiorrhiza stems and leaves through transcriptional metabolic analysis.
Our previous study found that treatment of potato tubers with low concentrations of DON also affected resistance to late blight caused by P. infestans [43]. However, the effects of fungal toxins on potato tuber resistance and the relationship between resistance and transcriptional and metabolic changes have not been reported, and the molecular basis and physiological and biochemical changes in potato-induced resistance under DON induction are not clear. Therefore, in this study, RNA-seq technology and non-targeted metabolomics were used to study the expression patterns of genes and overall changes in metabolites under DON stress using Atlantic potato tuber as material. Four time points were taken under DON treatment at 5 ng/mL to reveal the specific genes regulating the resistance response under different treatment times as well as the transcriptional network and major metabolic pathways in response to DON stress in order to further understand three things-the mechanism of DON stress response in potato, the DEGs in the transcriptional network, and major metabolic pathways in response to DON stress.

UPLC-MS/MS-Based Quantitative Metabolomic Analysis of the Atlantic Potato Cultivar Tuber
To better understand the dynamics of relevant metabolites during DON stress, total metabolites were extracted from each sample, and metabolite detection was performed using UPLC-MS/MS. The mass spectrometry data of each group of samples (CK0h, DON0h, DON4h, DON12h, and DON48h) were processed using Analyst 1.6.3 software, and a total of 676 metabolites were detected based on the local metabolic database. The integration of the peaks and the correction of the chromatographic peaks were performed on the downstream mass spectrometry results of each sample using MultiaQuant software (v3.0.3) to obtain the relative content of each substance. The 15 samples were clearly divided into 5 groups by heat map and correlation analysis ( Figure 1A), and there was a significant correlation between the samples within the groups ( Figure S1). Additionally, PCA showed clear separation between the groups ( Figure 1B). As shown in Figure 1C, 676 metabolites could be divided into 13 different categories, but most were concentrated in lipids (17%), flavonoids (12%), alkaloids (12%), amino acids and derivatives (12%), organic acids (9%), and others. Overall, the results indicate that DON stress in potato has different metabolite distributions at different times.

Metabolomic Changes at Different Time Points
To show the changes in the metabolite content of potato tubers under different DON stress durations, we used orthogonal partial least squares discriminant analysis (OPLS-DA) to screen for differential metabolites, and metabolites with fold change ≥2 or ≤0.5, along with VIP ≥ 1, were considered as significantly differential metabolites. Compared with CK0h, 75, 117, 180, and 197 differential metabolites were identified in DON0h, DON4h, DON12h, and DON48h, respectively, and more upregulated metabolites were specific at 48 h ( Figure 2). Among them, six metabolites were upregulated and 69 metabolites were downregulated in the DON0h treatment compared with the control. In DON4h, 75 and 42 metabolites were upregulated and downregulated, respectively. A total of 102 and 78 metabolites were upregulated and downregulated, respectively, in DON12h. In addition, 144 upregulated metabolites and 55 downregulated metabolites were identified in DON48h.

Metabolomic Changes at Different Time Points
To show the changes in the metabolite content of potato tubers under different DON stress durations, we used orthogonal partial least squares discriminant analysis (OPLS-DA) to screen for differential metabolites, and metabolites with fold change ≥2 or ≤0.5, along with VIP ≥ 1, were considered as significantly differential metabolites. Compared with CK0h, 75, 117, 180, and 197 differential metabolites were identified in DON0h, DON4h, DON12h, and DON48h, respectively, and more upregulated metabolites were specific at 48 h ( Figure 2). Among them, six metabolites were upregulated and 69 metabolites were downregulated in the DON0h treatment compared with the control. In DON4h, 75 and 42 metabolites were upregulated and downregulated, respectively. A total of 102 and 78 metabolites were upregulated and downregulated, respectively, in DON12h. In addition, 144 upregulated metabolites and 55 downregulated metabolites were identified in DON48h.  To further understand the in vivo interactions of the differential metabolites in plants, we annotated and enriched the functions of the differential metabolites based on the KEGG database and found that all the differential metabolites were mainly annotated in metabolic pathways and biosynthesis of secondary metabolites. The enrichment results showed that nine pathways-metabolic pathways, biosynthesis of secondary metabolites, ABC transporters, biosynthesis of amino acids, 2-oxocarboxylic acid metabolism, aminoacyl-tRNA biosynthesis, phenylpropanoid biosynthesis, arginine and proline metabolism, and cysteine and methionine metabolism-were enriched at various time periods under DON stress ( Figure 3C).
Taken together, our results indicate that potatoes under DON stress mainly exhibit phenolic acid and alkaloid metabolite accumulation patterns, with significant differences in defense responses for approximately 50% of the metabolites that are directly or indirectly attributed to the phenylpropane metabolic pathway. metabolism, and cysteine and methionine metabolism-were enriched at various time periods under DON stress ( Figure 3C). Taken together, our results indicate that potatoes under DON stress mainly exhibit phenolic acid and alkaloid metabolite accumulation patterns, with significant differences in defense responses for approximately 50% of the metabolites that are directly or indirectly attributed to the phenylpropane metabolic pathway.

Transcriptional Analysis of Potato Tubers under DON Stress
To study the gene expression profile of potato tubers under DON stress at different times, we performed RNA-seq on the same samples as those used in metabolomics. Fifteen samples from five groups of treatments were sequenced with libraries, and each li-

Transcriptional Analysis of Potato Tubers under DON Stress
To study the gene expression profile of potato tubers under DON stress at different times, we performed RNA-seq on the same samples as those used in metabolomics. Fifteen samples from five groups of treatments were sequenced with libraries, and each library of Clean Data was greater than 6 GB, Q30 base percentage was greater than 90%, average length was 150 bp, and GC content of each sample was greater than 43% (see Table S1 for raw sequencing data composition statistics and Table S2 for quality control details), indicating that the transcriptome sequencing results were reliable and suitable for further analysis. We mapped each sample (at least 47 million clean reads) using hisat2 software (v2.2.1), in which more than 95.85% of the reads were compared to the reference genome, with a Unique mapped ratio between 51.78% and 52.81% (Table S3), and a Unique mapped ratio of approximately 50% may be due to the reference genome species being tetraploid resulted. Unique-mapped assemblies were quantified using the StringTie software (v2. 10.5) to determine the count of each gene, which was reported as transcripts per kilobase per million (TPM, transcripts per kilobase of transcripts per million mapped reads). After normalization of expression, a total of 98,701 genes in at least one treatment with TPM > 0 were considered to be expressing genes, and were further analyzed ( Figure S2A).
To assess the reproducibility of the RNA-seq data, we calculated Pearson correlation coefficients between the three biological replicates for each sample. The results showed that the Pearson correlation coefficients between the three biological replicates ranged from 0.88 to 0.99, with large differences between groups, and PCA graph also showed good sample consistency ( Figure S2B,C), indicating that the transcriptomic data had a high reproducibility.
We identified DEGs between the control (CK0h) and other treatments using |log2 (fold change)| ≥ 1 and Padj < 0.05 as screening criteria. Compared with the control, 1053, 20,943, 20,347, and 21,172 were identified at 0 h, 4 h, 12 h, and 48 h under DON stress, respectively ( Figure 4A). We found that most DEGs were identified under DON48h treatment, followed by DON4h, DON12h, and DON0h, compared with CK0h. The fewest DEGs were identified at DON0h, and the number of upregulated genes was twice as many as the number of down-regulated genes. The other three treatments (DON48h, DON12h, and DON0h) had approximately 20 times more DEGs than at DON0h. The number of upregulated DEGs and the number of downregulated DEGs were similar, indicating that DON induced significant changes in potato tuber gene expression after 4 h ( Figure 4B).
Analysis of DEGs at different treatment times showed that each treatment time had its specific DEGs, and the highest number of specific DEGs was found at DON4h, followed by DON48h treatment, whereas 97 and 30 consistently upregulated and downregulated, respectively, and expressed DEGs were found in the four sets of comparisons, indicating that these genes may play a crucial role in response to DON stress in potato ( Figure 4C,D). Expression pattern analysis showed that most DEGs under DON stress showed specific expression over time.
Gene Ontology (GO) and KEGG enrichment analyses of DEGs were performed. All protein sequences were annotated based on the EGGNOG website, and the Orgdb package was constructed to functionally annotate and enrich each group of DEGs according to GO classification. In the GO analysis of the CK0h-DON0h group, 420, 328, and 391 DEGs were classified as biological processes, molecular functions, and cellular components, respectively. In the GO analysis of the CK0h-DON12h group, 7854, 6783, and 7846 DEGs were classified into biological processes, molecular functions, and cellular components, respectively. In the GO analysis of the CK0h-DON48h group, 8013, 6928, and 7991 DEGs were classified as biological processes, molecular functions, and cellular components, respectively. Compared with the CK0h control, DEGs under DON stress from 0 to 48 h were significantly enriched in the secondary metabolite biosynthetic process (GO:0044550), phenylpropanoid biosynthetic process (GO:0009699), and phenylpropanoid metabolic process (GO:0009698) ( Figure 5, Tables S7-S10).
Analysis of DEGs at different treatment times showed that each treatment time had its specific DEGs, and the highest number of specific DEGs was found at DON4h, followed by DON48h treatment, whereas 97 and 30 consistently upregulated and downregulated, respectively, and expressed DEGs were found in the four sets of comparisons, indicating that these genes may play a crucial role in response to DON stress in potato ( Figure 4C,D). Expression pattern analysis showed that most DEGs under DON stress showed specific expression over time. Gene Ontology (GO) and KEGG enrichment analyses of DEGs were performed. All protein sequences were annotated based on the EGGNOG website, and the Orgdb package was constructed to functionally annotate and enrich each group of DEGs according to GO classification. In the GO analysis of the CK0h-DON0h group, 420, 328, and 391 DEGs were classified as biological processes, molecular functions, and cellular components, respectively. In the GO analysis of the CK0h-DON12h group, 7854, 6783, and 7846 DEGs were classified into biological processes, molecular functions, and cellular components, respectively. In the GO analysis of the CK0h-DON48h group, 8013, 6928, and 7991 DEGs were classified as biological processes, molecular functions, and cellular components, respectively. Compared with the CK0h control, DEGs under DON stress from 0 to 48 h were significantly enriched in the secondary metabolite biosynthetic process (GO:0044550), phenylpropanoid biosynthetic process (GO:0009699), and phenylpropanoid metabolic process (GO:0009698) ( Figure 5, Tables S7-S10). DEGs were annotated in the KEGG database to understand their biological functions and gene interactions. In the CK0h-DON0h group, 423 of the 1053 DEGs were enriched in 27 pathways, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and phenylalanine metabolism ( Figure S3A). In the CK0h-DON4h group, 5495 of the 20,943 DEGs were enriched in 79 pathways, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and photosynthesis antenna proteins ( Figure S3B). In the CK0h-DON12h group, 5928 of the 20,943 DEGs were enriched in 91 pathways, including photosynthesis- DEGs were annotated in the KEGG database to understand their biological functions and gene interactions. In the CK0h-DON0h group, 423 of the 1053 DEGs were enriched in 27 pathways, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and phenylalanine metabolism ( Figure S3A). In the CK0h-DON4h group, 5495 of the 20,943 DEGs were enriched in 79 pathways, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and photosynthesis antenna proteins ( Figure S3B). In the CK0h-DON12h group, 5928 of the 20,943 DEGs were enriched in 91 pathways, including photosynthesis-antenna proteins, phenylpropanoid biosynthesis, and photosynthesis ( Figure S3C). In the CK0h-DON48h group, 6036 of the 21,172 DEGs were enriched in 100 pathways, including phenylpropanoid biosynthesis, photosynthesis antenna proteins, and flavonoid biosynthesis ( Figure S3D). Notably, phenylpropanoid and flavonoid biosynthesis consistently occupied an important position in the enrichment of DEGs in KEGG during 0 to 48 h of DON stress, which is consistent with GO and metabolomic analyses. Therefore, phenylpropanoid and flavonoid biosynthesis may be the main pathways in response to DON stress.

Weighted Gene Co-Expression Network Analysis (WGCNA)
Genes with very low expression (average TPM < 1) were removed from this analysis to avoid the inclusion of spurious edges in the networks. We screened 52,387 highly expressed DEGs from 63,883 DEGs using co-expression network analysis. The power values were first filtered to ensure that the gene distribution conformed to a scale-free network. Hierarchical clustering of all the samples showed no outlier samples. When the power value (β value) was 20, several genes with similar expression profiles belonging to the same subnetwork were clustered into the same co-expression modules (Figures 6A and S4). Finally, 14 modules were generated, with module sizes ranging from 60 to 13,318 genes, of which 1461 genes were not clustered into co-expression modules ( Figure 6B). In addition, 400 genes were selected for visualization to determine the correlation within the module genes ( Figure S5).
We investigated the correlation between differential metabolites and 14 co-expression modules and found that three modules were significantly correlated with 15 substances (|r| > 0.80, p < 0.01): darkseagreen1 (6 positive correlations), darkseagreen3 (3 positive correlations), and darkgoldenrod3 (1 positive correlation), and seven modules showed a clear trend of positive expression. Among all modules, the darkseagreen1 module was significantly correlated with seven differential metabolites, suggesting that genes in this module are highly correlated with the resistance mechanism of potato plants under DON stress. KEGG enrichment of three highly correlated modules revealed that the phenylpropanoid biosynthesis, oxidative phosphorylation, cysteine and methionine metabolism, necroptosis, and phenylalanine metabolism pathways were significantly enriched ( Figure 6C), suggesting that these pathways may play an important role in the stress response of potato to DON.
A scale-free network is only linked to a few nodes whose connectivity can represent the importance of the nodes. Therefore, we used the MCC algorithm in Cytoscape 3.9.1 software cytohubba plug-in to identify all hub genes from the three highly correlated modules separately and took the top 20 genes as hub genes ( Figure 6D-F, Table S4). addition, 400 genes were selected for visualization to determine the correlation within the module genes ( Figure S5). We investigated the correlation between differential metabolites and 14 co-expression modules and found that three modules were significantly correlated with 15

Validation of Single-Gene Expression
To further verify the accuracy of the data analysis, we randomly selected 15 of these 60 hub genes for qRT-PCR validation. The results showed that the expression of each gene obtained by qRT-PCR analysis correlated with the TPM values obtained from the RNA-seq data, confirming the reliability of the RNA-seq data (Figure 7). lism, necroptosis, and phenylalanine metabolism pathways were significantly enriched ( Figure 6C), suggesting that these pathways may play an important role in the stress response of potato to DON.
A scale-free network is only linked to a few nodes whose connectivity can represent the importance of the nodes. Therefore, we used the MCC algorithm in Cytoscape 3.9.1 software cytohubba plug-in to identify all hub genes from the three highly correlated modules separately and took the top 20 genes as hub genes ( Figure 6D-F, Table S4).

Validation of Single-Gene Expression
To further verify the accuracy of the data analysis, we randomly selected 15 of these 60 hub genes for qRT-PCR validation. The results showed that the expression of each gene obtained by qRT-PCR analysis correlated with the TPM values obtained from the RNAseq data, confirming the reliability of the RNA-seq data (Figure 7).

Discussion
Late blight caused by P. infestans occurs in nearly all potato-growing regions and is a major devastating disease for potato production of which the extent of damage is determined by the local climate of the year [44]. Prevention is made more difficult by the rapid mutation of the pathogenic fungus to adapt to the environment in which it survives, by

Discussion
Late blight caused by P. infestans occurs in nearly all potato-growing regions and is a major devastating disease for potato production of which the extent of damage is determined by the local climate of the year [44]. Prevention is made more difficult by the rapid mutation of the pathogenic fungus to adapt to the environment in which it survives, by the development of resistance, or by the resulting weakened resistance of resistant potato varieties [45]. Plant activation of PTI by conserved PAMP molecules recognized by PRRs and induction of ETI by pathogen effectors recognized by resistance proteins (R proteins) can enhance plant resistance and adaptation to biotic or abiotic stresses, which involve genes and pathways associated with different mechanisms [46]. Recently, a high-quality genome of the Atlantic variety of potato has been reported, which may facilitate our understanding of the role of stress-induced resistance in potato [47]. In our laboratory, we preliminarily found that low concentrations of DON could induce resistance to P. infestans in potatoes. To understand the metabolic changes and regulatory mechanisms of active compounds during DON induction, we constructed transcriptomic and metabolomic profiles of potato tubers at different time points under low concentration treatment combined with broad target metabolomic profiles to investigate DEGs and differentially accumulated metabolites to explore the mechanisms associated with DON-induced resistance in potatoes.
Studies have shown that DON and nivalenol (NIV) convert toxins to their glucosides (DON3G and NIV3G) via xenobiotic detoxification [48]. DON treatment induces ROS production in wheat. DON treatment rapidly induced the transcription of some defense genes in a concentration-dependent manner [49,50]. Additionally, it causes H 2 O 2 production in the grain, which triggers antimicrobial defense but may also stimulate cell death. Cloning DON resistance genes or transgenic detoxification-related genes can improve wheat blast resistance [51]. In 5 ng/mL DON treatment, we found that DEGs at all four time points were enriched in plant stress-related pathways, such as phenylpropane metabolism and synthesis and secondary metabolite synthesis, suggesting that the mechanism of potato defense induced by DON stress is mainly in the phenylpropane pathway. Multiple KEGG pathways were significantly enriched, with phenylpropanoid biosynthesis being the most abundant pathway, followed by flavonoid biosynthesis. These changes in the transcriptome correlated with metabolomic data. However, DON is still more toxic, and the number of virulence factors is a key element in inducing disease resistance in plants using excitons. Some studies have shown that both 3-keto-DON and 3-epi-DON show lower toxicity in converting DON to less toxic derivatives, and exhibit similar gene expression profiles to DON under 3-keto-DON treatment, which provides crucial information for future transgenic plant applications and development [52].
JA and SA are the major signaling pathways in plants, and these two signaling pathways are usually antagonistic; however, synergistic or additive effects have also been reported [53,54]. In this study, key genes associated with JA transduction, such as JAR1 and JAZ, were upregulated, while COI1 and MYC2 were significantly downregulated, and key genes associated with SA transduction, such as NPR1, TGA, and PR-1, were significantly upregulated ( Figure S6). Substance accumulation was also reflected in the metabolic data. As reported in existing studies, DON induction can enhance plant defense by regulating the gene expression of SA and JA pathways, causing SA accumulation and activation of pathogenesis-related genes (PR) [55]. Therefore, we hypothesize that DON induces plant defense through the interaction of SA and JA signaling, and that this process is regulated by DNA methylation [43,56].
Many changes occur in the primary metabolisms of plants, such as amino acids, nucleotides, and their derivatives, in response to stress, and the composition and levels of amino acids affect the resistance of plants [57,58]. Many plant defense compounds are derived from amino acid precursors, such as secondary metabolites and thioglucosides [59]. In the present study, 31 amino acids and their derivatives were significantly altered in potato tubers following DON induction. Among them, upregulation was mostly seen in toxic amino acids or sulfur-containing amino acids, such as L-tyramine, S-(5 -adenosyl)-Lhomocysteine, S-adenosyl-L-methionine, and L-alanyl-L-phenylalanine, which may play a key role in DON-induced plant defense. Consistent with this finding, genes related to these amino acids and their derivatives were also induced.
Plant secondary metabolites such as alkaloids, phenolic compounds, and flavonoids are typical plant antitoxins that play an important role in protection against diseases [60]. However, the content and distribution of secondary metabolites in different species differ greatly [61]. In our metabolomic analysis, the metabolites varied at different times, with changes in phenolic acids, flavonoids, terpenoids, lignans, and coumarins. Notably, phenolic acids, alkaloids, and flavonoids accounted for the highest percentage of the differential metabolites in each comparison group. Similarly, Isayenka et al. [62] found that the phytotoxin thaxtomin A (TA) induced browning of potato tuber flesh due to the accumulation of phenolic compounds, which may play a role in protecting potato tubers from Streptomyces scabiei. In addition, it is possible that because the family of tyrosinederived metabolites is flavonoids, DON induction caused significant downregulation of L-tyrosine, resulting in a reduction of 84% of all flavonoids in DAMs, suggesting that flavonoid levels may be negatively correlated with DON-induced plant defense. Combined with the transcriptome enrichment results, the phenylpropane synthesis pathway is a rich source of plant metabolites and a starting point for the production of many other important compounds such as alkaloids, coumarins, and lignans [31]. Therefore, the accumulation of gene expression and metabolites associated with phenylpropanoid biosynthesis suggests that the defense mechanisms after DON induction may be related to phenylpropanoid metabolism. This is consistent with the findings of Yogendra et al. [63,64] that metabolic pathways associated with the biosynthesis of phenylpropanoids, flavonoids, and alkaloids were strongly induced.
WGCNA divided all DEGs into 14 modules, and three modules were highly correlated with 15 shared DAMs. KEGG enrichment analysis showed that genes of related modules were mainly enriched in phenanthrene biosynthesis, which is consistent with the above results. The annotation of hub genes was mainly for functional protein-encoding genes such as ATTPS6 and COMT [65]. Most hub genes in the three modules were related to the regulation of plant immune responses. For example, COMT in the darkseagreen3 module is a key enzyme in phenylpropanoid biosynthesis, catalyzing a variety of phenylpropanoids and flavonoids, including caffeic acid, caffeol, and 5-hydroxyferulic acid, and catalyzes the production of melatonin from N-acetyl-5-hydroxytryptamine [66]. The ASMT encoded N-acetylserotonin methyltransferase in the darkseagreen1 module is the final enzyme in the biosynthetic pathway that produces melatonin. ASMT is induced in rice by both ABA and methyl jasmonic acid treatment [67,68]. Studies have shown that melatonin administration ameliorates the toxic effects of DON in mice [69]. MYB48 encoded transcription factor plays an important role in regulating stress response in response to JA signaling [70], and a previous study showed that MYB48 plays a central role in improving tomato tolerance [71] and a positive response of MYB48 to stress treatment in rice [72,73]. The cytochrome P450 family is not only involved in various metabolic pathways, such as alkaloids, terpenoids, and phenylpropanoids to form resistance substances in plants and improve plant resistance, but also acts as an antitoxin detoxification enzyme in plants [74]. CYP71B36 and CYP76C2 belong to the cytochrome P450 family, of which CYP71B has been shown to be involved in terpene biosynthesis (nerol and geraniol). For example, CYP71BJ1 in rose, which forms part of the pathway leading to 19-O-acetyl horhammericine, will help elucidate how this branch point is controlled [75]. CYP76C is a sesquiterpene plant antitoxin hydroxylase that hydroxylates rishitin in potato to form rishitin-M1 [76]. ABCG22 encodes an ABC transporter protein and the ABC transporter protein TaABCC3.1 has been reported to contribute to DON tolerance [77]. qRT-PCR assays showed altered expression of these hub genes at different times under DON treatment, suggesting that these hub genes may play an important role in DON-induced plant immune responses. This study has important implications for the design of defense strategies against P. infestans in potatoes.

Plant Materials and Treatments
In this study, Atlantic variety potato tubers were used as experimental material (provided by Heilongjiang Academy of Agricultural Sciences, Harbin, China). To study the dynamic response of the whole potato genome and metabolome to DON stress, potato tubers were treated with 5 ng/mL Vomitoxin (DON, purchased from FERMENTEK, Jerusalem Israel) and sampled at 0 h, 4 h, 12 h, and 48 h (recorded as DON0h, DON4h, DON12h, DON48h), and 0 h under sterile water treatment as control (recorded as CK0h), and all samples were frozen in liquid nitrogen and stored at −80 • C. Three biological replicates were performed for all five groups of samples.

Metabolite Profiling Using UPLC-MS/MS
The preserved samples were freeze-dried under vacuum using a lyophilizer (Scientz-100F, SCIENTZ Biotechnology Co., Ningbo, China) and the freeze-dried samples were crushed to powder form using a grinder (MM 400, Retsch, Verder Shanghai Instruments and Equipment Co., Ltd., Shanghai, China) at 30 Hz for 90 s. 100 mg of sample powder was weighed and dissolved in 1.2 mL of 70% methanol extract at 4 • C, and vortexed every 30 min for 30 s for a total of 6 times, and the samples were placed in a refrigerator at 4 • C overnight. After centrifugation at 12,000 rpm for 10 min, the supernatant was aspirated, and the sample was filtered through a microporous membrane (0.22 µm pore size) and stored in the injection vial for subsequent UPLC-MS/MS analysis.
For each sample, three biological replicates were analyzed independently. Chromatography mass spectrometry acquisition was performed using ultra performance liquid chromatography (UPLC, SHIMADZU Nexera X2, Kyoto, Japan) and tandem mass spectrometry (MS/MS, Applied Biosystems 4500 QTRAP, Foster, CA, USA).LIT and triple quadrupole (QQQ) scans were performed by a triple quadrupole linear ion trap mass spectrometer (Q TRAP), AB4500 Q TRAP UPLC/MS/MS system. The chromatographic column was an Agilent SB-C18 (1.8 µm, 2.1 mm × 100 mm); mobile phase was ultra-pure water in phase A (with 0.1% formic acid): acetonitrile in phase B (with 0.1% formic acid); and elution gradient: 5% in phase B at 0.00 min, linearly increasing to 95% in phase B at 9.00 min and maintaining at 95% for 1 min (10.00-11.10 min). The B-phase ratio decreased to 5% and equilibrated at 5% until 14 min, with a flow rate 0.35 mL/min, a column temperature of 40 • C, and injection volume 4 µL. For the mass spectrometry analysis [78], the ESI source operating parameters were as follows: ion source, turbo spray; source temperature, 550 • C; ion spray voltage (IS), 5500 V (positive ion mode)/−4500 V (negative ion mode); ion source gas I (GSI), gas II (GSII), and curtain gas (CUR) were set to 50, 60 and 25.0 psi, respectively, and the collision-induced ionization parameters were set to high. The instrument was tuned and mass calibrated in QQQ and LIT modes with 10 and 100 µmol/L polypropylene glycol solutions, respectively. QQQ scans were performed using MRM mode with the collision gas (nitrogen) set to medium. Further DP and CE optimization was accomplished for individual MRM ion pairs, and a specific set of MRM ion pairs was monitored in each period based on the metabolites eluted within each period.
Metabolites were characterized and quantified based on Metware Biotechnology Ltd.'s own database MWDB and public metabolite databases (MassBank, KNAPSAcK, HMDB, MoToDB, ChemBank, PubChem, NIST Chemistry Webbook and METLIN) [79]. The MS data were characterized by a precise comparison of precursor ion (Q1), product ion (Q3) values, and retention time (RT). The metabolites were quantified by MRM analysis using QQQ mass spectrometry. The characteristic ions of each substance were screened by triple quadruple rods, and the signal intensities (CPS) of the characteristic ions were obtained in the detector. The integration and calibration of the peaks were performed on the sample off-board mass spectrometry files using MultiaQuant software (v3.0.3), the peak area (Area) of each peak represented the relative content of the corresponding substance, and the integrated data of the area of all peaks were obtained [80].

Metabolome Data Analysis
Principal component analysis for unsupervised pattern recognition was performed on the filtered data using the prcomp function (parameter scale = TRUE) in R software (v4.2.1) [81]; orthogonal partial least squares discriminant analysis (OPLS-DA) was performed by the OPLSR.Anal function in MetaboAnalyst 4.0 software, based on the OPLS-DA results, combined with the difference fold change values; and the VIP values of OPLS-DA model were used to screen out the differential metabolites based on the OPLS-DA results, and the screening criteria were: metabolites with fold change ≥2 and fold change ≤0.5 and VIP ≥ [82]. Metabolites with significant differences were normalized by unit variance scaling (UV), and heatmaps were drawn by the pheatmap package. The differential metabolites were annotated and enriched using the KEGG database [83].

RNA-Seq Analysis
Total RNA was extracted from each sample, and RNA-seq libraries were constructed using Illumina Stranded Total RNA Prep (three biological replicates per treatment) based on the Illumina Novoseq 6000 system, and the fragment size and concentration of the libraries were detected by Agilent 2100 Bioanalyzer. QC of the RNA-seq raw data was performed using FastQC (v0.11.9, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 27 October 2022) and Trimmomatic (v0.36) to obtain clean reads [84], which were mapped to the Atlantic potato reference genome [47] by Hisat2 (v2.2.1) [85]. Expression level quantification was performed using FeatureCounts (v2.10.5) and expression normalization was performed with edgeR software (v3.36.1) using the TPM (transcripts per kilobase per million mapped reads) algorithm [86,87]. Differential gene identification was performed using the default parameters of the DESeq2 (v1.38.0) package in R (v4.2.1), and DEG identification was based on a cutoff value of |log2(fold change)| ≥ 1 as a screening criterion [88]. GO annotation was performed through the EGNOG online website (eggNOG-mapper (embl.de) for functional annotation of all protein sequences in the Atlantic potato genome [89], the AnnotationForge package built the protein annotation result files into the OrgDB package, and the self-built package was used to perform GO ORA enrichment analysis of DEGs. The enricher package was used to perform KEGG enrichment analysis on the DEGs. All results were visualized using R.

Weighted Gene Co-Expression Network Analysis (WGCNA)
Gene co-expression networks were constructed using R software version 4.2.1 with WGCNA (v1.71) [90]. Lowly expressed genes (TPM < 1) in DEGs were removed and coexpression networks were constructed by WGCNA analysis. Based on Pearson correlation coefficient calculation, the signed function was transformed to construct the correlation matrix between genes (threshold 0.8). The correlation matrix was transformed into an adjacency matrix by the power adjacency function. To characterize the nonlinear relationship between genes, a topological overlap matrix (TOM matrix) was constructed from TOM correlation coefficients, through R software package blockwise Modules function, to build the coexpression network by dividing modules and merging similar modules to the module feature vector (MEs matrix) and character matrix in order to analyze the correlation between the module feature genes and DON processing different time/control processing, which would estimate the correlation of the module character [91]. The hub genes were calculated by Cytoscape software (v3.7.2) and network visualization was completed [92].