Mining Heat Stress Associated Genes in Tomato Fruit (Solanum lycopersicum L.) Through RNA-seq

Zhang Ying-ying, Liu Ya-hui, Zhang Hui, and Zhu Wei-min

Shanghai Key Laboratory of Protected Horticulture Technology, The Protected Horticulture Institute, Shanghai Academy of Agricultural Sciences, Shanghai 201403, China

Abstract: Tomato (Solanum lycopersicum L.) is a thermophilic vegetable crop, but sensitive to high temperature stress, especially under the greenhouse conditions. Due to global climate changes, heat stress has now become a great threat to tomato production and fruit quality. Many studies have been conducted to determine the functions of genes in tomato responsive to abiotic and biotic stresses, but transcriptomic information on heat stress responses of tomato fruit is still limited. To investigate heat stress associated genes in tomato fruit, a cDNA library was constructed using fruit harvested from tomato cv. P19-9 plants grown under 42℃ for 0, 1, 2 and 4 h and the expression profiles of heat stress responsive genes in tomato fruit were analyzed through RNA-seq. A total of 632 224 558 clean high quality paired-end reads were obtained and then mapped to reference genome for RNA-seq analysis. After quality control analysis, alignment analysis and transcript assembly, a total of 55 457 RNA transcripts were obtained with functional annotations. Overall, 6 869 differentially expressed genes (DEGs) were identified with a significant response to one or more of the three heat stress treatment times. Based on GO enrichment analysis, 22 genes potentially involved in tomato thermo-tolerance were selected and validated for their expressions through qPCR. The expression profile of tomato fruit genes obtained in this study could shed light on the mechanism and gene editing breeding projects for tomato thermo-tolerance. These findings could also benefit improvement of harvest and storage of tomato in greenhouse.

Key words: tomato fruit, heat stress, transcriptomic analysis, gene expression profiling

Introduction

Tomato (Solanum lycopersicumL.) is a popular vegetable crop with great economic values. Tomato is also a model plant for functional genomics. Although cultivated tomato prefers warm weather (25-33℃ during the day time and 15-20℃ at the night time) (Camejoet al., 2005), it is susceptible to high temperature stress, especially during greenhouse productions. In addition to high temperature stress, tomato is susceptible to other abiotic stresses, such as salinity, drought and chilling (Yuet al., 2019). It was reported that the global warming-induced heat stress has become a severe threat to tomato industry worldwide (Tubielloet al., 2007). For example, high temperature will inhibit tomato development and growth, leading to a significant reduction of tomato yield and quality (Wahidet al., 2007). During summer tomato harvest season of Shanghai City, greenhouse temperature always reaches over 40℃, and after several hours the fruit will be transported in low temperature storage after being picked. As plant heat tolerance can be influenced by many environmental factors, management practices and especially tomato genotypes (Martínet al., 2014), identifications of genes responsible for tomato heat tolerance and their expression profiles can benefit the targeted tomato breeding projects for heat resistance.

In the past few years, the negative impacts of climate change on agriculture have encouraged plant researchers to investigate the mechanism controlling plant tolerance to heat stress (Bokszczanin and Fragkostefanakis, 2013; Crameret al., 2011). It is well known that high temperature stress can regulate the expressions of heat shock protein (HSP) genes. To date, manyHSPgenes, includingHSP100,HSP90,HSP70,HSP60 and smallHSPgenes (smHSPs) (Young, 2010), have been identified in plants and the proteins encoded by theseHSPgenes acting as chaperones are important for protein stability and function. In addition, transcription factors associated with heat stress responses are known to play a key role in signaling pathways (Baniwalet al., 2004).

