Comparative Transcriptome Profiling of Two Tomato Genotypes in Response to Potassium-Deficiency Stress

Tomato is a crop that requires a sufficient supply of potassium (K) for optimal productivity and quality. K+-deficiency stress decreases tomato yield and quality. To further delve into the mechanism of the response to K+-deficiency and to screen out low-K+ tolerant genes in tomatoes, BGISEQ-500-based RNA sequencing was performed using two tomato genotypes (low-K+ tolerant JZ34 and low-K+ sensitive JZ18). We identified 1936 differentially expressed genes (DEGs) in JZ18 and JZ34 at 12 and 24 h after K+-deficiency treatment. According to the Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analyses, the DEGs that changed significantly primarily included transcription factors, transporters, kinases, oxidative stress proteins, and hormone signaling-and glycometabolism-related genes. The experimental results confirmed the induced expression of the responsive genes in the low-K+ signaling pathway. The largest group of DEGs comprised up to 110 oxidative stress-related genes. In total, 19 ethylene response factors (ERFs) demonstrated differential expression between JZ18 and JZ34 in response to K+-deficiency. Furthermore, we confirmed 20 DEGs closely related to K+-deficiency stress by quantitative RT-PCR (qRT-PCR), some of which affected the root configuration, these DEGs could be further studied for use as molecular targets to explore novel approaches, and to acquire more effective K acquisition efficiencies for tomatoes. A hypothesis involving possible cross-talk between phytohormone signaling cues and reactive oxygen species (ROS) leading to root growth in JZ34 is proposed. The results provide a comprehensive foundation for the molecular mechanisms involved in the response of tomatoes to low K+ stress.


Introduction
K + is required in most plants as a cationic mineral nutrient [1,2]. This element is essential for many physiological processes, including photosynthesis, enzyme activation, protein synthesis, osmoregulation, cell turgor, and ion homeostasis in plant cells [3]. Plant growth requires K + , but most plants can absorb only a minimal amount of soluble K + from the soil [4]. K + -deficiency can directly disturb physiological activities and restrict assimilation and partitioning into fruits from the nutrient source, resulting in decreased crop growth, and fruit production and yield [5]. Increasing the utilization of the K + nutrient or the resistance of the crop itself to K + deficiency is the major method that can solve the K + -deficiency problem. However, the molecular mechanisms responsible for these differences are not clear. The experiments were repeated three times. K + represents normal K + (4 mM); K − represents K + deficiency (0.5 mM). For each line, different lowercase and uppercase letters indicate significant differences (p < 0.05 and p < 0.01 respectively) among the treatments. K + can increase crop yields and improve crop quality. The K + concentrations in the K + -deprived seedlings of both genotypes declined after 7 d under K + -deficiency treatment. Nevertheless, JZ34 maintained a higher K + content under K + -deficiency stress, indicating that this cultivar was better at maintaining a high K + concentration under K + -deficiency than was JZ18. The K + content of JZ34 declined slightly under these conditions, with a relative coefficient of 0.89. However, the K + content in JZ18 dropped significantly, with a relative coefficient of 0.57. K + accumulation still trended similarly. To determine the K + uptake ability of both genotypes, K + was depleted in the uptake solution, which reflected its net uptake by the roots monitored over 30 h. External K + was depleted to a greater extent by JZ34 than by JZ18 (Figure 2). We determined the maximum ion absorption rate (I max ) and Michaelis-Menten constant (K m ) in the presence of 0.2 mM external K + . These two parameters quantitatively described the ability of the plants to absorb nutrient ions (JZ34: I max = 0.10 mmol·g −1 ·min −1 ; K m = 0.107 mmol/L; JZ18: I max = 0.06 mmol·g −1 ·min −1 ; K m = 0.125 mmol/L). The K + -uptake ability of JZ34 was better than that of JZ18 [20]. The K + uptake ability was compared using K + depletion between two different genotypes. The experiments were repeated three times. * and ** denote significant differences at p < 0.05 and p < 0.01, respectively.

