Comparative transcriptome analysis reveals gene network regulation of TGase-induced thermotolerance in tomato

Transglutaminase (TGase), the ubiquitous protein in plants, catalyzes the post-translational transformation of proteins and plays a vital role in photosynthesis. However, its role and mechanism in tomato subjected to heat stress still remain unknown. Here, we carried out a transcriptomic assay to compare the differentially expressed genes (DEGs) between wild type (WT) and TGase overexpression (TGaseOE) plants employed to high-temperature at 42 °C and samples were collected after 0, 6, and 12 h, respectively. A total of 11,516 DEGs were identified from heat-stressed seedlings, while 1,148 and 1,353 DEGs were up-and downregulated, respectively. The DEGs upon high-temperature stress were closely associated with the pathways encompassing protein processing in the endoplasmic reticulum, carbon fixation, and photosynthetic metabolism. In addition, 425 putative transcription factors (TFs) were identified, and the majority of them associated with the bHLH, HSF, AP2/ERF, MYB, and WRKY families. RNA-seq data validation further confirmed that 8 genes were linked to protein processing and photosynthesis, and the mRNA level of these genes in TGaseOE was higher than that in WT plants, which is consistent in transcriptome results. In conclusion, these results reveal the transcriptional regulation between WT and TGaseOE in tomato under heat stress and shed light on a new dimension of knowledge of TGase-mediated thermotolerance mechanism at the molecular level.


Introduction
The temperature of the earth surface has gradually increased over the last few decades as a result of anthropogenic activities, and it is predicted that it may rise by an average of 2 to 4.5 °C by 2100 (Bita and 2 Gerats, 2013;Wei et al., 2019). It is well known that high temperature severely retardation of plant growth by alteration of the photosystem, production of excess reactive oxygen species (ROS), and break-down of the membranes stability and proteins (Jahan et al., 2019a;Jahan et al., 2019b;Kim and Lee, 2019;Tao et al., 2020).
It has been projected that crop yields in Africa and Asia will be decreased by 15-35% and in the Middle East by 25-35% if the temperature increases by 3-4 °C (Ortiz et al., 2008). Therefore, high temperature is a severe threat to food security to all over the world. For this reason, it is essential to study the mechanisms of plant adaptation to high temperature stress.
Transglutaminases (TGases) are ubiquitous in all plant tissues, catalyzing the post-translational alteration of proteins via constant intra-and inter-molecular connections between glutamine residue and/or free polyamine or lysine residue (Serafini-Fracassini et al., 2009). This plays a major role in the regulation of photosynthesis, fertilization, senescence, and response to environmental stimuli (Aloisi et al., 2016). AtPng1p, the first identified TGase in plants, contains the Cys-His-Asp triad of the TGase catalytic domain and widely distributes in plants (Della Mea et al., 2004). TGases identified in maize and rice, and their expression depends on light and relevant to plastidial protein-mediated photoprotection and the thylakoid electrochemical gradient (Campos et al., 2010;Campos et al., 2013). Proteomic analysis showed that TGase interacts with LHCII, ATPase, and PsbS proteins (Campos et al., 2010;Campos et al., 2013). TGase is a vital component of the photosystem II (PSII) protein complex and regulating photosynthesis. Indeed, overexpression of TGase from maize in Arabidopsis extends the stroma thylakoid network, increases linear electron flow and threshold of photoprotection (Ioannidis et al., 2016). The initial and activation status of RuBisCO, RuBisCO activity and FBPase activity (Fructose-1, 6-bisphosphatase) increases with an increase of TGase activity, resulted in elevating the net photosynthetic rate (Zhong et al., 2019a). Furthermore, TGase interacts with cytoskeleton regulates the apical growth of pollen tube and mediates the process of self-incompatibility (Del Duca et al., 2013;Mandrone et al., 2019). Importantly, under saline condition, the activity of TGase is elevated, and the salt-tolerant variety has higher TGase activity than the salt-sensitive variety in cucumber . Overexpression of TGase in tomato enhanced salt tolerance through photosynthesis coordination . However, the exact mechanism in other abiotic stimuli, in particular, high temperature is still not clear. In the current study, a comparative transcriptome analysis was used to reveal the gene network regulation of TGase-induced thermotolerance in tomato.