Numerous studies have been conducted to identify genes involved in various mechanisms, such as oxidative burst, salicylic acid, abscisic acid, ethylene and calmodulin in higher plants grown under heat stress conditions (Kotaket al., 2007; Larkindale and Huang, 2004; Larkindale and Knight, 2002; Suzukiet al., 2011; Zenget al., 2015). Several heat stress related genes have been cloned and their functions have been validated in tomato. For example, the tomato mitogen-activated protein kinase 1 (SlMPK1) gene has been identified, cloned and confirmed to negatively regulate tomato antioxidant defense signaling. Meanwhile, tomato serine-proline-rich protein homolog (SlSPRH1), a substrate ofSlMPK1, has been shown to involve in multiple signaling pathways, including the pathway controlling the response to heat stress (Dinget al., 2018; Dinget al., 2019). SlMAPK3 is a negative regulator of tomato thermo-tolerance, through regulation of expressions of antioxidant enzymes andHSPs/HSFs(Yuet al., 2019). More recently, Wanget al. (2019) have found that over-expression of tomato DnaJ protein 20 (SlDnaJ20) in transgenic tomato plants enhances tomato thermotolerance, whereas suppression ofSlDnaJ20 gene expression in tomato plants increases tomato heat sensitivity. To date, how tomato fruit response to heat stress remains unknown.

Transcriptome analyses are widely applied to be used to systematically investigate the molecular reactions by which plants respond and adapt to complex ambient changes. By RNA-seq analyses, many heat stress responsive genes have been identified inBrassica juncea(Bhardwajet al., 2015),Triticum aestivumL. (Kumaret al., 2015),Brachypodium distachyon(Chen and Li, 2016),Zea maysL. (Shiet al.,2017) andOryza sativaL (Gonzalez-Schainet al., 2016). Although the mechanisms of heat stress response are well studied few information is available regarding heat stress response in tomato fruit. The genetic basis and molecular response to heat stress of tomato fruit have not been investigated.

In this study, a cDNA library was constructed from tomato cv. P19-9 fruits subjected to heat stress and mapped to reference genome for RNA-seq analysis. In response to a time course of heat stress exposure (0, 1, 2 and 4 h), expression profiling analysis of tomato fruit by RNA-seq was performed. Expression changes in response to the HS treatment were confirmed for a subset of the genes. Then multiple differentially expressed genes (DEGs) encoding heat shock factors (HSFs), heat shock proteins (HSPs), ROS-scavenging enzymes and signal transduction- and carbohydrate metabolism-associated enzymes were selected for validations of their expressions through qPCR analysis. The heat stress responsive genes identified in this study could be useful in future tomato breeding programs targeted for heat stress tolerance in fruits.

Materials and Methods

Growth conditions and sample preparation

Seedlings of cherry tomato (Solanum lycopersicumL.) cv. P19-9 were grown in pots inside a greenhouse maintained at 28℃ and 14/10 h (light/dark) photoperiod. Tomato plants at the ripening stages were shifted to 42℃ and ripening tomato fruit is collected from these plants at 0 (refers to HS_0h), 1 (HS_1h), 2 (HS_2h) or 4 h (HS_4h). The three biological replicates were applied at each time point. The heat stress treatments (42℃, 0, 1, 2 and 4 h) were set to simulate environmental conditions in greenhouse of tomato harvest season in Shanghai City. The total RNA was extracted from collected tomato fruit samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). The resulting total RNA samples were checked for RNA concentration and purity using a Nanodrop2000 instrument (Thermo Fisher Scientific, Waltham, MA, USA) and followed by electrophoresis in agarose gels. After analysis of RIN values using an Agilent 2100 Bioanalyzer Lab-on-Chip system (Agilent, Santa Clara, CA, USA), RNA samples were stored at -80℃ for further experiments.

Construction of cDNA library and sequencing