Overview of the Gene Expression Profile Sequence Data and Mapping Results
To obtain a global overview of the transcriptome relevant to different K + -deficiency stress conditions in the tomato, we separately sequenced complementary DNA (cDNA) libraries from root samples of JZ34 and JZ18 at different time points under K + -deficiency (0.5 mM) using the BGISEQ-500 platform. After the removal of low-quality reads, an average of 2.4 × 10 7 clean reads was obtained for each sample, and the total length of the clean reads reached 1.2 × 10 9 nt. For each sample, 99% of the clean reads were mapped to the tomato reference transcriptome (Table 2). The expression levels of the mapped genes were normalized using FPKM values. To confirm the quality of the RNA-seq results, we evaluated the expression of the 12 highest-ranking housekeeping control genes in tomatoes at both 12 and 24 h, including GAPDH, catalase, Cys protease, α-tubulin, ubiquitin, actin, DNAJ, and translation initiation factor 5A ( Table 3). None of these reference genes were significantly differentially expressed in the JZ34 and JZ18 lines by pairwise comparison after the induction of K + deficiency (Table 3), indicating that the obtained sequences and transcript levels were suitable for further transcriptome analysis.  (Figure 3F-T). The gene function annotations are shown Table S1. The selection criteria for DEGs (F-T) were a difference in the expression levels between two genotypes that was greater than 5, and a number of DEGs in this group that was greater than 5. The resulting regulation model was consistent between the transcriptome and the qRT-PCR results (Figure 3), suggesting that our transcriptome analysis was reliable.  Table S1.Gene-specific primers used for real-time PCR are listed in Table S5.

Global Analysis of DEGs
To investigate the molecular responses of the tomatoes under the K + -deficient condition, we identified up-and downregulated genes at 12 and 24 h post-treatment using the t-test (false discovery rate (FDR) ≤ 0.001) and an absolute value of log 2 (fold change) ≥ 1 (treatment/control). We detected 1936 DEGs, including 966 upregulated and 970 downregulated genes. JZ34 had more upregulated genes than did JZ18, whereas JZ18 had more downregulated genes than did JZ34. The number of upregulated genes in JZ34 was almost twice as high as that in JZ18 (Figure 4). More DEGs were upregulated at 24 h than at 12 h. The gene expression patterns differed between JZ34 and JZ18.

DEGs Encoding Transporters and Kinases
Five genes encoding K + channels and ion transporters demonstrated significantly differential expression in this study, including three upregulated K + channel genes (Solyc11g011500, Solyc03g097930, and Solyc01g010480) and two downregulated K + transporter genes (Solyc12g009540 and Solyc07g014680). The expression of genes encoding nitrate and ammonium transporters changed greatly under K + -deficiency stress ( Table 4). As the largest subfamily of receptor-like kinases (RLKs), leucine-rich repeat receptor-like kinases (LRR-RLKs) regulate the growth, development, and stress responses of plants. A total of 24 LRR-RLK DEGs were found, with 19 downregulated genes and two upregulated genes detected in JZ18. None downregulated genes and eight upregulated genes detected in JZ34. The overall expression trends of the LRR-RLKs in the two genotypes were opposite. Table 4. Genes encoding protein transporters and kinases showed genotypic differences in response to K + -deficiency (0.5 mM) stress.