Materials and Methods
Plant materials and treatment Tomato (Solanum lycopersicum L. cv. Ailsa Craig) seedlings were used as test material in this experiment. Two independent homozygous TGase overexpressed (TGaseOE) lines were obtained from our laboratory (Zhong et al., 2019a). Germinated seeds were sown in a 50-holes seedling tray filled with peat and vermiculite mixture (2:1, v:v) and transferred to it artificial growth chamber. After the second true leaf fully developed, the seedlings were shifted into plastic pots (250 cm 3 ) filled with same growth medium as described above and placed in the growth chamber, where diurnal temperature 25 °C/20 °C (day/night), relative humidity, 65-75%, 14 h photoperiod at 300 μmol m -2 s -1 photosynthetic photon flux density (PPFD) were maintained. The seedlings were watered with half-strength of Hoagland nutrient solution every two days.
Six-week-old tomato seedlings were transferred into growth chambers and exposed to high temperatures (25 °C or 42 °C). Leaf samples were collected at 0, 6, and 12 h from control and heat-treated plants and immediately kept in liquid nitrogen, followed by storage at -80 °C for subsequent gene expression and RNAseq (For RNA-seq, the collected samples were sent to company in the second day).
3 Electrolyte leakage measurement Electrolyte leakage was estimated following Hong et al. (2003) described method. The electrolyte leakage was calculated with the following formula: L1 is the initial electrical conductivity after keeping in 2 h at room temperature, and L2 is the final electrical conductivity after 15 min boiled the leaves at 95 °C.
Net photosynthetic rate determination Tomato seedlings exposed to heat-stressed for 24 h used to determine the net photosynthetic rate (Pn). The Pn measurement was done using portable photosynthesis system (LI-6400; Li-COR, Lincoln, NE, USA), where CO2 concentration maintained at 380 μmol mol -1 , and photosynthetic photo flux density (PPFD) at 800 μmol m -2 s -1 .
RNA isolation, library construction and transcriptome sequencing For RNA extraction, leaf samples were collected from WT and TGaseOE-1# plants at different time points and homogenized into liquid nitrogen. The total RNA was isolated using Trizol ® Reagent (Invitrogen, USA) following the manufacture's reference manual. RNA concentration was determined employing the NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA) and the integrity of RNA was verified applying the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The Transcriptomic sequencing libraries were created using the NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) in accordance with the manufacturer's reference manual, and index codes were included to mark up sequences for an individual sample. The library quality was verified on the Agilent Bioanalyzer 2100 system. The clustering of index-coded samples was executed on a cBot Cluster Generation System employing the TruSeq PE Cluster Kit v3-cBot-HS (Illumia) in accordance with the company's protocol and the library preparations were sequenced on an Illumina Hiseq platform (Beijing Biomics Biotechnology Co., Ltd., Beijing, China) followed by constructed paired-end reads.

