Deciphering of benzothiadiazole (BTH)-induced response of tomato (Solanum lycopersicum L.) and its effect on early response to virus infection through the multi-omics approach

One of the preventive methods used to limit the losses caused by viruses is the application of synthetic immunity inducers, such as benzo(1,2,3)-thiadiazole-7-carbothioic acid S-methyl ester (BTH). This study aimed at explaining how the BTH treatment affects the defence and developmental processes in tomato plants (Solanum lycopersicum L.) as well as plant response to virus infection. The comparative multi-omics analyses concerning tomato plants treated with BTH were performed, including transcriptomics (RNA-Seq), proteomics (Liquid Chromatography-Mass Spectrometry), and metabolomics (targeted hormonal analysis). To confirm the priming effect of BTH on tomato resistance, the plants were infected with tomato mosaic virus (ToMV) seven days post-BTH treatment. The combined functional analysis indicated the high impact of BTH on the plant's developmental processes and activation of the immune response early after the treatment. In the presented experimental model, the increased level of WRKY TRANSCRIPTION FACTORS, ARGONAUTE 2A, thiamine and glutathione metabolism, cell wall reorganization, and detoxification processes, as well as accumulation of three phytohormones: abscisic acid, jasmonic acid-isoleucine (JA-Ile), and indole-3-carboxylic acid (I3CA), were observed upon BTH application. The immune response activated by BTH was related to increased expression of genes associated with the cellular detoxification process, systemic acquired resistance, and induced systemic resistance as well as post-transcriptional gene silencing. Increased levels of I3CA and JA-Ile might explain the BTH's effectiveness in the induction of the plant defence against a broad spectrum of pathogens. For the first time, the BTH involvement in the induction of the thiamine metabolism was revealed in tomatoes.