Group
Gene ID  The number of DEGs associated with oxidative stress metabolism was highest among all categories analyzed in the present study (Table S2). Among these genes, 42 peroxidase, 44 cytochrome P450, and 24 glutathione S-transferase genes were identified, including 46 upregulated and 56 downregulated genes. Approximately eight times as many genes were downregulated as upregulated in JZ18, whereas six times as many genes were upregulated as were downregulated in JZ34. The fluorescence signal was particularly strong in the guard cells and roots of JZ18 ( Figure 6). The determination of ROS concentrations in the plants revealed significant accumulation in the leaves and roots of JZ18 at 12 h and 24 h (Table 5). These results show that stronger ROS scavenging activity is induced by K + -deficiency stress in JZ34 than in JZ18.    DEGs involved in hormone signaling were analyzed (Table 5), including auxin (5), gibberellin (12), cytokinin (8), ethylene (3), and salicylic acid (2). Genes encoding abscisic acid (ABA) -related genes were not found among the DEGs. Most DEGs were downregulated in JZ18 and upregulated in JZ34. The number of auxin, gibberellin, and cytokinin genes was greater than five (Table 6). Plant hormones are very closely related to plant morphogenesis. To better explore the effects of phytohormones on tomato resistance under K + -deficiency and the effects on root architecture, the IAA and CK contents were determined and used to calculate the IAA/CK ratio under K + -deficiency (0.5 mM) treatment in the root systems of the two genotypes ( Figure 7). The IAA content in the roots of the two genotypes was higher after K + -deficiency treatment than in the control at 12 h. The content in JZ18 when the treatment was conducted for 24 h to 3 d was lower than that in the control group (Figure 7). At the same, the IAA content in JZ34 was higher than that in the control group. These results showed that the IAA contents in the roots of the two genotypes were differentially altered after K + -deficiency stress. The CK content in the roots of the two genotypes showed an almost opposite trend after K + -deficiency stress, with lower levels in JZ18 and higher levels in JZ34 compared with those in the control group, except at 12 h. Typically, the IAA/CK ratio is used to determine whether to promote bud growth or root growth. A lower ratio enables promotion of bud differentiation, whereas a higher ratio is more conducive for root formation.  The IAA/CK ratio was higher in the treatment group than in the control group at 12 h in both genotypes. With extended treatment, the ratio was lower in the treatment group than in the control group at 24 h, and the difference was significant in JZ18. At 3 d of treatment, JZ18 still showed a significantly lower ratio, whereas the JZ34 treatment group had a higher ratio than did the control group.
In this study, 14 DEGs were associated with glycometabolism and transport. These genes included members of the SWEET family, as well as a glucose transporter, sucrose phosphate synthase, and ATP-dependent 6-phosphofructokinase. Of these DEGs, 11 were upregulated, and two were downregulated in JZ34 ( Table 7). Most of the genes related to glycometabolism were detected in JZ34, and almost all were upregulated. Although a few genes were detected in JZ18, most of these genes were downregulated. In the transcriptome analysis, five DEGs directly related to root architecture were found. All these DEGs were all upregulated in JZ34 at 24 h (Table S3), with Solyc08g066590 was upregulated (2.31-fold) to the highest level compared to the expression of the other four genes.

GO and KEGG Analyses of K + -Deficiency Stress Tolerance-Related DEGs
Expression analysis by GO enrichment was conducted to identify DEGs that were significantly upregulated in JZ34, but downregulated or unchanged in JZ18, and DEGs that barely changed in JZ34, but were downregulated in JZ18. A total of 1210 DEGs met the above criteria under K + -deficiency stress in the two varieties ( Figure S2 and Table S4). A total of 551 DEGs was downregulated in JZ18 but unchanged in JZ34, whereas 517 DEGs were upregulated in JZ34 but unchanged in JZ18. Only 161 DEGs that were found in both genotypes were upregulated in JZ34 and downregulated in JZ18. More genes were differentially expressed at 24 h than at 12 h in the two accessions ( Figure S2 and Table S4). The DEGs were divided into 36 functional groups. GO functional enrichment analysis showed that genes associated with binding (GO: 0005488) and catalytic activity (GO: 0003824) were significantly enriched, accounting for as much as 86% of the molecular function domain (Figure 8). The GO terms "metabolic process", "cellular process", "response to stimulus", "single-organism process", and "biological regulation" accounted for the majority of biological processes (Figure 8). The top two dominant terms in the cellular components were cell and cell part. In total, 1210 DEGs encoding various enzymes were further matched with the enrichment of 117 KEGG pathways ( Figure S3). The pathways with the greatest numbers of unique sequences were all metabolic pathways, including those related to secondary metabolites, amino acids, nucleotides, lipids, carbohydrates, energy, and other metabolism. The significantly enriched KEGG pathways are shown in Figure S3. The comparative transcriptome analysis of JZ34 and JZ18 under K + -deficiency stress laid a foundation for further elucidation of gene functions and the metabolic pathways in tomatoes.

