Introduction
Global warming is causing climate change, including extremely high temperatures that lead to devastating reductions in crop production. Transcriptome studies on the effect of heat stress (HS) on several crops have been reported (Kumar et al., Reference Kumar, Goswami, Kala, Rai, Rai, Chinnusamy and Rai2015; Wu et al., Reference Wu, Taohua, Gui, Xu, Li and Ding2015; Chen and Li Reference Chen and Li2017). Heat-responsive genes are often implicated in pathways involving protein homeostasis, transcriptional regulation, antioxidant activities and gene expression (Wang et al., Reference Wang, Liu, Tian, Huang, Shi, Dai and Zhang2017).
Transcription factors (TFs) including heat shock factors (HSFs), MYB, WRKY and ethylene-responsive factors (ERF) are critical regulators of the plant response to abiotic stress (Liu et al., Reference Liu, Yu, Zhang, Wang, Cao and Cheng2013; Guo et al., Reference Guo, Liu, Ma, Luo, Gong and Lu2016; Bai et al., Reference Bai, Sunarti, Kissoudis, Visser Richard and van der Linden2018; Klay et al., Reference Klay, Gouia, Liu, Mila, Khoudi, Bernadac, Bouzayen and Pirrello2018). In Brachypoduim distachyon, 454 TF genes belonging to 49 gene families were differentially expressed upon HS, including ERF and WRKY genes (Chen and Li Reference Chen and Li2017). In plants, one of the most important functions of WRKYs appears to be the regulation of responses to various biotic and abiotic stresses (Li et al., Reference Li, Fu, Chen, Huang and Yu2011; Hu et al., Reference Hu, Dong and Yu2012).
ERFs are a class of TFs that play a role during plant growth, metabolism and in responses to various stresses (Pan et al., Reference Pan, Seymour, Lu, Hu, Chen and Chen2012). Notably, ERFs play a critical regulatory role in enhancing plant abiotic stress tolerance (Song et al., Reference Song, Li and Hou2013). The over-expression of ERF-like TFs can promote plant stress resistance (Gilmour et al., Reference Gilmour, Sebolt, Salazar, Everard and Thomashow2000; Hsieh et al., Reference Hsieh, Lee, Yang, Chiu, Charng, Wang and Chan2002). Although a functional link between ERFs and various abiotic stresses is known for several crops, the involvement of ERFs in HS has not yet been highlighted (Debbarma et al., Reference Debbarma, Sarki, Saikia, Boruah, Singha and Chikkaputtaiah2019).
Along with WRKY and ERF TFs, the MYB TFs represent one of the largest and functional TF families in the genome (Riechmann et al., Reference Riechmann, Heard, Martin, Reuber, Jiang, Keddie, Adam, Pineda, Ratcliffe, Samaha, Creelman, Pilgrim, Broun, Zhang, Ghandehari, Sherman and Yu2000). Diverse MYB proteins in plants have been shown to act in a wide range of stress-related responses (Mmadi et al., Reference Mmadi, Dossa, Wang, Zhou, Wang, Cisse, Sy and Zhang2017). In tomato, MYB TFs are involved in gene regulation (Li et al., Reference Li, Peng, Tian, Han, Xu and Yao2016), and reported to act in abiotic stress responses (Cao et al., Reference Cao, Zhang, Wang, Zhang and Hao2013; Liu et al., Reference Liu, Yu, Zhang, Wang, Cao and Cheng2013).
HSF TFs belong to a complex transcriptional network. They are part of a signal transduction chain that activates HS-candidate genes belonging to plant species as tomato (Song et al., Reference Song, Liu, Duan, Liu, Huang, Ren, Li and Hou2014; Xue et al., Reference Xue, Sadat, Drenth and McIntyre2014; Guo et al., Reference Guo, Liu, Ma, Luo, Gong and Lu2016). Plant HSFs are divided into the A, B and C classes (Nover et al., Reference Nover, Bharti and Döring2001). Class A HSFs act as transcriptional regulators during abiotic stress responses (Shim et al., Reference Shim, Hwang, Lee, Lee, Choi, An, Martinoia and Lee2009; Anfoka et al., Reference Anfoka, Moshe, Fridman, Amrani, Rotem and Kolot2016). HsfA1s are considered as master regulators of the response to environmental stresses, and are directly involved in the expression level regulation of HsfA2 and HsfB1 expression level (Yoshida et al., Reference Yoshida, Ohama, Nakajima, Kidokoro, Mizoi, Nakashima, Maruyama, Kim, Seki, Todaka, Osakabe, Sakuma, SchöZ, Shinozaki and Yamaguchi-Shinozaki2011). Tomato HsfA1-mutant plants are HS-sensitive (Mishra et al., Reference Mishra, Tripp, Winkelhaus, Tschiersch, Theresa, Nover and Scharf2002) and the Arabidopsis HsfA2-knockout is thermosensitive (Charng et al., Reference Charng, Liu, Liu, Chi, Wang, Chang and Wang2006). Class B HSFs modulate HsfA1 activity by acting either as repressors in Arabidopsis or co-activators in tomato (Shim et al., Reference Shim, Hwang, Lee, Lee, Choi, An, Martinoia and Lee2009; Ikeda et al., Reference Ikeda, Mitsuda and Ohme-Takagi2011; Zhu et al., Reference Zhu, Thalor, Takahashi, Berberich and Kusano2012). Another feature of HsfA2 as a key regulator of stress responses is that it promotes anti-oxidant genes expression, mainly via the upregulation of genes encoding ascorbate peroxidase enzymes (APXs) (Panchuk, Reference Panchuk2002; Charng et al., Reference Charng, Liu, Liu, Chi, Wang, Chang and Wang2006). The APXs are multigene family and act as ROS-scavenging enzymes by catalysing the conversion of H2O2 into H2O with ascorbate as an electron donor. APX1 is known to connect the antioxidant pathways of common cellular functions in response to HS (Storozhenko et al., Reference Storozhenko, De Pauw, Van Montagu, Inzé and Kushnir2002). Conversely, few data concerning the regulation of the second cytosolic isoform APX2 expression are reported (Panchuk, Reference Panchuk2002).
HS modulates the expression of multiple proteins including HSPs. They act as molecular chaperones and are involved in the maintenance and restoration of protein homeostasis (Kotak et al., Reference Kotak, Larkindale, Lee, Koskull-Döring, Vierling and Scharf2007). HSPs are divided into five major classes: HSP100, HSP90, HSP70, HSP60, HSP20 and the small HSPs (Iba, Reference Iba2002). HSP101 and HSP17 were known to protect plant cells from heat-induced programmed cell death (Rikhvanov et al., Reference Rikhvanov, Gamburg, Varakina, Rusaleva, Fedoseeva, Tauson, Stupnikova, Stepanov, Borovski and Voinikov2007). During all stages of HS, the expression of HSPs seems to be regulated by heat shock transcription factors (HSFs) on a transcriptional level (Lin et al., Reference Lin, Jiang, Chu, Tang, Zhu and Cheng2011).
In addition to TFs, tomato plants also deploy gene regulation at the post-transcriptional level to preserve cellular homeostasis. MicroRNAs (miRNAs) have emerged as ubiquitous regulatory factors (Bartel, Reference Bartel2004). They are a group of small endogenous non-coding RNAs (19–24nt) that modulate post-transcriptional gene expression by directing the cleavage of target genes or by inhibiting translation depending on the extent of the complementarity between the miRNA and its target (Reinhart et al., Reference Reinhart, Weinstein, Rhoades, Bartel and Bartel2002; Chen et al., Reference Chen, Li, Lodish and Bartel2004). In addition to their role in growth and development, the link between miRNAs and responses to environmental stress is well known (May et al., Reference May, Liao, Wu, Shuai, Mccombie, Zhang and Liu2013; Agarwal et al., Reference Agarwal, Mangrauthia and Sarla2015). Even so, the investigation of high temperature-responsive miRNAs has been limited to only a few crops, including tomato (Jeong and Green, Reference Jeong and Green P2013). A recent study indicates that HS imposes a wide range of effects on miRNA expression (Zhu et al., Reference Zhu, Zhang, Tang, Qu, Duan and Jiang2019). Understanding how miRNAs are expressed and their target genes regulated in response to elevated temperature will help to better understand the molecular pathways connected to HS responses in tomato.
In Tunisia, tomato (Solanum lycopersicum L.) is one of the most important agricultural vegetable crops. Because of the warming climate, tomato growth and metabolism can be severely impaired, leading to reduced yields and fruit quality. The objective of this study is to characterize the heat-responses of commercial tomato genotypes at both the genetic and transcriptional levels. To this end, we performed Conserved DNA-derived Polymorphism (CDDP) analysis to assess genetic traits linked to HS response of tomato cultivars. CDDP analysis is a well-established PCR-based technique, used in plant genomics because of wide genome coverage. CDDPs have the advantage to be low cost, easy to apply and to give rich polymorphic patterns (Collard and Mackill, Reference Collard and Mackill2009). Although these molecular markers were applied to explore genetic diversity and features of several plants as chrysanthemum (Li et al., Reference Li, Guo, Zheng, Xia and Xian2014), chickpea (Hajibarat et al., Reference Hajibarat, Saidi, Hajibarat and Talebi2014), Rosa rugosa (Jiang and Zang, Reference Jiang and Zang2018) and Elaeagnus macrophylla (Wang et al., Reference Wang, Ma, Jia, Wu, Zang and Yu2020), this is the first use of CDDP markers to screen tomato genome resources. DNA markers display broad applications for enhancing plant adaptation to abiotic stress such as genetic identification of donors in breeding strategies and identification of candidate genes that shape responses to the stress. In our work, we used CDDP markers targeting ERF, WRKY and MYB sites. Results allowed us to group tomato cultivars into two discriminating clusters. To validate genetic results, we subjected tomato cultivars to a progressive HS treatment (0, 2, 8, 24 h). Visual phenotyping combined with genetic polymorphism data permitted the identification of tolerant and sensitive genotypes. Then, we performed qRT-PCR to investigate changes in the transcript profiles of genes that participate in the HS response. HSF family genes (HsfA1, HsfA2, HsfB1) were expressed differentially in different tomato clusters, and were likely responsible for altered HSP gene expression (HSP101, HSP17, HSP90) between cultivars. As HsfA2 tightly regulates the APX gene expression family, accumulation of transcripts for the antioxidant enzymes APX1 and APX2 were also explored. To further increase our knowledge, we included proteins encoded by CDDP markers target-genes into a functional network with related putative heat-responsive genes. We identified three gene families corresponding to LACCASES, GRAS and SPL. These genes were integrated, with their corresponding regulating miRNA, into functional modules. Thus, we analysed miRNA397-LACCASE7, miRNA171d-GRAS24 and miRNA156d-SPL2 expression in tomato genotypes to determine whether these modules are HS and genotype-responsive. Our work provides new insights into the heat responses of tomato genotypes. In addition, CDDP-based traits and molecular data are integrated to highlight tomato genetic potential, supporting future introgression strategies to overcome HS.
Materials and methods
Plant material, genomic DNA extraction and CDDP amplification
In this study, eight tomato varieties, commonly commercialized in Tunisia, were used. They correspond to Heinz, Pomodoro, Romelia, Tomate Royale, Merveille des Marchés, Chebli, Wiem and Coeur de bœuf.
Total DNA was extracted from 0.2 g of mature leaves according to the DNAeasy Plant Mini Kit (Qiagen). DNA quantification was performed with ND-1000 spectrophotometer (Nanodrop Technologies, Wilmington, USA). Nine CDDP primers, designed by Collard and Mackill (Reference Collard and Mackill2009), based on conserved sequences of well-characterized genes from diverse plant species, were used (online Supplementary Table S1). Amplifications were optimized in a final volume of 25 μl according to the following: 25 ng of sample DNA, 2.5 μM of each set of primers, 10 mM of dNTP, 3 mM of MgCl2, 0.5 U Taq polymerase (5 U/μl, Promega) with 2.5 μl of Taq polymerase buffer (10×). PCR reactions were performed in a thermal gradient cycle TProfessional TRIO Thermocycler (Biometra, Jena, Germany), under the following cycle conditions: 3 min at 94°C followed by 35 cycles of 30 s at 94°C, then 1 min at the appropriate annealing temperature for each primer and 2 min at 72°C; the final extension was held for 5 min. For each CDDP marker, PCR-amplified bands were identified and scored by the Experion™ System (Automated Electrophoresis Station).
HS treatment and phenotype scoring
Tomato seeds were sterilized with 20% bleach [v/v] and two drops of Tween for 30 min and washed with sterilized distilled water for at least three times in a row. Sterilized seeds were germinated in Petri plates with autoclaved filters. The germinated seeds were transferred to pots with peat and grown in an environmentally controlled chamber at 25°C/20°C, day/night and a 16 h light/8 h dark cycle. Tomatoes at the stage of five leaves were exposed to HS treatment for 24 h at 44°C and then allowed to recover at 25°C. The experiment was carried out with three replicates. Biological control corresponds to a set of three plants, for each variety, grown in normal conditions. Leaves were harvested at 0, 2, 8, 24 h of HS treatment and immediately frozen in liquid nitrogen. Samples were stored at −80°C until use.
Three weeks later, heat-stressed tomato plants were evaluated, and their phenotype was recorded regarding control plants.
Total RNA extraction and qRT-PCR
Total RNAs were extracted from tomato leaf tissues using TRIZOL Reagent (Trizol RNA stabilization solution, Invitrogen; Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions.
First-strand cDNA was synthesized from 2 μg of total RNA with oligo (dT) and MMLV reverse transcriptase (200 U/μl, Invitrogen) according to the manufacturer's instructions. An ABI A Prism 7000 sequence detection system (Applied Biosystems, Beverly, USA) was used for quantitative real-time PCR (qPCR). Reactions were performed in a final volume of 25 μl under the following thermal profile: 50°C for 2 min, 95°C for 2 min, followed by 39 cycles, each consisting of 95°C for 15 s and 60°C for 1 min, followed by melting at a temperature between 65 and 95°C with an increment of 0.5°C for 10 s. Tomato β-actin gene was used as an internal reference gene (Løvdal and Lillo, Reference Lovdal and Lillo2009). Genes and their corresponding primers are listed (online Supplementary Table S2). Reactions were carried out in 96-well optical reaction plates (Applied Biosystems, Beverly, USA). The reaction mixture included 2 μl of 20-fold dilution of cDNA, each primer at 2.5 mM concentration and 12.5 μl of IGreen PCR master Mix-Rox (BIOMATIK, Wilmington, USA). Each qPCR assay was run in three technical replicates and three biological replicates. Relative quantification was performed by applying the comparative 2−ΔΔCt method (Livak and Schmittgen, Reference Livak and Schmittgen2001) to calculate the changes in gene transcript as a relative fold difference between an experiment and a calibrator sample.
Prediction, validation and qRT-PCR amplification of candidate genes
The functional interacting network of proteins encoded by CDDP markers-target genes (ERF, WRKY and MYB) was built according to STRING software with the confidence parameter set at 0.4 threshold (Franceschini et al., Reference Franceschini, Szklarczyk, Frankild, Kuhn, Simonovic, Roth, Lin, Minguez, Bork, von Mering and Jensen2013). Putative genes (SPL2, GRAS24 and LACCASE7) encoding proteins involved in this network were identified, retrieved from Kegg database (https://www.genome.jp/kegg/kegg2.html) and selected for further investigation. Primers for qRT-PCR reactions were designed from the retrieved sequence of each putative gene using Primer3 designing software (Premier Biosoft, USA) (online Supplementary Table S3). For each HS stage and for controls, reactions were performed using biological triplicate replicates and duplicate technical replicates.
miRNA selection and qRt-PCR amplification
A web-based plant small RNA target analysis server ie.psRNA Target http://plantgrn.noble.org/psRNATarget/ was used to identify miRNAs, regulating candidate genes that were previously selected. Cytoscapev3.0.2 (Dai et al., Reference Dai, Zhuang and Zhao2010) was applied to visualize the interactions between miRNAs and their corresponding target genes. Effective prediction of miRNA/target duplexes was based on minimum free energies (MFEs) (Rehmsheier et al., Reference Rehmsheier, Steffen, Chsmann and Giegerich2004). miRNAs that bind nearly perfectly to their target gene were selected.
Total RNAs were extracted following the protocol described above. miRNAs were purified using NucleoSpin® miRNA–Separation of small and large RNA (Mechery-Nagel GmbH & Co. KG). First-strand cDNA synthesis and qPCR were performed as described above. Tomato β-actin gene was also used as the endogenous control (Løvdal and Lillo, Reference Lovdal and Lillo2009). Primer sequences of three miRNA (miR171d, miR397 and miR156d) were designed based on miRNA sequences retrieved from miRNA Database (http://www.mirbase.org) (online Supplementary Table S4). Primer efficiency was calculated by a standard curve of a 1:5 serial dilution of a pool of cDNA of the respective miRNA. For each HS stage and for controls, reactions were performed using biological triplicate replicates and duplicate technical replicates.
Statistical analysis
CDDP markers analysis
CDDP markers were scored for the presence (1) or the absence (0) of the same band corresponding to each set of primers. The binary 1/0 matrix generated accordingly was involved to determine the following parameters: total number of amplified bands (TNB) and number of polymorphic bands (NPB). Resolving power (Rp) was calculated according to the equation described by Gilbert et al. (Reference Gilbert, Lewis, Wilkinson and Caligari1999). Polymorphism information content (PIC), per each marker, was estimated using the formula given by Botstein et al. (Reference Botstein, White, Skolnick and Davis1980). Principal component analysis (PCA) was performed using Xlstat Software (Addinsoft, 2020), to assess the partitioning of the genotypic variability at the population level.
QRT-PCR analysis
ΔCt data analyses were performed with Data AssistTM v3.0 Software (Applied Biosystems). The output data were analysed using two-way ANOVA with times and varieties as the two predictor variables. Differences at Tukey's test HSD (P=0.05) were considered statistically significant. Analyses were performed using GraphPad Software (version 8.0, CA, USA). A heat map was performed to set up the correlation of candidate genes expression, during HS. Only the comparisons with P < 0.05 were regarded as showing differential expression.
Results
Genetic clustering of tomato cultivars
Out of nine CDDP primers used across the genomes of eight commercial tomato cultivars, two were monomorphic and thus excluded from subsequent analysis. The seven remaining markers generated a total of 46 polymorphic bands out of 53 identified (Table 1).
Table 1. Genetic parameters inferred by CDDP primers
TNB, total number of amplified bands; NPB, number of polymorphic bands; PPB, percentage of polymorphic bands; Rp, resolving power; PIC, polymorphic information content.
The number of polymorphic bands ranged from 10 (WRKYR3 and ERF2) to 4 (WRKY R2) with an average of 7.57 bands per primers. The percentage of polymorphisms (PPB) reached 100% for ERF2 and WRKY (R2B, R3B) and 85.71 for MYB2 and 80% for WRKY R3. On the opposite, PPB was lowest at MYB1 (60%) and WRKY R2 (50%), with an average value of 82.24% indicating the effectiveness of this approach in revealing the molecular polymorphism in tomato. PIC values ranged from 0.71 to 0.87, with an average value of 0.76 per primer, whereas Rp values ranged from 1 to 6.15 with an average value of 4.096 with a total of 28.675 for all primers.
PCA was performed based on CDDP markers. The first axis explained 35.52% of the total variance whereas the second axis explained 21.57% of the total variance. Projection and dispersal of tomato cultivars revealed a scattered distribution. Even so, we worth noticed that Wiem cultivar was well differentiated from the rest of tomato cultivars forming thus a putative cluster. The remaining cultivars were assigned to a single cluster (Fig. 1).
Fig. 1. Two-dimensional plot of principal component analysis (PCA) based on CDDP data.
Phenotyping of tomato plants under HS
Based on CDDP marker data, two contrasting genotypes corresponding to the Wiem and Chebli cultivars were selected. Afterwards, these two cultivars were subjected to HS treatment and screened for HS tolerance. Compared to the plant control (online Supplementary Fig. S1(a)), stressed plants displayed a wide range of symptoms. In addition, phenotyping differences validated their genetic contrasting features. Chebli cultivar showed swollen leaf blades and intervein browning leaves associated with wilting and dwarfing (online Supplementary Fig. S1(b)). Consequently, it was considered as a heat-sensitive genotype. By contrast, Wiem cultivar seems to be more tolerant with slightly curling and yellowing leaves and a wilting phenotype. Interestingly, newly emerging leaves were green, indicating that the whole plant was recovering gradually from HS imposition (online Supplementary Fig. S1(c)). Combined with genetic traits, plant behaviour under HS validates the separation of Wiem and Chebli cultivars into heat-tolerant and heat-sensitive genotypes, respectively.
Molecular heat-response factors
Genetic and phenotypic analysis allowed tomato cultivars to be assigned into two contrasting heat-sensitive and heat-tolerant groups. To support such clustering, we explored well-known heat-responsive factors expression at the transcriptional level within tomato cultivars belonging to both groups. Transcript levels of HSF family genes showed significant upregulation patterns upon HS regardless of the genotype, except for HsfA2, which was drastically repressed within Chebli genotype. HsfA2 and HsfB1 transcripts displayed a significant increase within the tolerant Wiem genotype (Fig. 2(b) and (c)).
Fig. 2. Relative expression of HsfA1 (a), HsfA2 (b), HsfB1 (c), Hsp90 (d), Hsp17(e), Hsp101(f), APX1 (g) and APX2 (h) in response to HS (44°C) in different tomato genotypes. Total RNA was purified from tomato tissues of plants exposed to HS during three stages (0, 2, 8 and 24 h). Error bars show the standard error between three replicates performed. Bars are significantly different at P < 0.05 according to Tukey's test.
In response to HS, HSPs showed distinct changes in expression levels. HSP90 transcripts accumulated gradually in both tomato genotypes at 24 h of HS stage (Fig. 2(d)). HSP17 was downregulated in Chebli genotype (Fig. 2(e)). Meanwhile, accumulation of HSP101 transcripts occurred mainly at the later stages of HS within Wiem genotype while remaining repressed in the Chebli genotype (Fig. 2(f)). Both APX1 and APX2 transcript levels gradually increased during all stages of the HS, being particularly high within the tolerant Wiem genotype (Fig. 2(g) and (h)).
miRNA-target gene modules
Putative HS-responsive genes and corresponding miRNA
To give insights into crosstalk between proteins encoded by CDDP-markers genes, an interacting protein network was built based on regulating pathways experimentally determined or retrieved from data bases (Fig. 3). According to the predicted network, MYB and ERF TFs are co-expressed with LACCASE family genes. WRKY TFs seem to participate in several functional pathways by interacting with both GRAS family and SPL TFs. GRAS genes are reported to be key actors to regulate plant responses facing abiotic stresses (Galovic et al., Reference Galovic, Orlovic and Fladung2015; Huang et al., Reference Huang, Zhiqiang, Xia, Ning and Zhengguo2015). Interestingly, multiple genes encoding SPL proteins were implicated in abiotic stress, especially salinity constraint (De Paola et al., Reference De Paola, Zuluaga and Sonnante2016). Through this functional network, WRKY proteins seem to be involved in another pattern involving MAPK5 (mitogen-activated protein kinase), a positive regulator of plant defence response, which is co-expressed with LACCASE TFs group.
Fig. 3. Functional interaction network of WRKY, ERF and MYB proteins in tomato plants.
Based on putative proteins (GRAS, SPL and LACCASE), predicted to act through this functional network, corresponding gene sequences were retrieved from the Kegg databases (https://www.genome.jp/kegg/kegg2.html). In addition, large sequence databases allowed the prediction of potential miRNA targeting these genes. Afterwards, the selection of the appropriate target was achieved through the analysis of miRNA/target complex MFEs. Three miR156d, miR397, miR171d were selected as they displayed good matches with their corresponding SPL2, LACCASE7 and GRAS24 genes (Table 2). We assumed that energy <−25 Kcal is an indicator of the strength of hybridization. Indeed, energetically best hits were recorded for miR171d/GRAS24 and miR397/LACCASE7 with −42.2 and −40 Kcal/mol, respectively.
Table 2. Selected miRNAs and their targets
Module expression
At the last stage of the HS, SPL2 gene expression was drastically downregulated in the heat-tolerant genotype, while transcripts accumulated in the sensitive one. LACCASE7 gene expression was variably enhanced during HS treatment. Transcripts were upregulated gradually in the Chebli-sensitive cultivar, particularly in the latest HS stage (Fig. 4(b)). Besides, GRAS gene was highly expressed in both heat-stressed genotypes, although preferentially in Wiem-tolerant cultivar (Fig. 4(b)). According to their gene target patterns, miRNA expression displayed opposite trends. Exposure to HS increased the abundance of miR156d and miR397 transcripts in the Wiem genotype at the 24 h HS-stage, while miR171d behaved oppositely. miR171d exhibited an expression pattern, which significantly increased at 24 h of HS within the Chebli cultivar (Fig. 4(a)). Overall, miRNA transcripts were not expressed within non-stressed tomato plants, regardless of their genotypes. Despite their differential expression patterns, the dynamic changes in miRNAs accumulation deeply outlined that HS impacts their biogenesis.
Fig. 4. Heat maps of the relative expression of three putative miRNAs (a) and gene targets (b) miR156d, miR397, miR171d, SPL2, LAC7 and GRAS24 under HS (44°C). Transcription levels were determined by quantitative RT-PCR. Total RNA was purified from tomato leaves. Colours scaled from red to purple to indicate low to higher gene expression, respectively in the genotypes Wiem and Chebli.
Discussion
Genotypic and phenotypic assortment
Candidate genes involved in abiotic stress response were targeted using CDDP DNA-based markers. Several studies have demonstrated the efficiency of the CDDP technique in assessing genetic diversity (Collard and Mackill, Reference Collard and Mackill2009; Hajibarat et al., Reference Hajibarat, Saidi, Hajibarat and Talebi2014; Wang et al., Reference Wang, Ma, Jia, Wu, Zang and Yu2020). Out of nine sets of markers, seven were polymorphic and gave allele-rich patterns. We used commercial tomato cultivars, often characterized by a narrow genetic variability. The PIC value ranges from 0.713 to 0.817 with a mean value of 0.735. This is higher than the PIC value observed in our previous study (Gharsallah et al., Reference Gharsallah, Fakhfakh, Grubb and Gorsane2016). This could be explained by the use of 19 polymorphic co-dominant SSR markers and 20 tomato cultivars, which gave a PIC value between 0.82 and 0.22 with a mean value of 0.43 and was an indicator of low genetic diversity.
The approach reported here provided additional data since it allowed cultivars to be assigned to two contrasting clusters. Indeed, PCA analysis supported a first cluster consisting of the Wiem cultivar and a second one, regrouping the remaining cultivars. Based on genetic data inferred from CDDP, two contrasting genotypes corresponding to Wiem and Chebli were selected. Upon HS imposition, both cultivars behaved differently regarding their symptoms, validating the contrasting genetic features. The Wiem cultivar displayed normal growth, with a wilting phenotype. The oldest leaves exhibited yellowing and curling whereas newly emerging leaves seemed to be recovering from stress. Conversely, the Chebli cultivar showed much more pronounced symptoms, such as intervein browning in leaves and associated wilting and dwarfing. Considering the visual appearance and differential behaviour towards HS, the Wiem genotype is considered as a heat-tolerant genotype and the Chebli cultivar is a sensitive one (online Supplementary Fig. S1).
Molecular assessment for HS tolerance
The HSF gene family corresponds to highly conserved plant-specific TFs that play an important role in plant HS responses, with at least 17 members in tomato (Koskull-Döring et al., Reference Koskull-Döring, Scharf and Nover2007). Among these HSFs, three act as central regulators of the HS response in tomato. Indeed, HsfA1, HsfA2 and HsfB1 are part of a regulatory network that controls HS-induced gene responses (Mishra et al., Reference Mishra, Tripp, Winkelhaus, Tschiersch, Theresa, Nover and Scharf2002; Schramm et al., Reference Schramm, Ganguli, Kiehlmann, Englich, Walch and Koskull-do2006). HsfA1 regulates the initial response to HS while HsfA2 controls gained thermotolerance. HsfB1 is a transcriptional repressor and a coactivator of HsfA1 as well (Fragkostefanakis et al., Reference Fragkostefanakis, Simm, El-Shershaby, Hu, Bublak, Mesihovic, Darm, Mishra, Tschiersh, Theres, Scharf, Scheleiff and Scharf2019).
At the beginning of the current experiment, HsfA1 was not expressed in any non-stressed tomato cultivars. This is in agreement with reports that HsfA1 is repressed by interaction with HSP70, under normal conditions (Hahn et al., Reference Hahn, Bublak, Schleiff and Scharf2011; Liu et al., Reference Liu, Liao and Charng2011). By contrast, at the end stage of the HS treatment, the expression of HsfA1 increased significantly in both genotypes. HsfA1 is known to greatly contribute to HS tolerance by inducing heat shock-related pathways in plants. It is constitutively expressed, acting as a master regulator for the synthesis of HsfA2/B1 and ensures the synthesis of other HSPs including HSP101 (Mishra et al., Reference Mishra, Tripp, Winkelhaus, Tschiersch, Theresa, Nover and Scharf2002). Our findings outlined the enhanced expression of HsfA2, particularly in the thermotolerant genotype. Analysing HsfA2 transcript levels during HS indicates that its expression undergoes HSP network regulation. HsfA2 transcripts were reported to be detectable even 48 h after the stress and HsfA2 protein was stable during this period as well (Fragkostefanakis et al., Reference Fragkostefanakis, Mesihovic, Simm, Paupière, Hu and Paul2016). Although it has been argued that HsfA2 is functionally equivalent to HsfA1 in experimental situations, it cannot substitute for HsfA1 because its expression in whole plants is dependent on HsfA1 (Baniwal et al., Reference Baniwal, Bharti, Chan K, Fauth, Ganguli, Kotak, Mishra S, Nover, Port, Scharf, Tripp, Weber, Zielinski and Von Koskull-Döring2004). As for HsfA2, HsfB1 is HS-inducible, although these proteins are probably not indispensable for thermotolerance (Mishra et al., Reference Mishra, Tripp, Winkelhaus, Tschiersch, Theresa, Nover and Scharf2002). In our work, both tomato cultivars showed an increase in HsfB1 expression, with transcripts accumulating the most at the end of the HS treatment. This result is in line with previous studies reporting that HsfB1 transcription is triggered in response to heat, by acting as a co-activator of HSFs in tomato (Bharti et al., Reference Bharti, Koskull-Doring, Bharti, Kumar, Tintschl-Korbitzer, Treuter and Nover2004). Deviations of HsfB1 levels permit its dual role to regulate the balance between growth and HS adaptation (Fragkostefanakis et al., Reference Fragkostefanakis, Simm, El-Shershaby, Hu, Bublak, Mesihovic, Darm, Mishra, Tschiersh, Theres, Scharf, Scheleiff and Scharf2019). The activity of HsfA1 is regulated in different ways. In normal conditions, HSP70/90 interacts with HsfA1, to repress its function (Liu et al., Reference Liu, Liao and Charng2011). Upon HS, HsfA1 dissociates from HSPs and becomes activated (Scharf et al., Reference Scharf, Berberich, Ebersberge and Nover2012; Ohama et al., Reference Ohama, Sato, Shinozaki and Yamaguchi2017). In our work, prior to HS, tomato plants exhibited a basal level of HSP90 gene expression. Then, the release of HsfA1 was concomitant with the increase in HSP90 expression (24 h). Considering the tolerant genotype, the expression profiles of Hsp101 and HSP17 were either enhanced or maintained, respectively, during all HS stages. Without regard to their order of appearance, HSP90 and HSP101 proteins correspond to a set of upregulated tomato genes in response to HS (Arce et al., Reference Arce, Spetale, Krsticevic, Cacchiarelli, Las Rivas, De Ponce, Pratta and Tapia2018; Moshe et al., Reference Moshe, Rena, Yule and Henry2016). Otherwise, HSP17 and HSP101 are mainly expressed at the latest stage of HS-treatment (Fragkostefanakis et al., Reference Fragkostefanakis, Mesihovic, Simm, Paupière, Hu and Paul2016).
On the other hand, the on-going accumulation and maintenance of HsfA2 protein is concomitant with the increased expression of HSPs (Giorno et al., Reference Giorno, Wolters-Arts, Grillo, Scharf, Vriezen and Mariani2010). This argued that HsfA2 acts as a major coactivator of HS thermotolerance (Fragkostefanakis et al., Reference Fragkostefanakis, Mesihovic, Simm, Paupière, Hu and Paul2016). Liberation of HsfA2 requires the ATP-dependent HSP70 and/or HSP101 chaperone mechanisms (Baniwal et al., Reference Baniwal, Bharti, Chan K, Fauth, Ganguli, Kotak, Mishra S, Nover, Port, Scharf, Tripp, Weber, Zielinski and Von Koskull-Döring2004). Moreover, HsfA2 is involved in an epigenetic regulatory network that is a feature of HS memory, enabling plants to survive upon recurring high temperatures (Ding et al., Reference Ding, Fromm and Avramova2012; Lämke et al., Reference Lämke, Brzezinka, Altmann and Baurle2015).
ROS, an adverse effect of HS, accumulate in chloroplasts and act as a signal transducer. To alleviate oxidative damage, HSFs behave as molecular sensors in detecting ROS molecules and subsequently, activate ROS-scavenging genes coding antioxidant enzymes as APX family (Driedonks et al., Reference Driedonks, Xu, Peters, Park and Rieu2015). In this current work, the pattern of APX1 and APX2 expression appeared to follow a similar model since transcripts accumulated at the end of HS treatment. The highest accumulation was observed in the Wiem genotype which indicates that this cultivar might enhance the activity of protective enzymes under HS. It was reported that oxidative stress occurs in parallel with high temperature, suggesting a cross-talk between heat and oxidative stress signalling (Suzuki et al., Reference Suzuki, Miller, Sejima, Harper and Mittler2012). Hence, thermotolerance is accompanied by APX upregulation and increases plant ability to counter oxidative damage. It was assumed that ROS display a dual role during HS by engendering side effects and acting as a signal transduction molecule for HSF-induced gene responses (Suzuki and Mittler, Reference Suzuki and Mittler2006).
Heat-responsive module patterns
Regarding CDDP-site targets, a predicted functional protein network was built. We focused on related proteins and their corresponding putative miRNA. Based on databases, we retrieved gene sequences encoding proteins predicted by the network. Then, we selected miRNAs that bind perfectly to target genes according to the online software (Rehmsheier et al., Reference Rehmsheier, Steffen, Chsmann and Giegerich2004). miRNAs take part in plant responses to abiotic stresses, across a complex target-gene network regulation. According to the complementarity between miRNA and their corresponding target, miRNAs regulate gene expression through mRNAs cleavage or protein translation inhibition (Winter and Diederichs, Reference Winter and Diederichs2011; Iwakawa and Tomari, Reference Iwakawa and Tomari2013). In wild tomato, heat-responsive miRNAs were reported to be either up or downregulated (Zhou et al., Reference Zhou, Wang, Jiang, Cao, Sun, Liu and Wu2016). Authors speculated that miRNA could be down or upregulated depending on the plant species and the type of stress.
Three modules consisting of miRNA-target genes were selected to be explored in two contrasting tomato groups, under HS imposition. Opposite changes during HS stages outlined an inverse relationship between miRNA expression level and their corresponding target genes.
The first module consists of SPL2-miR156d. SPLs family members are reported to be targeted by miR156 since they present a putative-specific binding site (Salinas et al., Reference Salinas, Xing, Höhmann, Berndtgen and Huijser2012). miR156, the first miRNA identified in plants, is reported to display a pleiotropic function in plant biological mechanisms and to be highly expressed following various environmental stresses (Axtell and Bowman, Reference Axtell and Bowman2008). Transgenic plants over expressing miR156 exhibited enhanced and stable tolerance to HS, suggesting that miR156 was required for plant HS memory (Stief et al., Reference Stief, Altmann, Hoffmann, Pant, Scheible and Baurle2014). miR156 was reported to be heat-induced in many plants including Triticum aestivum L. and Brassica rapa L. (Xin et al., Reference Xin, Yu, Xie, Peng, Ni and Sun2010; Yu et al., Reference Yu, Wang, Lu, de Ruiter, Cariaso, Prins, van Tunen and He2012) as well as tomato (Zhou et al., Reference Zhou, Wang, Jiang, Cao, Sun, Liu and Wu2016). It is interesting to note that high expression of miR156 in plants leads to decreased leaf size and reduced apical dominance (Cao et al., Reference Cao, Li, Wang, Nan, Wang, Lu, Jiang, Li, Shi, Fang, Yuan, Zhao, Li, Liu and Kong2015). This is likely a strategy undertaken by plants to escape the adverse effects of high temperatures. In our experiments, upregulation of miR156d was noticeable at 24 h of HS treatment, only within the tolerant Wiem genotype. miR156d displayed a good match with its target SPL2. Consequently, our data showed that SPL2 gene exhibited low expression at 24 h of the HS stage, in agreement with miR156d abundance. SPL2 expression was enhanced within the sensitive Chebli genotype, at the end of the HS treatment. Interestingly, downregulation of SPL2 under HS is linked to HS memory (Stief et al., Reference Stief, Altmann, Hoffmann, Pant, Scheible and Baurle2014; Zhu et al., Reference Zhu, Zhang, Tang, Qu, Duan and Jiang2019). miR156 is also reported to maintain the expression of HS memory-related genes including HSPs and HsfA2 (Ohama et al., Reference Ohama, Sato, Shinozaki and Yamaguchi2017). Our results showed that HsfA2 was activated concomitantly with miR156 accumulation in the tolerant genotype, validating thus the involvement of this module to establish a durable HS tolerance.
The second module, consisting of LACCASE7-miR397, showed a dynamic expression pattern with an enhanced LACCASE7 expression during all stages of HS for both cultivars. Nevertheless, the Wiem genotype maintained moderate expression at the end of the HS stage while Chebli showed a remarkable accumulation of LACCASE7 transcripts. While LACCASE7 was gradually expressed during HS stages, miR397 accumulation was negligible within the sensitive genotype. miR397 is a conserved miRNA across plant species (Jones-Rhoades and Bartel, Reference Jones-rhoades and Bartel2004) and its expression varied depending on the considered stress (Li et al., Reference Li, Fu, Chen, Huang and Yu2011, Reference Li, Hu, Chen, Wang, Xiong and Zheng2017; Xia et al., Reference Xia, Zhao, Li, Chen, Jiao, Wu, Zhou, Yu and Fan2018). miR397 displays a perfect match with the fifth exon of LACCASE7 gene, contributing to its downregulation (Abdel-Ghany and Pilon, Reference Abdel-ghany and Pilon2008; Pan et al., Reference Pan, Zhao, Yu, Bai and Dong2017). LACCASE gene repression might limit some non-essential biological processes, allowing plants to conserve energy to cope with the adverse effect of HS and to be protected from oxidative damages (Abdel-Ghany and Pilon, Reference Abdel-ghany and Pilon2008). This may explain the intervein-leaf burning observed with heat-stressed Chebli genotype which displays a high expression of LACCASE7 gene. Our findings also underline that LACCASE7 expression was slightly enhanced in the Wiem cultivar at the late HS stage. This may indicate that tolerant genotype undergone a physiological pathway to preserve essential elements for physiological functions.
The third module consists of GRAS24-miR171d. GRAS family consists of proteins that regulate the features of growth, development and responses to biotic and abiotic environments (Huang et al., Reference Huang, Shiyuan, Zhiqiang, Dongbo, Guojian, Lu, Ren and Li2017). In tomato, GRAS24 has a conserved MIR-binding sequence that perfectly matches with miR171 (Huang et al., Reference Huang, Zhiqiang, Xia, Ning and Zhengguo2015). miR171 is known to play a critical role in abiotic stress within several plant species (Liu et al., Reference Liu, Tian, Li, Wu and Zheng2008; Kong et al., Reference Kong, Elling, Chen and Deng2010; Hwang et al., Reference Hwang, Shin and Yu2011).
In our work, miR171d displayed a perfect match with its target GRAS24, and opposite changes in transcript accumulation. Within both cultivars, GRAS24 exhibited an increased expression during all the stages of HS with a notable high expression at 24 h for the tolerant cultivar. Conversely, miR171d displayed highest expression within the sensitive genotype.
In summary, elevated temperatures activate stress-related candidate genes to induce a cascade of transcriptional activation and ROS-scavenging that contributes to optimal defence, normal growth and fitness. Integrating genetic and phenotyping data allowed the discrimination of two contrasting genotypes. The Wiem genotype displayed heat-thermotolerance requiring a coalition of several pathways that culminate in the (i) upregulation of HSFs with subsequent HSP activation, (ii) increase of HsfA2 accumulation for HS memory, and (iii) upregulation of APX1/2 ROS-scavenging enzymes to overcome HS-induced oxidative stress. In parallel, regulation of LACCASE7-miR397 and SPL2-miR156d modules, upon HS, seems to be significantly modulated in a genotype-specific manner, leading to better adaptation to such a stress. The Wiem genotype is therefore a valuable cultivar that could be used as a donor in breeding programmes to enhance tomato heat tolerance.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S1479262121000083
Acknowledgments
The authors thank Dr Bootheina Majoul, Associate Professor of English Literature & Studies (ISLT) and Dr Benjamin Field, University of Aix-Marseille for their conscientious reading and meticulous corrections of the manuscript.