Introduction
Tomato (Solanum lycopersicum L.) is one of the most widespread vegetable crops, essential in the human diet, and a rich source of minerals, antioxidants, carotenoids (mainly β-carotene) and vitamins (Yusufe et al. 2017). Exposed to a wide array of stressors, including abiotic and biotic factors, the tomato plants have developed various mechanisms of response to ubiquitous stimuli (Santner and Estelle 2009;Verma et al. 2016). Plant hormones (phytohormones) are responsible for maintaining metabolic balance in plants, starting from growth, through developmental processes, and ending with defence response. Crucial phytohormones associated with this equilibrium are, particularly, salicylic acid (SA), jasmonic acid (JA), abscisic acid (ABA), and auxin (Aux/IAA) (Denancé et al. 2013). The whole-plant signalling pathways are interconnected. The sophisticated interaction between hormones and genes related to each signalling pathway enables plants to survive unfavourable conditions.
The metabolic homeostasis of plants can be compromised by various environmental factors, including pathogens, among other viruses (Suzuki et al. 2014). Importantly, besides prevention, there is no other method to control virus spread. To this aim, a combination of chemical and biological methods supported by agricultural technologies is being used (Maksimov et al. 2020). The best approach to minimize the risk of crop losses caused by viruses is, inter alia, the use of the resistant plant varieties in the production, the eradication of diseased plants from the field, and control of the population of virus vectors, among others insects (Nicaise 2014). An important group of prevention agents reducing the negative effects of virus infection are natural or synthetic immunity inducers, such as chitosan (Chirkov 2002) or thiadiazoleacetamides, for instance, benzo(1,2,3)-thiadiazole-7-carbothioic acid S-methyl ester (BTH) (Zhao et al. 2006) or its derivatives. Resistance inducers mobilize the plant immune system, and most importantly, are not phytotoxic (Pospieszny 2017). However, activation of the plant defence by resistance inducers is often associated with a fitness cost (Trejo-Saavedra et al. 2013).
Pathogens trigger a wide range of defence mechanisms in the plant, for instance, changes in the cell wall composition, activation of oxidative burst, synthesis of resistance compounds, or up-regulation of defence-related genes, among others those encoding pathogenesis-related proteins (PRs) (Dodds and Rathjen 2010). The BTH, a functional analogue of SA, activates the same pathways as in the pathogen response and initiates the induction of systemic acquired resistance (SAR) in the treated plants (Smith-Becker et al. 2003). The first use of BTH in modern agronomy dates back to 1996 (Thomson et al. 1998). Recently, research on BTH has been intensified and focused on its mechanisms of action and effectiveness in various biological systems. Noteworthy, new chemically modified BTH derivatives are being synthesized to improve their effectiveness through better assimilation by plants, higher dissolution rate, biodegradability, and protective potential (Smiglak et al. 2016(Smiglak et al. , 2017. The BTH-mediated resistance induction was evaluated, among others, in pepper against pepper golden mosaic virus (Trejo-Saavedra et al. 2013) and rice against sheath blight disease (Sood et al. 2013), indicating that the developmental and maturation processes of treated plants have not been disturbed by the BTH exposure at the optimal dose and concentration level in different host plants. Also, the global transcriptomic analyses of wheat (Li et al. 2020), strawberry (Landi et al. 2017), and banana (Cheng et al. 2018) plants under BTH treatment indicated upregulated genes associated with immune response. There are also reports on the positive effect of BTH treatment of tomato plants against Botrytis cinerea (Achuo et al. 2004) and activation of basal response to tomato spotted wilt virus and citrus exocortis viroid (López-Gresa et al. 2016. So far, the performed research has tended to focus on BTH as a SAR activator, when a positive relationship between marker genes of SA and JA (both increased after inducer treatment) was observed possibly suggesting the synergistic mode of action of these phytohormones in defence induction (Frąckowiak et al. 2019). Although all the mentioned studies provide general information about the decrease in the accumulation levels of those pathogens, the detailed mode of BTH action in tomato plants is still poorly understood.
Previously, we showed that the BTH and its ionic derivatives enhanced the plant immune response in Nicotiana tabacum-tomato mosaic virus (ToMV) pathosystem. Plants activated general immune response a few hours after BTH application and have shown up-regulated expression of SAR and ISR marker genes, including PR-1b and PR-2, LIPOXY-GENASE, ETHYLENE RECEPTOR (Frąckowiak et al. 2019). Moreover, in tobacco plants treated with the BTH either before or after ToMV inoculation, the accumulation of viral RNA in systemic leaves was significantly reduced, indicating the priming effect of tested inducers on the analysed host plants (Frąckowiak et al. 2019). We have also observed that the priming effects in BTH-treated tobacco plants were maintained up to a month after its application to soil.
This study aimed to comprehensively reveal the BTH mode of action in tomato plants, which has not been analysed before. Additionally, to determine the strength of immune priming to a viral pathogen, BTH-treated and non-treated plants were infected with ToMV. To decipher the overall BTH-induced response of the plant, the data obtained from transcriptomic and proteomic analyses were combined and enriched with the results of metabolomic analyses. We found out that the response of tomato plants to BTH treatment was related to changes in the expression level of WRKY transcription factors, ARGO-NAUTE 2A (AGO2A), genes involved in thiamine and glutathione metabolism (with GLUTHATHIONE S-TRANSFERASE (GST) gene expression), cell wall reorganization, detoxification processes and differences in the amount of three phytohormones: ABA, JA-Ile, and Aux.

Tomato plants, chemical treatment with selected inducers and virus inoculation
The experimental model consisted of a host plant (S. lycopersicum cv. Betalux), a virus (SL-1 isolate of ToMV), and a resistance inducer-the benzo(1,2,3)thiadiazole-7-carbothioic acid S-methyl ester (BTH). The plants were maintained in greenhouse conditions with a 16-h day/ 8-h night cycle in 26°C day/23°C night.
BTH was synthesized at the Poznań Science and Technology Park, AMU Foundation (Poland). It was used for watering six-to eight-week-old plants (100 mL) in the concentration of 10 mg/L. Initially, a group of twelve plants was treated with water and another twelve plants were treated with BTH. Next, water-treated and inducer-treated plants were collected 1-day post BTH treatment (1dpt). After one week, 4 plants previously water-or BTH-treated, were inoculated with the purified virus inoculum (at a concentration of 7.2 × 10 −4 ng per leaf) by rubbing carborundum-dusted leaves. The second sampling of the plant material was performed at 8 dpt. Four replicates were performed for each analysed experimental condition. All collected samples were frozen in liquid nitrogen and stored at -80°C. Portions of the collected plants from 1 and 8 dpt were freeze-dried using a laboratory freeze dryer-Alpha 1-4 L.S.C. basic (Martin Christ, Germany) and prepared for targeted hormonal analysis.
Isolation of total RNA and cDNA synthesis For the RNA sequencing (RNA-Seq), the total RNA was isolated using the acidic phenol-chloroform method (Molnár et al. 2005) with some modifications. The isolated RNA was suspended in 30-50 µL of RNase-free water. Next, RNA was DNase-treated (Thermo Scientific, USA), and the purified RNA (between 3-5 µg) was stored at -80°C until use. The quality and concentration of the prepared RNA were assessed using a 2100 Agilent Bioanalyzer (Agilent, USA).
For validation of the obtained results, the total RNA was isolated using the RNeasy Plant Mini Kit (Qiagen, Germany). The purified RNA (2 µg) was used for cDNA synthesis with the RevertAid RT Reverse Transcription Kit (Thermo Scientific, USA) in the presence of 100 µM random hexamers (Thermo Scientific, USA).

RNA-Seq and data analysis
The purified total RNA samples from all 24 plants were used for cDNA libraries synthesis. The libraries and RNA-Seq were prepared by MGI Tech Co., Ltd. (China) as described before (Zhu et al. 2018). The quality check of the obtained raw data reads and adapter trimming was performed by Genomed S.A. (Warsaw, Poland). The filtered clean reads from each sample were aligned to the reference tomato genome of S. lycopersicum SL3.0 from the International Tomato Genome Sequencing Project (ftp:// ftp. solge nomics. net/ tomato_ genome/ assem bly/ build_3. 00/) with annotation file ITAG3.2 (ftp:// ftp. solge nomics. net/ tomato_ genome/ annot ation/ ITAG3.2_ relea se/). Mapping and subsequent analyses were made using the OmicsBox software (OmicsBox-Bioinformatics made easy. BioBam Bioinformatics (Version 1.2.4). March 3, 2019. www. biobam. com/ omics box) (Dobin et al. 2013;Anders et al. 2015;Anders 2018). The count table obtained after the mapping process was filtered using the edgeR packet in the R environment (Version 3.6.2) (Robinson et al. 2009).
The differential gene expression analysis was carried out through a comparison between the treated and control groups from each time point of the experiment using OmicsBox (Robinson et al. 2009). The number of clean reads that mapped to each annotated transcript after RNA-Seq analysis was calculated using the HTSeq algorithm followed by the TMM (Trimmed Mean of M values) normalization method. The GLM (Quasi Likelihood-F test) statistical test was used in further data analysis. The DEGs were filtered by p-Value ≤ 0.05 and log2 fold change -1 ≥ (FC) ≥ 1.

Proteomic analysis
For protein isolation, 0.2-0.4 g of plant material was homogenized using liquid nitrogen and in the presence of the extraction buffer (0.1 M Tris-HCl, pH 6.8; 2 mM EDTA-Na 2 ; 20 mM DTT and 2 mM PSMF with 1% SDS) followed by centrifugation (15,000 × g, 3 min, 4˚C) to remove plant debris. Then 100% chilled acetone was added to the supernatant at the ratio (5:1) and left overnight at -20˚C. The precipitate was centrifuged (5,000 × g, 5 min, 4˚C) and the resultant pellet was rinsed three times with 80% chilled acetone. The precipitate was then dried using a Vacuum Concentrator (Heraeus Instruments, Denmark) and stored at -80 ˚C until used.
A detailed description of the subsequent stages of protein analysis is presented in (Supplementary File 1A). To prepare a list of statistically significant DEPs, the results were filtered by using p-Value < 0.05 and log2FC ≤ 0.75 (decrease in abundance) or log2FC ≥ 0.75 (increase in abundance) as the threshold.

Functional analysis of DEGs and DEPs
The annotation file of S. lycopersicum genes was downloaded from the BioMart Ensembl database (https:// www. ensem bl. org/ bioma rt/ martv iew/ e0215 cf00b ebb5d 16399 1ca78 59c5e 80). To obtain functional annotation and identify putative biological pathways of statistically significant DEGs and DEPs separately, the NCBI Gene Ontology (GO) (Harris et al. 2004) and the Kyoto Encyclopaedia of Genes and Genomes (KEGG) (Kanehisa and Goto 2000) databases were used, all performed in the Omics-Box software. For functional analysis of DEPs, the protein ID was mapped to the Ensembl Genome Transcript IDs using the UniProt mapping service. The p-Value < 0.05 was used as a threshold to define significantly enriched GO terms (Fisher's exact test) and KEGG pathways. To determine pathways enriched with transcriptomic and proteomic data, the combined file of DEGs and DEPs for each experimental condition was generated (including common terms). The OmicsBox software and R software with GOplot (Walter et al. 2015), ggplot2 (https:// www. rdocu menta tion. org/ packa ges/ ggplo t2/ versi ons/3. 3.3) and heatmap (https:// github. com/ raivo kolde/ pheat map) packages were used to specify the enriched GO terms and to visualize the results of GO and KEGG analysis. The STRING application in the Cytoscape software was used (Doncheva et al. 2019) to determine the most statistically significant KEGG pathways common for all experimental conditions. To perform this analysis, all of the statistically significant DEGs and DEPs have been converted to STRING ID using the g:Profiler browser tool (g:Converter) (https:// biit. cs. ut. ee/ gprofi ler/ conve rt). Results were collected and presented using PowerPoint software. The Venn diagram was prepared using an online tool Venny 2.1 (Oliveros 2007).
Targeted hormonal analysis-liquid chromatography coupled to ESI mass spectrometry The phytohormones were analysed as described by Sánchez-Bel et al. (2018)  . A detailed description of the subsequent stages of targeted hormonal analysis is presented in the Supplementary File 1B. The Student's t-test and Wilcoxon test were conducted in an R environment to examine significant differences between BTH-and watertreated samples. The p-Value was calculated for each experiment (with cut-off set at p < 0.05).
Gene expression studies were performed for four biological replicates, and for each biological replicate, three technical replicates were analysed. Relative gene expression was calculated using the GenEx software (MultiD Analyses AB, Sweden) and REST software (Pfaffl et al. 2002). An independent non-parametric Mann-Whitney test was conducted to test for significant differences between inducer (BTH)-and water-treated samples in all experiments. The p-Value was calculated for each experiment (with cut-off set at p < 0.05 (*)).

Results
General transcriptomic (RNA-Seq) and proteomic data analysis Plants (leaves) were harvested at two-time points, on day 1 and day 8 post the BTH treatment (dpt). Additionally, seven days after the BTH or water (control/ mock) treatment, the plants were inoculated with ToMV and their leaves were harvested one day after virus inoculation, to evaluate the impact of BTH on the virus-infected plant (8dpt/1dpi) (Fig. 1). Analyses presented below, including datasets obtained for BTH-treated and/or virus-infected plants were compared to datasets obtained for their corresponding control plants as follows: 1 dpt BTH vs 1 dpt control, 8 dpt BTH vs 8 dpt control, 8 dpt BTH / ToMV 1dpi vs 8 dpt control, and water/ToMV 1dpi vs 8 dpt control.
In general, 21,834 annotated transcripts for all experimental conditions were obtained (Supplementary Table 1) with the highest number of statistically significant differentially expressed genes (DEGs) identified in BTH-treated plants at 1 dpt (1,469 upregulated, 1,173 down-regulated) and the lowest number of DEGs identified for ToMV-infected plants (1 dpi) (140 up-and 170 down-regulated) compared to control plants at each time point (Fig. 1). The list of all DEGs is presented in Supplementary Table 2.
All the possible logical relationships between DEGs and DEPs in each experimental condition were presented on the Venn diagram (Fig. 1B  ToMV-infected plants at 1 dpi, no common DEGs and DEPs were identified, which may be due to the fact that shortly after virus infection, changes in the proteome corresponding to changes in the transcriptome are not yet evident. Main biological processes and KEGG pathways involved in tomato response to BTH To evaluate the impact of BTH treatment on tomato plants, the functional analysis of transcriptomic and proteomic datasets was performed. The Gene Ontology (GO) terms were reduced to the most specific ones, followed by the selection of the 15 most enriched biological processes for each type of experimental condition (Supplementary Table 3). The KEGG pathway analysis results were also reduced to the 20 most enriched in plants at the early (1dpt) and late (8 dpt) phase of response to BTH. This analysis was also performed for ToMV-infected and noninfected plants (Supplementary Table 4).
At the early stage of tomato response to BTH treatment (1 dpt), the majority of up-regulated processes at the transcriptional level were associated with the detoxification, photosynthesis, glutathione metabolic process, defence induction, ribosomal assembly process, and auxin (Aux)-mediated pathway related to root development. The regulation of the root organization process was also presented in the up-regulation of DEGs related to root development and proteins connected to lateral root formation and auxin polar transport. Additionally, 5 up-regulated DEPs connected to the response to the ABA process were noted. Both ABA and Aux may play a key role at the first stage of tomato response to BTH application. Interestingly, majority of the transcripts associated with the protein phosphorylation process and positive regulation of transcription by RNA polymerase II was down-regulated (71 and 18 of tested transcripts, respectively) (Supplementary Table 3).
In contrast to the early response, in BTH-treated tomato plants at 8 dpt, the up-regulations of DEGs and DEPs associated with the protein phosphorylation and down-regulation of the photosynthesis were observed. Moreover, we noted changes in the cell wall reorganization process, namely up-regulation of cell wall macromolecule catabolic process (transcriptomic data) and chitin catabolic process (transcriptomic and proteomic results), together with down-regulation of cell wall organization and biogenesis processes, xyloglucan metabolic process, and cellulose catabolic process (transcriptomic datasets). Likewise as in the early stage, in the late phase of tomato response to BTH application the up-regulation of glutathione metabolic process, detoxification process, and defence response to fungus and bacterium were observed at transcriptomic and proteomic levels. Furthermore, at the transcriptome level, the up-regulation of SAR and plant-type hypersensitive response (HR) were noted (Supplementary Table 3).
After viral infection in mock-and BTH-pre-treated plants, we detected an increased number of DEGs and DEPs associated with the response to oxidative burst (response to oxidative stress, hydrogen peroxide catabolic process, cellular oxidant detoxification) and a reduced number of transcripts and proteins linked to photosynthesis process (light-harvesting in photosystem I, electron transport in photosystem II, photorespiration, and chlorophyll biosynthetic process) in contrast to results observed for BTH-treated tomato plants 8 dpt. Noteworthy, in the water-treated-ToMV-infected tomato plants we observed up-regulation of the ABA catabolic process (transcriptomic data) together with down-regulation of thiamine (vitamin B1) biosynthetic process (proteomic data), which play a role in plant developmental process and defence induction in plant response to ToMV infection (Supplementary Table 3).
The results described above were additionally confirmed by KEGG pathway analysis. In all experimental conditions, the commonly represented pathways were related with starch and sucrose metabolism (mostly down-regulated at transcriptome level of BTH-treated plants both for 1 dpt and 8 dpt, and up-regulated in both transcriptome and proteome levels of BTH-ToMV-treated plants 1 dpi) and purine metabolism (highly up-regulated in BTH-treated plants 1 dpt and BTH-treated-ToMV-infected plants 8dpt/1dpi). Interestingly, the thiamine metabolism process was greatly up-regulated at transcriptome level in BTH-treated plants 1 dpt, at proteome level in BTH-treated plants 8 dpt, and at both transcriptome and proteome levels in BTH-ToMV-treated plants 1 dpi. At the early response of tomato plants to BTH treatment and after virus infection, the up-regulation of drug xenobiotic biodegradation and metabolism were noted, which was associated with detoxification and regulation of oxidative burst processes. Changes in the carbon fixation in photosynthetic organisms pathways were up-regulated in BTH-treated tomato plants 1 dpt and down-regulated in BTH-treated mock-inoculated (8 dpt) or ToMV-infected (1 dpi) plants (Supplementary Table 4).

Transcriptomic and proteomic datasets comparison
All DEGs and DEPs records were unified and the combined data were subjected to functional analysis. The plant processes and KEGG pathways, statistically most significant, were defined using the GOplot package in the R environment and Cytoscape software with String Enrichment application. Functional analysis of tested DEGs and DEPs showed enrichment of many processes and signalling pathways in BTHtreated plants, including primary metabolism, transcription and translation, cell wall organization, phosphorylation, photosynthesis, detoxification, response to phytohormone, and defence response (Fig. 2).
Similarly, as mentioned above, increased levels of the DEGs and DEPs associated with the photosynthesis, glutathione metabolic process, and gibberellin catabolic process were observed for BTH-treated plants 1 dpt (Fig. 2).
At the late phase of response, in BTH-treated plants 8 dpt, the reorganization of the cell wall was also noted (by up-regulation of xyloglucan metabolic process, cell wall macromolecule catabolic process, chitin-binding, and chitinase activity, and The size of bubbles is proportional to the number of DEGs and DEPs unified IDs assigned to each GO term. The x-axis represents the z-score and the y-axis represents the negative logarithm of the adjusted p-Value calculated in the R package. The threshold was set to 2.5 of -log(adj p-Value). Each category of GO term is represented by a different colour: Biological Process -red; Cellular Component -green and Molecular Function -blue down-regulation of lignin biosynthetic process and chitin catabolic process). Comparison of the transcriptomic and proteomic datasets obtained for BTHtreated plants 8 dpt showed that the protein phosphorylation and defence response processes were mostly up-regulated (Fig. 2).
In the BTH-treated-ToMV-infected plants 8 dpt/1 dpi, the up-regulation of processes concerning protein phosphorylation, cell wall reorganization process, and defence response were observed. Interestingly, the photosynthetic electron transport processes in photosystems I and II were up-regulated, while the photosystem reaction centre I, photosystem II oxygen-evolving complex, and chloroplast envelope components were highly down-regulated (Fig. 2).
The KEGG analysis showed significant changes in pathways associated with thiamine metabolism (up-regulation in BTH-treated plants 1 and 8 dpt) and up-regulation in phenylpropanoid biosynthesis and phenylalanine metabolism in all experimental variants of BTH-treated plants. In healthy BTH-treated plants (1 and 8 dpt) a slight increase in stilbenoid, diarylheptanoid, and gingerol biosynthesis was observed. At the late response to BTH treatment of tomato plants, an increase in ubiquinone and another terpenoid quinone biosynthesis was indicated together with a decrease in glyoxylate and dicarboxylate metabolism (both, in healthy and ToMV-infected plants). There was also a difference between healthy BTH-treated plants 8 dpt and BTH-ToMV-treated plants in purine metabolism, namely this process increased in healthy plants, while decreased after virus infection (Fig. 3). Drug metabolism, glutathione metabolism, and flavonoid biosynthesis were noticeably up-regulated in BTH-treated plants 1 dpt (Fig. 3). What is striking, all processes presented in Fig. 3 for water-treated-ToMV-infected plants were down-regulated (Fig. 3).
The significant up-regulation of pathways associated with detoxification and oxidative stress regulation, defence induction (plant-pathogen interaction), and hormone signalling regulation were also noted in the BTH-treated tomato plants. Also, the ribosome pathway was significantly up-regulated in the Fig. 3 Statistically significant KEGG pathways, previously analysed in the OmicsBox software. The visualisation was generated in the R environment (GO plot R package). Each dot represents one DEG/DEP (red dot -up-regulated or blue dot -down-regulated). The outer radius of the circle contains the pathway ID (explained in each table), while the internal radius represents the statistical significance (trapezoid size) and the direction of regulation of the presented pathway (based on the Z-score) BTH-treated plants at 1 dpt (46 of all associated DEGs/DEPs) and protein kinase activity process and phenylpropanoid biosynthesis pathways (with phenylalanine metabolism pathway) were mostly up-regulated in the BTH-treated plants 8 dpt (Fig. 4).

Effect of BTH treatment on changes in the levels of phytohormones in tomato plants
To determine the impact of BTH on the levels of phytohormones, seven different hormones were selected: ABA, SA, three associated with the JA pathway: JA-Ile, JA, and cis-( +)-12-oxo-phytodienoic acid (OPDA), and two connected with the auxin pathway: indole-3-acetic acid (IAA) and I3CA (Table 1).
In BTH-treated plants 1 dpt, the SA and the OPDA accumulation levels were lower than in control plants (1 dpt) 1.32 times and 1.89 times, respectively. In BTH-treated, non-infected plants (8 dpt) statistically significant changes were observed in the levels of phytohormone ABA (1.10 times higher than in control plants), and an increase was observed in I3CA (1.55 times) relative to control plants. Also, a high increase in the level of JA-Ile (2.73 times) was noted in BTH-treated tomato plants at 8 dpt.
The ToMV infection affected the abundance of the tested phytohormones in both BTH-treated and watertreated tomatoes. In the infected plants, in comparison to the control plants, BTH caused slight changes in the amount of ABA (1.12 times lower), SA (1.26 times lower), and I3CA (1.27 times higher). Interestingly, the level of JA-Ile was 3.59 times higher than in control plants (p-Value = 0.072). In virus-infected plants, treated only with water, an increase in JA-Ile (3.35 times), JA (2.39 times), and I3CA (1.66 times), and a decrease in SA (1.15 times), relative to controls, was observed.

BTH mode of action in tomato
To describe BTH's possible mode of action in tomatoes, genes involved in rapid response (BTH-treated Fig. 4 Maps of selected pathways, enzymes, DEGs/DEPs associated with tomato plants' response to BTH treatment (1 dpt and 8 dpt) and/or virus-infected (treated with BTH [BTH-ToMV] or with water [ToMV]) received from Cytoscape application. The number of found genes associated with each term is presented in brackets. Nodes have been grouped by experimental condition, where the green arrows show the number of up-regulated genes and the red arrows show the number of down-regulated genes. The bottom part presents heatmaps generated from the logFC values of each process/enzyme/protein presented above from each experimental condition. Each colour bar represents each term that is shown above and indicates the range of genes observed for a given term on heatmaps Table 1 Results of differences in the amounts of synthesized phytohormones acquired from the analysed plants calculated on a dry weight The first value represents mean results from at least 3 plants from each experimental condition in ng/g DW (dry weight). The second value (after ±, in italics) is recalculated SD (standard deviation) from each result. ABA -Abscisic Acid; SA -Salicylic Acid; JA-Ile -Jasmonic Acid-Isoleucine; JA -Jasmonic Acid; OPDA-cis-( +)-12-oxo-phytodienoic acid; IAA-Indole-3-acetic acid; I3CA-Indole-3-carboxylic Acid; ND -not detected. Statistical analysis, with calculated fold change, were presented in the bottom panel of the  Table 5). During the first phase of plant response to the BTH treatment, significant changes in various processes were noted, concerning, among others, protein translation (of all 77 DEGs and 11 DEPs the up-regulation was observed for 73 DEGs and 8 DEPs), defence response (of all 60 DEGs and 8 DEPs -up-regulation of 43 DEGs and 7 DEPs), and photosynthesis (of all 43 DEGs and 3 DEPs --up-regulation of 42 DEGs and 1 DEPs) (Fig. 5 and Supplementary Table 5). The highest number of DEGs were identified for BTHtreated plants 1 dpt, and among them were the upregulated genes encoding transcription factors, particularly 8 genes belonging to the WRKY family (e.g. Solyc02g080890 with 11.1-fold change (relatively to untreated plants]), and Solyc08g062490 with 4.8-fold change), 12 genes belonging to ethylene-responsive transcription factor family (e.g. Solyc01g090320 with 122.4-fold change, and Solyc04g071770 with 25.1-fold change), and 2 genes belonging to bHLH transcription factor family (mostly down-regulated, e.g. Solyc03g119390 with 24.3-fold change, Solyc12g036470 with 9.3-fold change). Among the processes related to plant defence response in BTH-treated plants 1 dpt, a high expression level was observed for genes belonging to the F-box protein family (e.g. Solyc01g106520 with 56.0-fold change, Solyc01g106130 with 49.9-fold change, and Solyc01g106140 with 23.7-fold change), while a change in the expression of pathogenesis-related genes (e.g., PR1, PR2, PR3) was generally observed later (at 8 dpt). The changes in photosynthesis were associated with the early phase of plant response to BTH and manifested with the up-regulation of genes related to photosystem I (e.g. Solyc12g033000 with 30.6-fold change, Solyc12g033040 with 20.6-fold

Fig. 5 A The GO results (Biological Process -red, Cellular
Component -green and Molecular Function -blue) from the Omicsbox software present DEGs and DEPs identified for GO records significally associated with tomato plants response to BTH application and/or ToMV infection. An asterisk at the GO name indicates which processes have been enriched (p-Value < 0.05). The orange bars specify the number of up-regulated genes, and the blue bars specify the number of downregulated genes; B) Selected Biological Processes in which the expression of genes has significantly changed. (All data from each experimental condition were obtained from analysis of each treated plant compared to the corresponding "Control" plant) change, and Solyc02g020960 with 9.7-fold change), and photosystem II (e.g. Solyc05g041230 with 41.96fold change, and Solyc09g016930 with 35.4-fold change) (Supplementary Table 5). Considering genes associated with the response to phytohormones pathways (64 of all DEGs and 11 of all DEPs) at 1 dpt, we can divide them into those associated with response to Aux (17 DEGs/4 DEPs up-and 11 DEGs/ 1 DEP down-regulated records), ABA (3 DEGs/ 5 DEPs upand 9 DEGs/ 1 DEP down-regulated), JA (5 DEGs/ 1 DEP up-and 6 DEGs down-regulated), and ethylene (4 DEGs up-and 5 DEGs down-regulated). Noteworthy, up-regulation of the genes belonging to the GST family, associated with signalling, defence, and hormone response processes, was also observed in plants 1 dpt (e.g. Solyc09g011500 with 123.1-fold change, Solyc09g011590 with 5.7-fold change or 2.4fold change at proteome level, and Solyc07g056470 with 28.5-fold change) ( Fig. 5 and Supplementary  Table 5).
At 8 dpt in BTH-treated plants we observed down-regulation of the photosynthesis (from all 10 DEGs and 11 DEPs -down-regulation was observed for 9 DEGs and 8 DEPs), up-regulation of the phosphorylation (from all 146 DEGs and 37 DEPs -up-regulation of 113 DEGs and 19 DEPs), and dephosphorylation process (14 DEGs/ 8 DEPs upregulated). Like in the case of BTH-treated tomato plants at 1dpt, the enhanced expression of WRKY family genes (e.g. Solyc04g072070 with 11.5-fold change, and Solyc05g015850 with 10.6-fold change) and the ethylene-responsive transcription factor family (e.g. Solyc11g011750 with 15.8-fold change, and Solyc09g066350 with 4.9-fold change) were observed in BTH-treated plant at 8 dpt. At the late stage of tomato response to BTH treatment, an increase in the up-regulation of the defence response (from all 64 DEGs and 17 DEPs -the up-regulation of 56 DEGs and 12 DEPs was observed), including SAR response (up-regulation of 7 DEGs and down-regulation of 1 DEG) was observed. The enhanced expression of pathogenesis-related protein family genes have been mainly identified in tomato plants 8 dpt after BTH application, including PR1 (Solyc09g007010 with 29.1 fold change), BETA-1,3-GLUCANASE (PR2, Solyc01g097240 with 10.4-fold change), CHITINASE (PR3; Solyc10g074400 with 25.9fold change), and OSMOTIN-LIKE PROTEIN (PR5; Solyc11g044390 with 24.2-fold change). High levels of expression have also been found for the genes belonging to the NIM protein family (NIMIN2c PRO-TEIN; Solyc03g119590 with 15.7-fold change and NIM1-INTERACTING PROTEIN; Solyc05g009480 with 7.9-fold change), CYSTEINE PROTEASE (e.g. Solyc02g076910 with 69.9-fold change), and MAJOR ALLERGEN D1 PROTEIN (Solyc09g091000 with 14.7-fold change) in BTH-treated 8 dpt tomato plants ( Fig. 5 and Supplementary Table 5).
After viral infection of BTH-pre-treated plants 8 dpt/ 1 dpi, a detailed analysis showed changes in the regulation of photosynthesis (of 121 analysed DEGs and 5 DEPs all of them were down-regulated) and phosphorylation processes (of all 214 DEGs and 20 DEPs, -the up-regulation was observed for 141 DEGs and 15 DEPs). Additionally, the expression level of PR4 was increased (Solyc01g097240 with a 10.1-fold change), which was not observed in noninfected plants (Fig. 5 and Supplementary Table 5).
The expression level of the gene encoding methyl esterase (e.g. Solyc03g070380), one of the major components of SA synthesis and SAR induction, changed in time after BTH application to tomato plants. In the early stage of the tomato response to BTH treatment, its expression was down-regulated (24.6-fold change), while in the late stage it was upregulated (8.8-fold change in healthy plants and 6.2fold change in ToMV-infected plants) (Supplementary Table 5). The contribution of the peroxidase genes in the detoxification process is also high, both in the initial and late response processes of tomato plants to BTH (e.g. up-regulation of Solyc01g067870 with 20.2-fold change and Solyc11g018777 with 16.9-fold change in BTH-treated plants 1 dpt, Solyc02g090450 with 6.5-fold change and Solyc04g071890 with 6.2-fold change in the healthy BTH-treated plants 8 dpt, Solyc11g018775 with 19.9-fold change and Solyc11g018777 with 7.6-fold change in the BTHtreated-ToMV-infected plants (at 8 dpt/1dpi)) (

Early response of tomato to ToMV infection after BTH treatment
To identify DEGs, DEPs, biological processes, and pathways associated with tomato response to ToMV infection after BTH treatment additional analyses were performed, including a comparison of datasets obtained for BTH-treated-ToMV-infected 8 dpt/ 1 dpi to datasets obtained for water-treated-ToMV-infected 1 dpi (as a control) and also to the datasets obtained for BTH-treated plants 8 dpt (the second control) ( Fig. 6 and Supplementary Table 6).
The analyses showed 1,424 DEGs and 271 DEPs for BTH-treated-ToMV-infected plants at 8 dpt/ 1 dpi compared to water-treated-ToMV-infected plants at 1 dpi, and 554 DEGs and 412 DEPs, when datasets of BTH-treated-ToMV-infected plants at 8 dpt/ 1 dpi were compared to the datasets obtained for BTH-treated plants at 8 dpt (Fig. 6A). The GO analysis of enriched processes revealed up-regulation of defence response to fungus protein phosphorylation, galactose catabolic process via UDP-galactose and glutathione metabolic, as well as down-regulation of negative regulation of endopeptidase activity, response to wounding, protein hexamerization, response to light stimulus and adaxial/ abaxial pattern specification. Worth mentioning is that virus infection in BTH-treated plants is associated with the up-regulation of negative regulation of reductive pentose-phosphate cycle process. Also, down-regulation of photosynthesis and chlorophyll biosynthetic processes were observed (Fig. 6B).
Pathways induced by virus infection in BTH pre-treated plants were related to up-regulation of unified IDs assigned to each GO term. The x-axis represents the z-score and the y-axis represents the negative logarithm of the adjusted p-Value calculated in the R package. The threshold was set to 2.5 of -log(adj p-Value). The outer radius of the circle in KEGG analysis contains the pathway ID (explained in each table), while the internal radius represents the statistical significance (trapezoid size) and the direction of regulation of the presented pathway (based on the Z-score). The visualisation was generated in the R environment (GO plot R package) thiamine metabolism, phenylpropanoid biosynthesis, primary metabolism, and detoxification (drug metabolism -cytochrome P450 and metabolism of xenobiotics by cytochrome P450 pathways). The early response of tomato plants to viral infection was also associated with the up-regulation of alpha-linolenic acid metabolism (related to the JA pathway) and down-regulation of carbon fixation in photosynthetic organisms (related to photosynthesis) (Fig. 6B). The Venn diagram of DEGs and DEPs showed some common genes and proteins induced by BTH after ToMV infection with the same direction of expression profile (44 records, 32 up-regulated and 12 down-regulated) ( Table 2) and with the opposite direction of gene expression profile (14 records) (Table 3). Interesting in this case is the expression profile of Solyc09g083440 (PIN-I protein) and Solyc09g089505 (Proteinase inhibitor I), which in both additional analyses presented the down-regulation on transcriptome level, and up-regulation on proteome level ( Table 2).
The highest number of statistically significant results was identified for BTH-treated plants 1 dpt, opposite to BTH-treated and ToMV infected tomato plants (Fig. 7B) (Fig. 7B).

Discussion
Changes in differentially expressed genes (DEGs) and proteins (DEPs) upon BTH treatment BTH is known as a synthetic plant protection agent that mimics SA properties and enhances the plant immune system before the appearance of a pathogen, which has been confirmed for many plant species, among others, rice (Sood et al. 2013), strawberry (Landi et al. 2017), banana (Cheng et al. 2018) and potatoes (Brouwer et al. 2020). In the present study, BTH-induced resistance of tomato plants from in-depth transcriptomic, proteomic, and targeted metabolomics analyses were described. For the first time, early and later phases of S. lycopersicum response to BTH were analysed in detail, including the effect on the initial stage of viral infection. BTH application to tomato plants resulted in a high increase in the DEGs number at 1 dpt, while in BTH-treated plants 8 dpt fewer DEGs have been identified. The contrary state was observed for DEPs, where an increase of a higher number of DEPs was noticed at 8 dpt rather than 1 dpt.

BTH treatment causes up-regulation of WRKY transcription factors genes in tomato plants
Transcription factors (TFs) play one of the most important roles in the cellular machinery, including WRKY family transcription factors, which were identified amongst the most up-regulated TFs in plants upon BTH treatment. Li et al. (2020) have highlighted the induction of WRKY family gene expression in the plants after BTH treatment, postulating their indirect effect on immunity induction through the synthesis of genes from the PR family in an NPR1-independent way (Li et al. 2004(Li et al. , 2020. In our study, changes in the expression levels of WRKY genes were observed, of which 3 were related to the developmental process (WRKY6/35/75), while 5 associated with the plant's defence response (WRKY 33/50/51/80/81) (Rushton et al. 2010). The upregulation of two genes belonging to the WRKY family, namely: WRKY6 (Solyc02g080890) and WRKY50 (Solyc08g062490), was found in BTH-treated tomato at each analysed time point (with the highest expression change at 1 dpt). Those two genes are supposed to be activated by BTH application and may have a crucial role in maintaining the metabolic balance of tomato Impact of BTH treatment on post-transcriptional gene silencing process in tomato The plant defence response is also related to the process of gene silencing, including PTGS (Vaucheret et al. 2001). In Arabidopsis thaliana pre-treated with BTH, the enhanced gene expression of the AGO2 gene was observed due to cucumber mosaic virus infection (CMV) (Ando et al. 2021). The AGO proteins, together with small RNAs, are involved in the formation of the RNA-induced silencing complex (RISC), the core element in the PTGS (Grene 2002). In the presented study, a significant increase in the expression levels of AGO2A and DICER-LIKE 2D (DCL2D) (Solyc11g008530) genes in BTH-treated plants at 1 dpt was detected. In late response to BTH treatment of tomato plants, including virus-infected ones, the expression level of AGO2A gene was decreased, while an increase in the AGO2A protein level was observed. These findings may suggest that AGO2 and DCL2D may be also involved in the defence response to ToMV infection in tomatoes.

Regulation of immune response in BTH-treated tomato plants
Under pathogen infection in Nicotiana benthamiana, the GST attenuated the oxidative stress inside chloroplasts, enabling the synthesis of virus RNA minus strand (Budziszewska and Obrępalska-Stęplowska 2018). GST is associated with the glutathione metabolic process and is responsible for signalling processes in the plant, which was also observed in the presented results (e.g. the increased expression of  Solyc09g011500 and Solyc07g056510 in all BTHtreated plants). However, a large number of GST transcripts were also associated with the response to hormones (mainly Aux) (e.g. mentioned before Solyc09g011500) and defence response (Supplementary Table 5). The protein kinases play essential roles in the cellular signalling processes and control the mechanisms of protein modification and activation in plant-pathogen interactions (Kersten et al. 2009). The analysis of tomato plants' response to BTH (1 dpt and 8 dpt) indicated a contribution of some of the identified DEGs and DEPs in the phosphorylation process and kinase activity function ( Fig. 5 and Supplementary Table 5). The increased expression of protein kinases genes (among others MAP KINASE KINASE 2, Solyc03g123800) and protein (MAP KINASE KINASE 1, Solyc12g009020) in virus-infected plants and their important role in the pathogenesis process have been described before (Wrzesińska et al. 2018(Wrzesińska et al. , 2021, also for BTH-treated plants (Cheng et al. 2018). Such significant changes in protein kinase levels probably contribute to maintaining the balance between the growth and the plant's defensive functions, which is supported by the lack of visible phytotoxicity effects on BTH-treated plants.
At the early stage of tomato response to BTH, the increased expression of defence protein gene precursor, signal transduction protein (GST), and genes associated with hormone response (like JAZ genes Solyc12g009220 and Solyc08g036660) was observed. Our previous research has revealed that BTH and its derivatives have a positive effect on defence induction not only related to SAR-type response but also cause changes in the JA pathway marker genes expression (like JAZ gene), which is associated with induced systemic resistance (ISR) (Frąckowiak et al. 2019). A likely explanation is that both types of plant defence responses (SAR and ISR) are activated during BTH application to tomato and tobacco plants. In the late phase of tomato response to BTH treatment, the increased expression of PR genes (including the PR2) was also reported in tobacco plants (Frąckowiak et al. 2019), and increased expression levels of genes associated with SAR and programmed cell death (such as Solyc06g071050, HYPERSENSITIVE-INDUCED RESPONSE PROTEIN) were also noted in this study.
The response of tomato plants to the BTH treatment included changes in expression levels of DEGs and DEPs associated with the organization of the cell wall and its modification in all experimental conditions. The up-regulation of cell wall macromolecule catabolic process in BTH-treated tomato plants 8 dpt (both in healthy and ToMV-infected plants) together with down-regulation of the lignin biosynthetic process in BTH-treated tomato plants 8 dpt were observed. Also decrease in the expression level of GLYCOSYLOTRANSFERASE (Solyc07g043480), which is involved in lignin biosynthesis pathways (Scheible and Pauly 2004), was confirmed for all experimental conditions of BTH-treated tomato plants. Together with the up-regulation of PR2 gene expression, these findings indicate that the plant's response is more complex and possibly explains why BTH can act on a wider spectrum of pathogens, viruses (Frąckowiak et al. 2019), bacteria (Brisset et al. 2000), fungi (Araujo et al. 2015), and some pests (Inbar et al. 1998)).

Changes in photosynthesis and vitamin B1 metabolism upon BTH treatment
Initial results suggest that there may be a relationship between BTH application to tomato plants and changes in the photosynthesis process. Photosynthesis leads to the production of oxygen, energy (ATP), and carbohydrates, and the biosynthesis of many primary compounds, vitamins, phytohormones (SA, JA, ABA), and defence proteins. Chloroplasts participate in the immune response through the ROS/ redox system (Grene 2002) and changes in cytosolic Ca 2+ concentration (Lu and Yao 2018), but they are also associated with virus replication (Budziszewska and Obrępalska-Stęplowska 2018). At the early phase of tomato response to BTH, the photosynthesis-associated processes are highly up-regulated, while in the BTH-treated plants 8 dpt (both in virusfree and virus-infected plants), the photosynthesis was down-regulated. It may indicate that this process is crucial in the cell machinery stabilization after BTH action, which has been also observed for strawberries and bananas (Landi et al. 2017;Cheng et al. 2018). The up-regulation of the thiamine (vitamin B1) metabolism pathway, bringing a highly significant contribution to the plant's defence response and development (Bocobza and Aharoni 2014;Dong et al. 2016), was observed in BTH-treated tomato plants 1 and 8 dpt. What is interesting, two genes involved in the thiamine metabolism pathway, showed the opposite direction of the expression level in BTH-treated plants. The first one, HEAT SHOCK PROTEIN 90 (Solyc04g081630), is associated with plant response to biotic and abiotic stresses (including defence response to fungus and response to heat stress) (Xu et al. 2012). This gene was up-regulated in BTH-treated 1 dpt and BTH-treated-ToMVinfected plants, while the second gene, THIAMINE THIAZOLE SYNTHASE (Solyc07g064160), associated with thiamine biosynthesis, was highly downregulated in all analysed BTH-treated plants. A former study confirmed the same situation when thiamine metabolism was increased with a reduction in its availability in plant stress conditions (Lee et al. 2007). Vitamin B1 acts as an enzymatic coenzyme in the maintenance of health and metabolism balance in plants (Fitzpatrick and Chapman 2020) and regulation of the response to abiotic (via ABA) and biotic stresses (Goyer 2010). It has been postulated that the increase in the expression of the genes associated with the thiamine metabolism pathway is correlated with the activity of SA (Ahn et al. 2005) and ABA (Rapala-Kozik et al. 2012) in rice, Arabidopsis, and cucumber. This is the first time that changes in the thiamine metabolism were found to correlate with BTH application.
BTH application has an impact on levels of I3CA and JA-Ile The phytohormones analysed in this study (ABA, SA, JA, and Aux) are associated, among others, with growth, signalling, transport, and immunity. Interestingly, the SA level was decreased at the early stage of tomato response to BTH application and in the plants infected with ToMV (both water-and BTH-treated). The decrease in the amount of SA in the BTH-treated plants both at 1 dpt and 8 dpt may be associated with the functional properties of the BTH molecule. The virus also has an impact on SA metabolism, both in water-and BTH-treated tomato plants, because a decrease in the amount of free SA in the plants was observed (higher in BTH-pre-treated tomato plants). Interestingly, an increase in the genes expression associated with the JA metabolic pathway in the BTH-treated plants (Solyc12g009220 and Solyc08g036660) was detected and the amount of JA-Ile was also increased, with statistically significant results for the BTH-treated 8 dpt plants. JA-Ile is a JA derivative synthesized in the cytoplasm and together with the increase in JAZ genes expression, it regulates the plant immune response and developmental processes such as growth, seed germination, and root formation (Ghorbel et al. 2021). This observation confirms the results previously reported for N. tabacum -another plant from the Solanaceae family (Frąckowiak et al. 2019), in which the expression levels of the genes related to the JA pathway also increased after the treatment with BTH or its derivatives. It seems that the application of BTH to plants influences a synergistic relation between SA-and JA-mediated pathways making the plant's defence response probably more effective. In the late response to BTH (8 dpt), the amount of ABA and I3CA increased in tested plants. For the BTH-treated-ToMV-infected tomato plants (8 dpt/1 dpi) an increase in the amounts of IAA and I3CA hormones was noted, while the level of ABA was decreased. Additionally, an up-regulation of gene expression associated with Aux/IAA (Solyc10g018340) and I3CA (Solyc02g064830), with the highest expression for BTH-treated plants 1 dpt was shown. IAA is important for plant root development and growth; it also plays a key role in the induction of AUXIN RESPONSE FACTOR PROTEIN gene expression and triggers the cellular immune response (Bari and Jones 2009). It has been reported that an increase in the amount of I3CA also induces ABA production (Gamir et al. 2018), which corresponds to callose deposition due to BABA treatment (Gamir et al. 2012). The data reported here appear to support the assumption that BTH influences the JA and I3CA biosynthesis pathways, but in general, it requires more studies to confirm that hypothesis.

Early response to virus infection in tomato plants treated with BTH
After virus infection in BTH pre-treated tomato plants, the early response was associated with upregulation of defense, detoxification, and phosphorylation processes, and with down-regulation of photosynthesis. Together with up-regulation of thiamine metabolism, phenylpropanoid biosynthesis, and alpha-linolenic acid metabolism pathways, our observations confirmed the previously obtained results and indicated a significant contribution of the JA and vitamin B1 pathways in response to ToMV infection. Interestingly, the synthesis of two proteins was up-regulated, when their transcript levels were down-regulated, namely Solyc09g083440 (PIN-I protein) and Solyc09g089505 (proteinase inhibitor I). The first one is associated with Aux pathway (Křeček et al. 2009) while the second one is associated with the JA pathway and was also identified as an important gene in tobacco response to ToMV in our previous work (Frąckowiak et al. 2019). These two proteins may play a key role in tomato response to ToMV infection at an early stage, however, it must be supported by further research.

Conclusion
Concluding, in this study, a multi-omics analysis was carried out for the first time to decipher the tomato response to BTH administration in tomato plants. This study indicated that tomato plants' response at the early stage after BTH application manifests itself through an increase in the expression levels of DEGs and DEPs belonging to the WRKY family and GST family with high activity of protein kinases. In the presented studies, we also indicated the BTH impact on the activation of genes expression associated with the PTGS process in tomato plants. In the BTH-treated tomato plants, the regulation of the cellular detoxification process, which protects cells from apoptosis and DNA damage caused by ROS was activated. Moreover, for the first time, the BTH application impact on the thiamine metabolism pathway was presented in tomato plants, which may suggest that this vitamin plays a role in plants' growth-defence balance in a tested experimental pathosystem. Also, changes in the amount level of I3CA and JA-Ile, together with up-regulation of the expression levels of these hormones marker genes (Solyc02g064830, Solyc12g009220, and Solyc08g036660) may suggest that BTH is not only responsible for SAR induction but also activates other pathways, thus increasing the effectiveness of the plant defence against a broad spectrum of pathogens.