Discussion
K + -deficiency stress in soil is quite common and becomes more severe during crop production [21]. Plants have evolved different strategies to cope with K + -deficiency. In the present study, a transcriptome analysis was performed 12 h and 24 h after K + -deficiency stress in two tomato genotypes. The early sequencing of adversity stress is more conducive to a comprehensive analysis of the differences in resistance caused by changes in the molecular mechanisms within the two genotypes. At the same time, we observed apparent differences in the root systems of the two genotypes after K + -deficiency treatment for 7 d. The experimental results proved that plants start from sensory signals, transmit the signals for self-regulation, change their phenotype, and finally adapt or not adapt to stress.

K + Transporter Expression Results in Differences in Tolerance to K + Deficiency between the Genotypes
K + acquisition typically exhibits dual (high-and low-) affinity mechanisms in plants [9,22]. K + uptake permeases (KUPs) encode plasma membrane K + /H + symporters and catalyze K + influx into cells under low-apoplastic-K + conditions [23]. The transcript representing the KUP-related K + transporter is Solyc12g009540; its expression was downregulated (6.79-fold) in JZ18, implying that this gene might assist the KUP gene family in K + absorption under K + -deficient conditions. K + translocation between different organs and tissues inside a plant and K + secretion from root cortex cells into the xylem are mediated by outward-rectifying channels, such as SKOR [24]. Solyc03g097930 is a SKOR-like K + channel whose expression was upregulated in the low-K + -tolerant variety JZ34. Solyc01g010480 is a voltage dependent K + -channel that is expressed in root hairs, and it displayed upregulated expression in JZ34. KAT1 is a highly selective inward-rectifying voltage-gated K + channel that mediates long-term K + influx into guard cells, leading to stomatal openings [25]. Expression of the KAT1 gene transcript Solyc01g010480 was upregulated (1.60-fold) in JZ34, suggesting that this variety was better able to perceive K + deficiency, open stoma, and stimulate K + transport. In our research, these genes were verified by qRT-PCR (Table S1).

Root-Related Genes, the ROS Signaling Pathway, and ERFs Cause Changes in Root Growth under K + D Eficiency
Root morphological and physiological traits are important for the absorption of nutrients from the soil [26]. K + -efficient genotypes have a longer root system than do inefficient genotypes under K + -deficiency stress in rice [27], tomatoes [28], and Hordeum maritimum [29]. The differences in the root system of the two genotypes after K + -deficiency stress were very notable in our study. Similarly, K + accumulation was positively correlated with the root length and root surface area in cotton [30]. Ethylene and ROS were also essential for root hair elongation [31]. Ethylene signaling plays an important role in low-K + -induced plant responses, such as root hair elongation, ROS production, and HAK5 expression [32].
Our study showed that the expression of many ERFs was upregulated, including Solyc03g005520 (ERF1a) (8.16-fold) and Solyc03g005500 (ERF14) (8.37-fold) (Table S1), were significantly higher in JZ34 than in JZ18 under the same conditions. The root fresh/ dry weight of JZ34 was also higher than that of JZ18. Therefore, some ERFs may play an important role in root growth in tomatoes in response to low-K + signaling in tomato. Furthermore, ethylene is known to act upstream of ROS in response to K + deprivation in Arabidopsis [32]. K + -deprived plants rapidly accumulate ROS, and the expression of some K + transporters is thought to be dependent on ROS production [33].
Peroxidase, cytochrome P450, and glutathione S-transferase are critical for the regulation of ROS production, and for reducing cellular damage during oxidative stress in plants [34]. In the present study, the highest number of DEGs was associated with oxidative metabolism than any other category. Whether or not this effect is related to ethylene needs further validation, but the oxidative stress response is clearly a key physiological activity in the response to short-term K + -deficiency stress. Our results showed differences in the ROS scavenging capacities of the two genotypes under K + -deficiency stress conditions, which was most likely a key factor underlying the different tolerances of these genotypes to low-K + . Peroxidase activity and ROS signaling are specifically required during LR emergence [35]. GPX was a peroxidase-encoding gene. In transcriptome analysis, Solyc04g071890 (GPX1) and Solyc07g056480 (GPX7) were upregulated (2.07-fold) in JZ34 and downregulated (1.64-fold) in JZ18, respectively (Table S2). Overexpression of the two genes in Arabidopsis could enhanced root growth [36]. Solyc07g052370 (CYP) is a cytochrome P450 gene that was upregulated (6.13-fold) (Table S1) in JZ34 in the transcriptome analysis; this gene is involved in root defensive mechanisms in tomatoes [37]. Therefore, DEGs associated with oxidative stress are most likely key factors leading to differences in the root systems of the two genotypes under K + -deficiency stress.