The cDNA library was constructed by Illumina TruseqTM RNA sample prep Kit. The mRNA was isolated by magnetic beads with Oligo (dT) and fragmented randomly by adding fragmentation buffer, of which small fragment of about 300 bp could be separated through magnetic beads screening. Under the action of reverse transcriptase, mRNA was used as template by random hexamers to reverse synthesize a strand of cDNA, followed by two-strand synthesis to form a stable double-stranded structure. Added the End Repair Mix to make double-stranded cDNA flat End, and then add "A" base at the 3' End to connect the Y-shaped connector. The resulting cDNA library was amplified through PCR (15 cycles) and then checked through electrophoresis in 2% agarose gels prior to PCR fragment recovery. TBS380 micro fluorescence (QuantiFluor ST/P, Promega, Madison, WI, USA) was used to visualize cDNA library. The clustering of samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, United States) following the manufacturer's procedure. Then, the constructed cDNA library was sequenced on the Illumina Novaseq 6000 platform using the paired-end protocol. The cDNA library construction and Illumina sequencing were performed by Shanghai Majorbio Bio-pharm Biotechnology Co (http://www.majorbio.com, Shanghai, China).

Data analysis

After sequencing, the raw sequencing reads were filtered using the SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle software (https://github.com/najoshi/sickle) as instructed. Sequencing saturation and coverage for cDNA library was determined using the RSeQC-2.3.2 software as described (Wanget al., 2012). Sequencing adapter sequences, lowquality reads, sequences with high number of N and very short sequences were removed. The remaining high-quality reads (with two or less base mismatches) were mapped to tomato reference genome (ftp://ftp.solgenomics.net/tomato_genome/assembly/build_3.00/, version SL3.0) using the Tophat software (http://tophat.Cbcb.umd.edu/) as described by Trapnellet al(2009). The mapped reads were assembled using the Cufflinks software (http://cufflinks.cbcb.umd.edu/) (Trapnellet al., 2010). The expression quantity of each gene (Transcripts Per Million reads, TPM) was estimated using the RSEM software (http://deweylab.github.io/RSEM/). "False Discovery Rate (FDR)≤0.05 (Benjamini and Hochberg, 1995) and |log2FC (Fold Change)|≥1" were used to determine the significance of each differentially expressed gene (DEG). Gene Ontology (GO, http://www.geneontology.org/) and functional enrichment analysis were conducted to identifyDEGsusing the Goatools software (https://github.com/tang haibao/goatools)(p≤0.05) (Tanget al., 2008).

Validation of gene expression through qPCR

The total RNA was extracted from various tomato fruit samples using MiniBEST Universal RNA Extraction Kit (TaKaRa, Dalian, China). The firststrand cDNA was synthesized with a PrimeScript RT reagent kit (TaKaRa) and used for qPCR. qPCR was conducted on the ABI Prism 7500 Fast Real-time PCR System (Applied Biosystem, Waltham, USA) using SYBR Premix ExTaqII kit (Tli RNaseH plus) as instructed (TaKaRa, Dalian, China). The expression level of tomato eukaryotic initiation factor gene (eiF,Solyc12g096000) was used as an internal control. The relative expression levels of the assayed genes were calculated using the 2-ΔΔCTmethod (Schmittgen and Livak, 2008). All the primers used in this study are listed in Table 1.

Table 1 List of primer sequence

Results

Basic statistical analysis and quality control of RNA-seq

In this study, a non-normalized cDNA library was constructed using an equally mixed pool of four mRNA populations isolated from four fruit samples harvested from tomato plants treated at 42℃ for 0, 1, 2 or 4 h. Fruit harvested from plant grown under 42℃ at 0 h were used as controls. Clean reads were obtained after filtering and conversion of raw reads (Table 2).

Alignment analysis and assembly of transcripts

Samples from three biological replicates representing four temperature treatments, HS_0h (42℃ for 1 h, control), HS_1h (42℃ for 1 h), HS_2h (42℃ for 2 h)or HS_4h (42℃ for 4 h), were sequenced for transcript profiling. Then, these high quality clean reads were aligned and mapped to reference genome (International Tomato Genome Sequencing Project, https://solgenomics.net/). The search results showed that about 95% clean reads were mapped to the coding regions of many genes (Table 2). Gene-body coverage, sequencing saturation, reads distributions of different regions and different chromosomes were analyzed. The results showed that the mapping rates were high and the mapping coverage was satisfactory. Genome-based assembly using the filtered reads data was performed using Cufflinks, identified 55 457 transcripts. Therefore, RNA-seq data were sufficiently high quality to merit further investigation.

Table 2 Summary of RNA-seq results

Identification of differentially expressed genes (DEGs) in response to heat stress (HS)

To identify genes whose expressions were altered after the HS treatments, gene expression profiles from HS_1h, HS_2h, and HS_4h treated plants were individually compared with those from the control plants (HS_0h). According to the criteria ofP-adjust ＜0.05 and the absolute value of log2FC≥1, a total of 6 869 genes (DEGs) were found to be differentially regulated after at least one of the three conditions. For example, 899 genes were found to be significantly up-regulated (P＜0.05) after HS_1h treatment, while 2 571 genes were significantly down-regulated genes. The number ofDEGsreached a peak at HS_1h and was gradually reduced at HS_2h and HS_4h (Fig. 1A). Among the identifiedDEGs, 716DEGswere found at all the three HS time points, when compared to HS_0h treatment (Fig. 1B).

Strongly up-regulated genes (log2FC≥4) and down-regulated genes (log2FC≤-4) were identified from expression profile analysis at the three HS treatment time points. These genes mostly comprised stress response HSPs, transcription factors, ROSrelated genes, metabolism-related genes and signal transduction genes. Genes encoding heat shock proteins includedSolyc03g115230, annotated as heat shock protein and expression of which was 68-fold higher at HS_1h and 56-fold higher at HS_2h than in the 0 h control, andSolyc09g065660, putative heat shock transcription factor A7, log2FC expression value of which reached approximately 6.3 at HS_1h.

A hierarchical clustering (HCL) algorithm was used to generate a heat map visualization of the global clustering pattern of the strongly regulated genes (Fig. 1C). Based on this clustering, most strongly up- and downregulated genes (|log2FC|≥4) appeared after HS_1h treatment followed by gradual declines.

Validation of RNA-seq identified DEGs through qPCR

To validate RNA-seq identifiedDEGs, the 22DEGswere selected and checked for their expressions through qPCR using gene-specific primers (Table 2). These selected genes were involved in heat-shock protein, transcription factor, metabolism-related genes and signal transduction genes. Results of qPCR (Fig. 2) agreed with the results from RNA-seq for theseDEGs(R2=0.777), indicating that the results obtained through RNA-seq were reliable (Fig. 3).

Fig. 1 Analysis of differentially expressed genes in response to heat stress

Gene ontology analysis of significantly regulated genes

Gene ontology (GO) classification of the significantly regulated genes was performed to identify genes regulating tomato fruit HS tolerance. The most prominent 50 GO sub-classifications for the three HS treatments are shown in Fig. 4. Based on these functional classifications, the significantly regulated genes were compared at each of the three time points. Overall, most regulated genes were associated with the categories of biological process, abiotic stimulus response, protein folding and metabolic activity. For example, the expression changes of some genes, including the proteins response to heat (GO: 0009408), proteins response to reactive oxygen species (GO: 0000302), showed similar expression patterns at all the three HS time points. Most of these genes showed an up-regulated expression response during the three HS time periods, while some showed the reverse patterns. The gene group for regulation for biological process (GO: 0050789) demonstrated these diverse patterns. There were 372, 211 and 186 significantly regulated genes at the treatment time points HS_1h, HS_2h and HS_4h, respectively. In this group, 118 genes showed no significant differences in expression levels, when compared at HS_1h and 279 did the same at HS_2h, either. Meanwhile, 84 genes were significantly regulated at all the three time points. These observations demonstrated that the diverse gene expression patterns defined the response to heat stress.

Fig. 2 Quantitative real-time PCR of differentially expressed genes

Fig. 3 Linear regression analysis between quantitative realtime PCR and RNA-seq results (R2=0.777) for 22 genes

Fig. 4 Gene ontology (GO) enrichment analysis of all differentially expressed genes

HSPs in response to HS

Based on the results obtained in this study, 19 upregulatedHSPs(log2FC≥4) were identified in at least one of the three HS treatments. Of these genes, 10 encoded various smallHSPs, including class I and class II DnaJ, five encodedHSP70 forms, two encodedHSP90, one encodedHSP110 and one encodedclpB/HSP110 chaperone. Of these genes, many showed their highest inductions after HS_1h treatment followed by gradual declines, although their expression levels after HS_4h treatment were still much higher than those in HS_0h control samples.

Transcription factors affected by HS treatments

In this study, the expressions of many transcription factors (TFs) of different families were significantly affected by HS treatments. The strongly up- or downregulated TFs (log2FC≥4/≤-4) were mainly in the families of HSF, WRKY, MYB, bHLH, NAC and AP2/ERF. Among these TFs, TFs in WRKY family displayed the greatest changes in their expressions. For example, 18 WRKY members, including WRKY11, WRKY16, WRKY22, WRKY33, WRKY40, WRKY41 and WRKY54, were significantly downregulated, and WRKY33 (Solyc09g014990) and WRKY54 (Solyc08g082110) exhibited the biggest down-regulations with log2FC values of -7.074 and-9.828, respectively. The threeHSFgenes were found to be significantly regulated. The expressions ofHsfA7 (Solyc09g065660) andHsfA2 (Solyc08g062960) were up-regulated after all the three HS treatments, while the expression ofHsf4b(Solyc07g055710) was downregulated by 7.110-fold after HS_1h treatment, and then by 4.023- and 3.067-fold down-regulations after HS_2h and HS_4h treatments, respectively.

ROS-scavenging genes in response to HS

In this study, the genes related to GO term of reactive oxygen species (GO: 0000302) were enriched at each of the three HS treatment time points, and several genes involved in the production or scavenging of ROS were up-regulated in samples from at least one of the treatment times. The expression of Tubby-like F-box protein gene (Solyc09g074510) was up-regulated by 5.457-, 5.410- and 5.059-fold at the three HS treatment time points, respectively. Meanwhile, the expression of ascorbate peroxidase gene (Solyc06g005160) was up-regulated at HS_1h and HS_2h. This induction was not observed at HS_4h. In contrast, the expression of catalase gene (Solyc01g100630) exhibited a significant downregulation through heat stress treatment.

Signal transduction-associated genes affected by HS treatment

HS-regulated genes corresponding to GO category of signal transduction (GO: 0007165) were analyzed. These genes were enriched at all the three HS time points, and 122 genes involved in signal transduction were found to show significant responses. Most genes in this GO term responded to heat stress immediately and showed significant differences of expression levels at the early stage. The auxin-responsive protein IAA2 (Solyc06g084070) showed immediate response to heat stress and the highest degree of down-regulation, whose expression level went down to 6.681-fold at HS_1h and exhibited no significant difference at latter HS time points.Solyc09g065850 (IAA3) andSolyc06g053840 (IAA4) were also found to be downregulated at HS_1h.

It was known that Ca2+/calmodulin (CaM) signaling pathway mediated the heat stress (HS) response and acquisition of thermotolerance in plants. Calciumdependent protein kinases (CDPKs or CPKs) were the unique and key calcium-binding proteins, which acted as a sensor for the increase and decrease in the calcium (Ca2+) concentrations. It was reported that CPKs could regulate plant adaptations to drought, salinity, heat and cold stress environments (Atifet al., 2019). In this study, several genes encoding CaM showed immediate response to heat stress. The gene encoding calmodulin 5 (Solyc12g099990) was found to be up-regulated at HS_1h and HS_2h. In contrast, the expression ofSolyc03g097100 was down-regulated at all the three HS treatment time points and the expression ofSolyc06g053930 was downregulated at HS_1h. The expression of calciumdependent protein kinase gene (Solyc05g056570) was up-regulated, while its homolog calcium-dependent protein kinase gene (Solyc03g033540) was downregulated at HS_1h. Genes encoding MAPKK2 (Solyc03g123800), MAPK3 (Solyc06g005170) and MAPK5 (Solyc01g094960) were down-regulated after HS treatment. In addition, many phytohormone signaling related genes, including ethylene response factor A.2 (Solyc03g093610), auxin response factor 7A (Solyc07g016180) and abscisic acid receptor (Solyc06g050500), were down-regulated.

Carbohydrate metabolic process associated genes affected by HS treatment

GO analysis also indicated that GO term of carbohydrate metabolic process (GO: 0005975) was enriched at all of the three HS treatment time points. For example, many sucrose-related genes were significantly affected by heat stress. Glucan endo-1, 3-beta-glucosidase gene (Solyc03g044085), UDP-glucuronate 4-epimerase 4 gene (Solyc10g018260) and sucrose-phosphate synthase gene (Solyc07g007790) were up-regulated after heat stress treatment, while Glucan endo-1, 3-beta-glucosidase gene (Solyc06g076170), xyloglucan endotransglucosylase/hydrolase gene (Solyc06g083400) and myo-inositol-1-phosphate synthase one gene (Solyc04g050860) were downregulated. One member in the Pectinesterase gene family (Solyc09g075330) was strongly downregulated at all the three HS time points, especially at HS_1h. Two members in Pectinesterase gene family(Solyc03g123620 andSolyc03g078090) were strongly down-regulated at HS_1h and then recovered, while other two members (Solyc01g099960 andSolyc01g098940) were up regulated at all the three HS time points. In this study, the expression of 6-phosphogluconate dehydrogenase gene (Solyc12g040570) exhibited the highest down-regulation level, at log2FC value of -8.860 at HS_4h. It was noteworthy that both strongly up-regulated genes and strongly down-regulated genes associated with carbohydrate metabolic process appeared at HS_1h.

Discussion

Recent release of tomato genome sequence by the Tomato Genome Consortium (Tomato Genome Consortium, 2012) promoted us to investigate tomato fruit responses to heat stress using the tomato reference genome database (SL3.0) and the tomato genome annotation platform (ITAG3.2). Through alignment, more than 95% of the clean reads had perfect matches, allowing us to identify genes in tomato fruits responsive to heat stress.

Heat stress was known to cause protein denaturation, leading to formation of cytotoxic protein aggregates (Suzukiet al., 2011). It was also known that heat shock proteins and other molecular chaperones could prevent or reverse inactivation and aggregation of proteins caused by abiotic and biotic stresses (Doyleet al., 2013). Expressions of heat shock proteins and other molecular chaperones could be influenced by various stresses, including heat stress (Jacobet al., 2017). Here the expressions of many heat shock protein genes in tomato fruit were up-regulated by heat stress. FiveHSP70 genes and twoHSP90 genes were found to be up-regulated and these genes were known to participate in protein folding and unfolding, thus supporting thermo-tolerance, when exposed to heat stress. Meanwhile, there was crosstalk betweenHSP90 andHSP70 chaperones and heat stress transcription factors in tomato (Hahnet al., 2011). In this study, oneHSP110 gene and oneclpB/HSP100 gene were also found to be up-regulated. These two genes had been reported to regulate stress responses in plant (Mishra and Grover, 2016). Many smallHSPswere also upregulated after HS treatments. These smallHSPshad been shown to be important for plant thermotolerance and could protect the functions of protein translation factors during heat stress (McLoughlinet al., 2016).

InArabidopsis, AtHsfA1 had been shown to play a pivotal role in heat stress-induced transcriptions of various genes, includingsmallHSPs,HSP70s,HSP101sand heat shock transcription factors (e.g., HsfA2, HsfA7a, HsfB1 and HsfB2a) and genes encoding heat stress-induced metabolic enzymes (e.g., Ips2 and GolS1) (Liuet al., 2011). In this study,SolycHsfA2 andSolycHsfA7 were up-regulated after HS treatment. In addition, the results suggested the that other transcription factor genes were involved in tomato fruit heat stress responses. Among these transcription factor genes, the expressions of genes in WRKY, MYB and bHLH families were mostly downregulated after HS treatment.

Numerous studies had shown that during heat stress, excessive ROS could accumulate in cells and cause damages to cell membrane, DNA and proteins. In response to the increased amount of ROS in cells, plant antioxidant enzyme synthesis pathways could be activated (Milleret al., 2008; Suzukiet al., 2011). Consistent with these reports, the expressions of ROS-scavenging genes were founded to be induced in this study. For example, the expressions of Tubby-like F-box protein gene (Solyc09g074510) and an ascorbate peroxidase gene (Solyc06g005160) were substantially up-regulated after HS treatments.

Signal transduction was also known to play important roles in plant heat stress responses. Several signaling pathways could activate a similar set of heat stress responsive genes to enhance plant thermotolerance (Kotaket al., 2007; Mittleet al., 2012). For example, transient changes of intracellular Ca2+concentrations had been shown to function as signals to trigger certain physiological responses to allow plant to couple with various environmental stimuli. Zenget al. (2015) reported that calmodulin (CaM) and calcium-dependent protein kinases (CPKs) were major Ca2+sensors that could recognize encrypted Ca2+signals. Plant plasma membrane could sense temperature changes and induce a transient opening of Ca2+channels during heat stresses (Garget al., 2014). The activation of CaM could stimulate the DNAbinding activity of heat shock transcription factors, leading to the productions of more HSPs to confer a thermo-tolerance (Jiaet al., 2014). In this study, multiple heat stress-affected CaM and CPK genes were known to involve in plant heat stress responses. For example, the expressions ofSolyc12g099990 andSolyc03g097100 were up- and down-regulated, respectively. Two CPK genes (Solyc05g056570 andSolyc03g033540) were also oppositely affected by HS treatments.

Mitogen-activated protein kinases, which belonged to multigene families, could effectively transmit specific stimuli and regulate the antioxidant defense pathway during responses to stress signaling (Sinhaet al., 2011). The repressions ofHAMKandStMPK1 in alfalfa and potato could be affected by heat stresses (Sangwanet al., 2002). In rice, the expression ofOsMSRMK2 could be rapidly induced during heat stresses (Agrawalet al., 2002). Over-expression of maizeMAPK1 (ZmMAPK1) in transgenicArabidopsisplant enhanced its heat stress resistance (Wuet al., 2015). Recent reports had indicated that SlMPK1 and SlMAPK3 were negative regulators of tomato heat stress tolerance (Dinget al., 2018). In this study, the expressions ofMAPKK2 (Solyc03g123800),MAPK3 (Solyc06g005170) andMAPK5 (Solyc01g094960) were decreased after HS treatments, indicating that genes associated with MAPK cascades were likely to participate in tomato fruit heat stress tolerance.

High temperature over 40℃ in green house become a severe problem for harvesting tomato in Shanghai City. The picked fruits could be temporally stored in greenhouse for hours, depending on the level of management and labor. These findings obtained through this study were consistent with previous reports, and were thus useful for the mechanism and future molecular marker-assisted breeding and gene editing breeding for tomato fruit heat stress tolerance. In addition, this study could help to improve harvest and storage of tomato fruits in greenhouse.

Conclusions

In this study, the fruits harvested from tomato cv. P19-9 plants grown under 42℃ for 0, 1, 2 and 4 h were used for RNA-seq analyses for expression profiling of heat stress responsive genes. More than 632 Mb clean high quality paired-end reads were obtained and mapped to reference genome for RNA-seq analysis. After quality control analysis, alignment analysis and transcript assembly, around 55 kb RNA transcripts were obtained with functional annotations. Overall, 6 869 differentially expressed genes (DEGs) showed a significant response at least one of the three heat stress treatment times. Based on GO enrichment analysis, 22 genes potentially involved in tomato thermo-tolerance were selected and validated for their expressions through qPCR. Genes encoding heat shock proteins and transcription factors showed significant expression changes in response to heat stress. In addition, through GO enrichment, GO term of reactive oxygen species, signal transduction and carbohydrate metabolic process were enriched at the three time points. This study provided valuable information for studies on heat stress response in tomato.