RNA sequence data analysis
The clean reads were generated from the raw reads (containing ploy-N) after removing the low-quality reads, mismatches, and adaptor sequences. The reference genome sequence Solanum lycopersicum (SL3.0) and annotation files (ITAG3.20) were used and downloaded from the tomato genome (http://solgenomics.net/) website. Tomato reference genome index was performed employing Bowtie v2.2.3 software (Langmead and Salzberg, 2012). The alignment of paired in clean reads were made through TopHat v2.0.12 software (Trapnell et al., 2009), and the transcript abundance was computed according to fragments per kilo-base of exon per million fragments mapped (FPKM). Differential genes expression (DEGs) analysis was performed using the DESeq R package (1.10.1) (Wang et al., 2010). To analyze the RNA-Seq data, the resulting P values were modified using the Benjamini-Hochberg approach for controlling the false discovery rate (FDR) and DEGs were calculated according to log 2 fold change≥2, a FDR < 0.01 and adjusted p-value <0.01. The DEGs were employed for gene ontology (GO) analyses, and all DEGs were mapped to GO terms using the GO seq R package. The statistical enrichment of DEGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with DEGs were analyzed using KOBAS 2.0 (https://www.biostars.org/p/200126/) database at the threshold of FDR ≤ 0.05.

Quantitative real-time polymerase chain reaction (qRT-PCR) used for the validation of gene expression
For the confirmation of the RNA-seq data, the transcript abundance of selected up-and down-regulated genes was further validated through qRT-PCR. The total RNA was isolated using Trizol® Reagent (Invitrogen, 4 USA) following the manufacture's reference manual. The total RNA was reverse transcribed into cDNA using the HiScript II Q RT SuperMix (+gDNA wiper) Kit (Vazyme, Nanjing, China) for qRT-PCR assays. The qRT-PCR was carried out following the ChamQ SYBR qPCR Master Mix (Vazyme, Nanjing, China), and 20 μL reaction system was used for amplification using the LightCycler® 480 II Real-Time PCR detection system (Roche, Basel, Switzerland). Gene-specific primers were designed (primer premier 5.0 software) based on the nucleotide sequences retrieved from NCBI (the National Center for Biotechnology Information) database and tomato genome sequence (http://solgenomics.net/) database (Table S1). Actin was used as an internal control for the normalization of transcript abundance levels. The relative transcript abundance levels were estimated using the 2 −ΔΔCT method in accordance with Livak and Schmittgen (2001).

Statistical analysis
All treatment data were analyzed using SPSS statistical software (IBM SPSS 20, Chicago, IL, USA). Significant differences (P < 0.05) between treatments were estimated according to Tukey's test, and at least three independent biological replications were performed for each determination.

Results
Overexpression of TGase enhanced heat tolerance in tomato To elucidate the role of TGase-induced heat tolerance in tomato, wild type (WT) and TGaseoverexpressed (TGaseOE) plants were used. Apparently, under normal growth conditions, there was no statistical variation observed in electrolyte leakage (EL) between WT and TGaseOE plants (Figure 1a). In case of WT, plants exposed to heat stress for 24 h, the EL value was markedly increased by 151.63%, while the TGaseOE-1# and TGaseOE-2# plants slightly increased by 23.98-and 75.42%, respectively (Figure 1a), suggesting that TGaseOE plants were more ability to heat tolerance than WT plants. In addition, TGaseOE plants had higher net photosynthetic rate (Pn) than WT under the control conditions ( Figure 1b). Although, heat stress significantly inhibited the Pn value in TGaseOE plants, the Pn value was still more than that in WT plants ( Figure 1b). Taken together, the above results suggest that overexpression of TGase enhanced heat tolerance by maintaining membrane integrity and holding higher Pn value. To investigate the detailed mechanism of TGase-induced heat tolerance in tomato, a total of 18 cDNA libraries were constructed from WT and TGaseOE-1# leaves at 0, 6, and 12 h after heat treatment and also identified genes responsible for heat tolerance in tomato. After removing reads containing adapter, low quality reads, and ploy-N from raw data, a total of 382,856,771 sequence reads and 114.84 Gb of sequence data were generated from the designed sample (Table S2). The ratio of clean reads to raw reads was over 96.5%, the Q30 was more than 91.2%, and the GC content was approximately 42% for each library (Table S2). Therefore, the quality and integrity of the sequence data were sufficient for further analysis. Moreover, to get similarity index of the reads, we performed gene alignment using tomato genome database, and the results indicated that over 88.9% of the reads were identical to the tomato genome, and more than 99.9% of the reads were uniquely mapped (Table S3).  Table S4). Venn diagram assay presented 6 that 11,516 DEGs were detected, while, 1,148-and 1,353 unigenes were up-and down-regulated, in WT and TGaseOE-1# plants, respectively under heat stress treatment for 0, 6, and 12 h (Figure 2b, c, d).

Gene ontology and pathway enrichment analysis of heat stress modulate DEGs
We applied to GO enrichment pathways analysis to categorize the functions of the 1,148 up-regulated and 1,353 down-regulated DEGs. The enriched genes were annotated into three fundamental GO classifications in particular, biological process, cellular component and molecular function. GO enriched pathways results showed that the top GO terms (18) were found belongs to the biological process, 13 GO terms associated with cellular components and 14 GO terms related to molecular function (Figure 3). The significantly enriched GO terms involved in the biological process were associated in the metabolic process (GO: 0008152), cellular process (GO: 0009987), single-organism process (GO: 0044699), biological regulation (GO: 0065007) and response to biotic and abiotic stress (GO: 0050896) (Figure 3). Cellular compartment GO enrichment analysis showed that they are mostly related to cell part (GO: 0044464), membrane part (GO: 0044425), organelle (GO: 0043226), organelle part (GO: 0044422) and macromolecular complex (GO: 0032991) (Figure 3). Furthermore, molecular function encoded classified GO terms highly enriched in catalytic activity (GO: 0003824), transporter activity (GO: 0005215), binding (GO: 0005488), and nucleic acid binding transcription factor activity (GO: 0001071) (Figure 3). To identify the biological functions and signal transduction pathways of the commonly regulated DEGs under high temperature stress, the highly enriched pathways were identified through applying KEGG analysis. The results have shown that 1,148 up-regulated DEGs were involved in 73 typical KEGG pathway, encompassing in protein processing in the endoplasmic reticulum, oxidative phosphorylation, plant-pathogen interaction, spliceosome, circadian rhythm-plant, ubiquitin-mediated proteolysis and glutathione metabolism (Figure 4a and Table S5). A total of 102 distinct KEGG pathways were involved in the down-regulation of DEGs (Table S6) and the highly enriched pathways included are carbon metabolism, carbon fixation in photosynthetic organisms and DNA replication (Figure 4b). Furthermore, 23 DEGs were identified, which induced high temperature and associated with protein processing in the endoplasmic reticulum pathway, including HSP20, HSP70 and HSP90 (Table S5). Importantly, the transcript abundance of most up-regulated genes associated with protein processing in the endoplasmic reticulum pathway and their expression were higher in TGaseOE-1# plants than WT (Table S7). Conversely, heat stress significantly inhibited the carbon metabolic-related gene expression, and the repression effect was prominent in WT plants compared to TGaseOE-1# plants (Table S8). Transcription factors (TFs) have a significant role in response to high-temperature stress in plants Ohama et al., 2017). To gain a better understanding of the TGase-mediated heat tolerance mechanism, we analyzed differentially expressed TFs. A total of 186 (100 up-regulated, 86 down-regulated), 200 (110 up-regulated, 90 down-regulated), 191 (80 up-regulated, 111 down-regulated) and 234 (108 upregulated, 126 down-regulated) differentially expressed TFs were detected in the sample of WT-0 h vs WT-6 h, WT-0 h vs WT-12 h, OE1-0 h vs OE1-6 h and OE1-0 h vs OE1-12 h, respectively (Figure 5a). The number of suppressed TFs in 6 h was lower than 12 h, which might be due to heat stress duration begin to inhibit plant response. Amongst the TFs, most of the heat stress induce isolated genes were belonging to AP2/ERF, bHLH, HSF, MYB and WRKY families ( Figure S1). Venn diagram analysis revealed that a total of 55 (29 unigenes upregulated and 25 unigenes down-regulated) TFs were observed in WT and TGaseOE-1# plants at 0, 6, and 12 h after heat stress (Figure 5b, c, d). These results collectively suggested that the differentially expressed TFs might be involved in regulation of heat stress in response to tomato.

Validation of DEGs via quantitative real-time PCR
For the verification of reliability of the expression trends of DEGs found from RNA-seq data, we selected highly expressed 8 DEGs for qRT-PCR analysis. Among these 8 DEGs, 4 genes were associated with protein processing in the endoplasmic reticulum, and another 4 DEGs were involved in carbon fixation in photosynthetic organisms. The transcript abundance results were shown that in the case of protein processing in endoplasmic reticulum pathway-related genes, including HSP90 (Solyc03g007890, Solyc06g036290), 9 HSP70 (Solyc03g117630) and HSP17.4 (Solyc03g123540) were significantly induced by heat stress, but the transcript abundance pattern of those genes in TGaseOE-1# plants were higher than that observed in WT plants ( Figure 6). In contrast, the expression pattern of carbon fixation associated DEGs, including Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (Solyc02g020940, Solyc04g009030), Ribulose bisphosphate carboxylase small chain (Solyc02g085950) and FBPase (Solyc05g052600) were significantly suppressed both in WT and TGaseOE-1# plants under heat stress, where, TGaseOE-1# plants still maintained higher expression level ( Figure 6). In general, the qRT-PCR results confirmed that these findings were in line with the RNA-seq findings, implying that the results of RNA-seq data was verified and reliable. Figure 6. Verification of the differentially expressed selected DEGs by qRT-PCR Each histogram denotes a mean±SE of three biological replicates (n=3). Means with the difference represent significant differences at P< 0.05 according to Tukey's test. 1#, 1 line of TGaseOE plants.

Discussion
Plants are sessile organisms, coping with harsh environments in particular heat stress evolved different mechanisms at physiological, biochemical and molecular levels. Heat stress-induced growth retardation is common and causes a profound loss of yield (Katano et al., 2018). Herein, we extensively studied the heat stressinduced mechanism at the transcriptomic level using WT and TGaseOE tomato plants. The results revealed that heat stress significantly increased EL level and decreased Pn value, and this trend is severe in WT than TGaseOE plants, indicating that TGaseOE plants are more capable of coping up under heat stress conditions ( Figure 1). To investigate the molecular mechanism of TGase-mediated thermotolerance of tomato, we used transcriptomic analysis to compare WT and TGaseOE plants tolerance mechanisms under heat stress at various time points and a total number of 11,516 DEGs were identified. GO analysis revealed that DEGs induced by heat stress and these enriched GO associated with protein processing in the endoplasmic reticulum, where most of the DEGs were up-regulated. The carbon metabolism-related DEGs, where genes were highly enriched down-regulated DEGs (Figure 4), which suggested that these two pathways might play a critical role in TGase-induced heat tolerance.
Heat shock proteins (HSPs) act as a molecular chaperones under heat stress environment to prevent the denaturation and aggregation of proteins and promotion of refolding and degradation of misfolded proteins as well (Izumi, 2019;Liu et al., 2019;McLoughlin et al., 2019). KEGG pathway analysis showed that under heat stress, numerous DEGs were up-regulated during protein processing in the endoplasmic reticulum pathway, and genes encoding HSP were induced by heat stress, including HSP20, HSP70 and HSP90 (Figure 4 and Table S7). Importantly, the transcript abundance of HSPs in TGaseOE plants was highly up-regulated than WT plants (Table S7). The present findings were consistent with the previous transcriptomic analysis and 10 noted that heat stress markedly induced HSPs gene expression, including HSP70, HSP90 and small HSPs (Frank et al., 2009;Liu et al., 2012;Qian et al., 2019). Furthermore, HSP70 and HSP90 are two important heat response factors that are essential for inhibiting aggregation and supporting the refolding of denatured proteins under stress (Li et al., 2015;Liu et al., 2012). Indeed, overexpression of HSP70 enhanced thermotolerance in Arabidopsis (Sung and Guy, 2003). These results suggested that TGase might promote the expression of HSPs to maintain the integrity of proteins and refold the denatured proteins to alleviate toxicity caused by misfolded and aggregated proteins in cells, thereby increasing heat tolerance in tomato plants.
Photosynthesis machinery and photosystem are distorted under high-temperature stress (Jahan et al., 2019b;Kim and Lee, 2019). In agreement with previous studies, heat stress significantly decreased Pn value both in WT and TGaseOE plants, while the value Pn in WT plants was lower than that in TGaseOE plants in presence or absence of heat stress (Figure 1b). TGase positively regulates photosynthesis and some photosynthesis-related proteins, indicating that TGase is a key regulator on photosynthesis (Campos et al., 2010;Campos et al., 2013). Moreover, the GO term enriched analysis showed that down-regulated DEGs were enriched in carbon metabolism, fixation in photosynthetic organisms, photosynthesis-antenna proteins, including GAPDH and FBPase Chlorophyll a/b binding protein (CAB) (Figure 4b and Table S8). All these carbon metabolisms and photosynthesis-related genes were suppressed under heat stress, their expression in TGaseOE plants were still higher than that of WT plants. GAPDHs are ubiquitous proteins that mediate carbon fixation through catalyzing glycerate-3-phosphate conversion to glyceraldehyde-3-phosphate to synthesize starch in chloroplasts (Chang et al., 2015;. The overexpression of GAPDH led to improve stress tolerance in Arabidopsis (Chang et al., 2015). FBPase enzyme contributes a key role for the regeneration of ribulose-1,5-bisphosphate (RuBP) in the Calvin cycle. Our recent research showed that TGaseOE plants significantly elevated both FBPase activity and its encoding gene expression and the speed of electron transport for the production of RuBP is more in comparison to WT plants under optimal growth conditions (Zhong et al., 2019a). Similarly, our results showed that TGaseOE plants had higher expression level of FBPase both in control and heat stress (Table S8). Furthermore, CAB is considered as the lightharvesting Chlorophyll a/b binding protein (LHC) in higher plants (Silva et al., 2016). CAB proteins strongly linked with photosystems (PS) I and PSII and they are recognized as Lhca (or LHCI) and Lhcb (or LHCII), respectively, and they also act as the core proteins in the light harvest complex for harvesting light and transfer energy to photosynthesis system (Engelken et al., 2012;Velez-Ramirez et al., 2014). Therefore, TGaseOE plants might maintain higher Pn by ameliorating the decline of GAPDH, CAB and FBPase to promote light harvesting and electron transport, resulting in elevating heat tolerance.
TFs are essential key elements played a predominant role in response to environmental stimuli in plants (Wani et al., 2018;Castelan-Munoz et al., 2019). A total of 425 putative TFs were isolated from our study, and most of them belong to AP2/ERF, bHLH, HSF, MYB and WRKY families, and these findings also supported by the previous studies (Li et al., 2015;Qian et al., 2019). Heat shock factors (HSFs) are the key regulators in thermotolerance, and these are binds with heat shock elements (HSE) of a desire gene promoters to regulate their expression, such as HSPs (Guo et al., 2016;Wang et al., 2015). RNA-seq data showed that HSFs were induced by heat stress, including HSFA2, HSFA3, HSFA4a, HSFA4c and HSFB2a ( Figure S1 and Table S4). Strikingly, the expression of heat-induced HSFs in TGaseOE plants was obviously greater compared to WT plants, which was consistent with the expression patterns of HSPs, suggesting that HSFs might execute leading roles in enhancing heat resistance in TGaseOE plants. In addition, AP2/ERF family is plant-specific TFs that are facilities optimum plant growth and development and play a pivotal role in response to various environmental stimuli (Thirugnanasambantham et al., 2015;Klay et al., 2018). Overexpression of DREB6, a type of AP2/ERF TF, induces the transcript of HSFA4, HSP90 and ROS scavenging genes and increases thermotolerance in chrysanthemum (Du et al., 2018). Moreover, ERF74 directly regulated the expression of RbohD to induce the production of ROS, which act as a signal molecule to active stress marker genes in Arabidopsis (Yao et al., 2017). Transcriptomic analysis showed that numerous of AP2/ERF TFs were regulates by high temperature ( Figure S1). Collectively, our RNA-sequence data inferred that TFs might contribute a significant role in TGase-induced heat tolerance in tomato plants.

Conclusions
In summary, to elucidate the molecular mechanism of TGase-induced thermotolerance in tomato, we performed a comparative transcriptomic study between WT and TGaseOE plants. Transcriptome data analysis identified a total of 11,516 DEGs, where 1,148 and 1,353 DEGs were commonly up-and downregulated, respectively. In addition, KEGG pathway evaluation revealed that these DEGs were predominantly linked with protein processing in endoplasmic reticulum, oxidative phosphorylation, carbon metabolism, carbon fixation in photosynthetic organisms. All these enriched pathways and their related genes may be involved in TGase-induced heat tolerance of tomato. Therefore, these results reveal the transcriptional regulation between WT and TGaseOE in tomato exposed to heat stress, which confers a comprehensive understanding at the transcript level mechanism of TGase-mediated thermotolerance.

Authors' Contributions
YW conceived and designed the experiments. MSJ, ZS, and MZ conducted the experiments and wrote the original manuscript. YZ, RZ and MME contributed supplied materials, data collection, analysis tools, analyzed the data. JS SS, and SG supervised and revised the manuscript.
All authors read and approved the final manuscript.