K + Deficiency between the Genotypes Leads to Changes in Hormonal Responses and Biosynthesis-Related DEGs
Plant hormones are continuously involved in the regulation of plant morphological changes under adverse conditions. Auxins play an important role in the development of plants, LRs, adventitious roots and root hairs. The PIN-FORMED (PIN) family serves as auxin export carrier proteins [38]. Many PIN family members participate in the transport and distribution of auxins, and thus play an important role in the establishment of auxin distribution patterns. The present study identified five auxin-related DEGs, four of which belong to the PIN family. Solyc01g068410 (PIN5) was downregulated (2.65-fold) in JZ18 and up-regulated (2.83-fold) in JZ34. Studies have suggested that PIN5 is mainly located in the endoplasmic reticulum; its expression in root hair cells can stimulate root hair growth to a certain extent, and PIN5 enhances the availability of internal IAA [39]. Solyc10g080880 (PIN7) was downregulated in JZ18, and Solyc07g006900 (PIN2) was upregulated in JZ34. Both genes are involved in IAA transport, and they regulate root growth and development [40]. This observation shows that genes associated with IAA transport are more active in JZ34 than in JZ18 and that the IAA content is increased in JZ34 and is reduced in JZ18 after K + -deficiency stress. Therefore, under K + -deficiency stress conditions, the PIN family showed differential expression patterns in the two genotypes, leading to differences in the IAA content to a certain extent, which might be one cause of the root system differences between the two genotypes. In addition, ROS have the ability to decrease IAA signaling [41], and therefore ROS accumulation in the roots of JZ18 may inhibit IAA production and transport.
In our study, eight genes related to CK metabolism were identified. Solyc10g017990, (CKX3) was downregulated (7.65-fold) in JZ18 (Table S1). In Arabidopsis thaliana, CKX3 promotes the growth of the main roots and LRs, and has a higher root biomass than does the wild type. Similar effects have also been found for CKX1, CKX2, and a family member of CKX in tobacco [42]. Therefore, we believe that DEGs that are associated with CKs are probably key factors leading to differences in the root systems of the two genotypes under K + -deficiency stress. We also measured the CK content. The JZ18 CK content decreased after K + -deficiency stress, whereas that of JZ34 increased. Existing studies suggest that CK has a negative regulatory effect on root growth and development [43]. However, the physiological processes of plants are often achieved by the coordinated regulation of various hormones. CKs can regulate the expression of auxin and PINs, and the ratio of IAA and CK can be used to explore synergistic on plant configurations more reasonably. Studies have suggested that a higher ratio indicates a greater utilization of root development. In terms of the IAA/CK ratio, JZ34 gradually surpassed the control group at 3 d. At this time point, the root systems of the two genotypes began to change gradually in configuration (Figure 1). Notably, DEGs related to salicylic acid metabolism were also identified in the present study. Salicylic acid is a new type of plant hormone that regulates some important metabolic processes in plants [44]. Studies have shown that exogenous salicylic acid can induce increased growth in soybean roots [45] and act as a signaling molecule to induce sustained resistance to abiotic or abiotic stresses in various plants. Salicylic acid induces the production of many enzymes related to plant resistance, and regulates their activities [46]. After the initiation of K + -deficiency stress in the present study, salicylic acid related genes were downregulated in JZ18 and upregulated in JZ34, with significant differences observed between the two genotypes. This hormone may be one key factor affecting the different tolerances of the two genotypes.

Glycometabolism and Transport Involved in K + -D Eficiency Stress
Sugar is not only is the energy source for plants, but also functions in signal transduction, regulates the expression and enzyme activity of related genes, and participates in the regulation of physiological activities [47,48]. Some studies have suggested that soluble sugars are involved in the production of NADPH, such as oxidation of the pentose phosphate pathway (OPP), to promote the clearance of ROS [49,50] and can also be used directly to scavenge hydroxyl radicals [51]. In the present study, the low-K + -tolerant strain showed no obvious ROS accumulation during the K + -deficiency treatment. This effect may be due to not only its strong ability to clear ROS, but also differences in sugar metabolism. SWEET family, glucose transporter, and solute carrier family genes involved in transporting sugar from the source to sink were also detected among the DEGs [52]. These DEGs were upregulated in JZ34 and downregulated in JZ18. These genes promote the transport of sugar from the leaves to the roots, and provided a basis for root growth in JZ34; however, the sugar transport capacity was decreased, and root growth was inhibited in JZ18 under K + -deficiency stress. Of these DEGs, Solyc03g097560 (SWEET14) was most downregulated (6.82-fold) in JZ18 (Table S1), which include cis-acting elements that respond to auxin. Thus, this gene is very likely to have a negative impact on IAA production. Based on our analysis, sugar metabolism and transport-related genes appear to play a positive role in improving plant resistance to K + -deficiency stress. Genes directly related to root development were also screened by the transcriptome analysis, but no further studies were found for these genes in tomatoes. Among the 20 genes identified in this study, the genes discussed include ERF, PIN, CKX, CYP, and SWEET14, which we believe were more likely to have an effect on root configuration under low K + stress, and could provide a good basis for subsequent in-depth experiments.
Two tomato genotypes were exposed to short-term K + -deficiency stress in this work. Multiple physiological metabolic pathways are involved in responses to K + -deficiency stress (Figure 9). Combining the transcriptome analysis and physiological analysis results suggested that a potential interaction occurred between them. ERF and ROS metabolism-related enzymes are probably involved in ROS accumulation and clearance. Hormone signaling is instructive for changes in the IAA and CK contents. Glycometabolism may play a very positive role in changes in the root architecture. No investigation found that LRR-RLKs affected the root architecture under K + -deficiency stress conditions, but these results indicated that LRR-RLKs might be more important in the tomato response to K + -deficiency. The changes in the ROS, IAA, and CK contents at the physiological level explained the changes in root architecture, although the mechanism underlying K + -deficiency stress resistance requires further study. In addition, the mechanisms regulating K + absorption, transport, and utilization are important for tomato root system growth and plant yields with the increasingly widespread issue of K + -deficient soil. Our current study has provided valuable materials and data for further research into K + -deficiency stress.

Plant Materials and K + Deficiency
Seeds of two tomato genotypes (JZ34, low-K + -tolerant, and JZ18, low-K + -sensitive) were surface sterilized in 1.0 mL 10% (v/v) NaClO for 15 min and rinsed five times with distilled water. The seeds were incubated to allow them to germinate in plastic Petri dishes in the dark for 3 d at 25 • C on two layers of sterile filter paper that were previously soaked in water. Tomato seedlings were transferred to a tray containing a 1:1:1 (v/v/v) mix of peat: perlite: vermiculite and grown in a greenhouse from March to April with day/night temperatures of approximately 26/18 • C in a daily average irradiance of 350 µmol·m −2 ·s −1 and at 70% relative humidity (RH). Seedlings at the vegetative growing stage (25 d) were cleaned with water, washed three times with distilled water and transferred to a pot ( (30 d), K + deficiency treatment was applied by reducing the concentration of KNO 3 from 4 mM (normal K + ) to 0.5 mM (K + deficiency) in the nutrient solution. A concentration of 4 mM KNO 3 was used as a control. After 7 d of K + -deficiency stress, different parts of the plant were sampled and dried to assess plant length, dry weight, and K + content. K + accumulation = K + content × dry weight; K + use efficiency (KUE) = dry weight (g) × 1000/K + content (mg/g). To more clearly observe the apparent difference between the two genotypes after K + -deficiency stress, K + -free nutrient solution was also used. Differences in the leaf ( Figure S1) and root systems of the two genotypes were examined at 3 d and 7 d, respectively. The experiments were repeated at least three times.

Measurement of K + Uptake Kinetics
The nutrient solution of four-leaf-stage (25 day) tomato seedlings was replaced with a solution without K + for 3 d with day/night temperatures of approximately 26/18 • C. Three seedlings with consistent growth were washed three times with 0.2 mM CaSO 4 and placed into a black triangular flask with 100 mL of absorption liquid (0.2 mM KCl + 0.2 mM CaSO 4 ). This process was repeated three times with each strain. Next, 1 mL of absorption liquid was removed from the triangular flask at different time points, and 1 mL of distilled water was added; this process was repeated until the absorption liquid K + concentration was stable or close to zero. The K + uptake kinetic parameters (K m and I max ) were calculated according to the method of Drew et al. [53]. I max is the maximum uptake rate, and K m is the Michaelis-Menten constant.

RNA-seq Sampling and RNA Isolation
Seeds of JZ34 and JZ18 were germinated for RNA-seq sampling under the conditions described above and placed into a plant growth chamber. Four-leaf-stage seedlings at the vegetative growth stage (30 d) were exposed to K + -deficiency stress (0.5 mM) or the full-strength solution as control for 12 h and 24 h. Three whole independent seedlings were collected and mixed together at each time point to reduce the differences between individual plants. Three biological replicates were included for each treatment. The samples were snap-frozen in liquid nitrogen and stored at −80 • C for the RNA-seq analysis.
RNA was isolated according to the instructions of the RNeasy mini kit (Qiagen, Hilden, Germany). The RNA purity was determined using the NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA). The RNA concentration was measured using the Qubit ® RNA Assay Kit in the Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). The RNA integrity number (RIN) was assessed using the RNA Nano 6000 Assay Kit on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) with 2.2 ≥ A260/A280 ≥ 2.0, 28 S/JZ18 S ≥ 1.5 and RIN ≥ 7.5. The RNA abundances and purity were tested to confirm that they met our requirements.

Library Construction, Sequencing, and Data Processing
Library preparations were sequenced on the BGISEQ-500 platform (BGI, Shenzhen, China), and 50-bp single-end (SE) reads were generated. BGISEQ-500 is powered by combinatorial probe-anchor synthesis (cPAS) and improved DNA Nanoballs (DNBs) technology. The cPAS chemistry works by incorporating a fluorescent probe into a DNA anchor on the DNBs, followed by high-resolution digital imaging. This combination of linear amplification and DNB technology reduces the error rate while enhancing the signal. In addition, the size of the DNBs is controlled so that only one DNB is bound per active site. This patterned array technology not only provides sequencing accuracy but also increases chip utilization and sample density.

Identification of DEGs
High-quality clean reads were mapped to the reference transcriptome, and the transcript levels of the unigenes were identified by TopHat (http://ccb.jhu.edu/software/tophat/index.shtml) and Cufflinks (http://cole-trapnell-lab.github.io/cufflinks/) [54] and normalized by the Fragments Per Kilobase of exon model per million mapped reads (FPKM) approach [55] according to the following formula: FPKM = 10 6 C/NL/10 3 . Given the expression of gene A, C indicates the number of fragments aligned to gene A, N represents the total number of fragments aligned to all genes, and L denotes the number of bases in gene A.
Sequencing was used to identify genes that were differentially expressed in the two different varieties [56]. FDRs were used in multiple hypothesis testing to correct the p-values. The FDRs were statistically preset to values less than 0.05, and the differential gene expression in different samples was calculated based on the FPKM values [57]. The expression ratio of 12 h/0 h or 24 h/0 h is presented as the fold change in the present study. An FDR ≤ 0.001 and an absolute value of log 2 (fold-change) ≥ 1 were used as the thresholds to screen DEGs. The Multiexperiment Viewer (MeV; http://www.tm4.org/mev.html) [58] was used to delineate heatmaps based on the DEG results.

Gene Annotation, GO Enrichment and KEGG Analysis
The GO and pathway enrichment analyses of the DEGs were performed using the PANTHER classification system (http://www.pantherdb.org/data/) [59] and KEGG (http://www.kegg.jp/) [60]. A basic local alignment search tool (BLASTn) programme was used to identify homologues in the tomato genome (https://solgenomics.net/).

qRT-PCR Analysis
For quantitative qRT-PCR analysis, first-strand cDNA was synthesized using the PrimeScript 1st Strand cDNA Synthesis Kit (TIANGEN, Beijing, China). The qRT-PCR was performed in triplicate for each sample using the SYBR Green Real Master Mix, according to the manufacturer's instructions. qRT-PCR amplification was performed using the Bio-Rad CFX Manager 3.1 Real Time PCR System (Applied Biosystems, Foster City, CA, USA) and software (Applied Biosystems). The 2 −∆∆Ct method was used to analyze the relative changes in the gene expression levels from three biological replicates. The data were analyzed based on the mean values of triplicates. The specificity of the reactions was verified through melting-curve analysis. Actin was used as an internal control. The primers used for the qRT-PCR are listed in Table S5.

Determination of ROS
The cultivation method was the same as described above, and K + -deficiency treatment was used to detect ROS in tomato leaves and roots at 0, 12 h, 24 h, 3 d, and 7 d. Three biological replicates were conducted for each treatment. ROS accumulation was detected using the Reactive Oxygen Species Assay kit (Beyotime Biotechnology Shanghai, Shanghai, China) on the ZEISS Axio Observer A1 inverted fluorescence microscope. The H 2 O 2 and O 2 − concentrations in the leaves and roots were determined according to a modified method [61,62].

Endogenous Hormone Content Determination
At 12 h, 24 h, and 3 d of K + -deficiency treatment, 0.1 g of root material was collected from each genotype. The experiments were repeated at least three times. The hormone contents were analyzed using an enzyme-linked immunosorbent assay [63].

Statistical Analysis
The data were analyzed using one-way analysis of variance (ANOVA) followed by Tukey's honestly significant difference (HSD) multiple comparison tests in Statistical Productions and Service Solutions 17.0 (SPSS, Chicago, IL, USA). p < 0.05 was considered significant. Asterisks (* and **) indicate a significant difference between the controls and transgenic plants at p < 0.05 and 0.01, respectively. Author Contributions: X.Z. contributed to the experimental design, tomato planting, and sampling, bioinformatic analysis and biological and data interpretation; Y.L., X.L. contributed to sampling and data interpretation. J.J. defined the work objectives and technical approach, and contributed to the experimental design, bioinformatics analysis and data interpretation.

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