Expression Proﬁling of Heat Shock Protein Genes as Putative Early Heat-Responsive Members in Lettuce

: High temperatures due to global warming can cause harmful effects on the productivity of lettuce, a cool-season crop. To identify lettuce heat shock protein ( HSP ) genes that could be involved in early responses to heat stress in plants, we compared RNA transcriptomes between lettuce plants with and without heat treatment of 37 ◦ C for 1 h. Using transcriptome sequencing analyses, a total of 7986 differentially expressed genes (DEGs) were identiﬁed including the top ﬁve, LsHSP70A , LsHSP70B , LsHSP17.3A , LsHSP17.9A and LsHSP17.9B, which were the most highly differentially expressed genes. In order to investigate the temporal expression patterns of 24 lettuce HSP genes with a fold-change greater than 100 under heat stress, the expression levels of the genes were measured by qRT-PCR at 0, 1, 4, 8, 14, and 24 h time points after heat treatment. The 24 LsHSP genes were classiﬁed into three groups based on the phylogenetic analysis and/or major domains available in each protein, and we provided a potential link between the phylogenetic relationships and expression patterns of the LsHSP genes. Our results showed putative early heat-responsive lettuce HSP genes that could be possible candidates as breeding guides for the development of heat-tolerant lettuce cultivars.


Introduction
Lettuce (Lactuca sativa L.) is one of the most economically important leafy vegetables in the world and its production is currently increasing due to the growing interest in healthy food [1]. Lettuce is a cool-season crop, with an optimal growing temperature range of 17 to 28 • C, and floral initiation generally occurs between 21 and 27 • C [2]. Heat stress refers to a temperature environment exceeding the optimal limit for a certain species. When exposed to heat stress, physiological disorders such as tip burn, premature bolting, puffy heads and rib discoloration are commonly observed in lettuce [3].
Climate change is expected to increase the frequency and severity of extremes in temperature, resulting in significant damage to vegetable crop production [4]. According to the global climate report issued by the National Oceanic and Atmospheric Administration, USA. (NOAA), 2020 was ranked as the second warmest year in 141 years of recording, with the global land and ocean surface temperature increasing by an average of +0.98 • C [5]. In addition to the continuous increase in average temperature due to the greenhouse effect, short periods of abnormally high temperatures, resulting in a significant reduction in crop growth and yield are also a result of climate change [6]. Peng et al. (2004) found a decline of approximately 10% in yield for each 1 • C increase during the dry season [7]. Even a few hours of temperatures that are higher than optimal can induce dormancy in lettuce seeds [8]. Spinach and lettuce may produce stalks (bolt) and cucumber tends to produce a higher ratio of male to female flowers in response to high temperature, which can reduce yield and commercial value [9].
Unlike animals, plants are unable to change their position when directly exposed to persistently changing environmental stresses [10]. Thus, plants cannot avoid the exposure to heat stress, which poses a serious threat to plants [11]; however, plants have a sophisticated adaptive system at the cellular and molecular levels [12,13]. Heat shock proteins (HSPs), also called stress proteins, play a central role in the protection of cells exposed to different types of stress in plant tissue [14]. In plants, these proteins are grouped into five principal classes based on their approximate molecular weight: (1) Hsp100 family, (2) Hsp90 family, (3) Hsp70 family, (4) Hsp60 family, and (5) small heat shock proteins (sHsps) family [15]. Under heat stress, HSPs, as molecular chaperones, bind to heatdenatured proteins and mediate refolding, assembly for repairing denatured proteins or degradation of misfolded proteins to maintain homeostasis of proteins [16][17][18].
It is important to develop heat tolerant lettuce varieties because there is a very low consumer acceptance threshold for lettuce with symptoms typically associated with exposure to heat stress. So far, studies on lettuce transcriptomics under cold or diverse light conditions are limited [19,20] since the genetically validated reference genome of lettuce was first reported in 2017 [21]. There is a need to identify the genes that respond to heat stress and investigate their expression patterns during the stress. In this study, we focused on identifying putative early heat stress-responsive HSP genes in lettuce by analyzing transcriptomes and their temporal expression patterns during heat stress and the phylogenetic relationships among them. This information will be helpful in developing the next generation of lettuce cultivars with enhanced heat tolerance through molecular breeding.

Plant Material and Growth Conditions
Korean blue leaf lettuce (Lactuca sativa cv. LB Sang 169) seeds obtained from the National Institute of Horticultural and Herbal Science (NIHHS), Korea were sown in 50 cell trays and grown in the NIHHS greenhouse, which was maintained at an optimal temperature of 22/20 • C (day/night). At 18 days after sowing, seedlings were transferred into a temperature-controlled LED growth chamber (VS-91G09M-HL. VISION Co., Daejeon, Korea). Seedlings were grown under white light at an intensity of 120 umol m−2 s−1 and a 16/8 h day/night photoperiod, with the temperature maintained at 20/18 • C (day/night), and 60-65% humidity. The 25-day-old plants were exposed to a heat stress treatment of 37/35 • C (day/night), without any change in the light intensity. At different time points (0 h, 1 h, 4 h, 8 h, 14 h and 24 h after heat treatment), the aerial part of the lettuce was sampled, frozen in liquid nitrogen and stored at −80 • C. After 24 h of heat treatment, plants were transferred back to the optimal temperature for observation of phenotypic alteration (Supplementary Materials, Figure S1).

Morphometric Analysis
The maximal length (maximal length from the apical to the basal parts of the leaf) and width of the 4th leaf (measured at the midpoint) of each plant were measured at 3, 6, 10 and 13 days after heat treatment. At 13 days after heat stress treatment, the fresh weight of each plant was measured and the dry weight was determined by oven-drying it at 70 • C for 72 h.

RNA Preparation and Transcriptome Sequencing Analysis
Total RNA was extracted from the frozen leaf samples using a Plant RNeasy Extraction Kit (Qiagen, Hilden, Germany). The RNA samples were treated with DNase I (Qiagen, Germany), and cDNA was synthesized using the iScript cDNA synthesis kit (Bio-Rad, Hercules CA, USA). The integrity of the RNA was examined by agarose gel electrophoresis. RNA-Seq libraries were generated with a TruSeq Stranded Total RNA LT Sample Prep Kit (Illumina, San Diego CA, USA) and sequenced by a paired-end-sequencing method using Illumina's HiSEquation 2000 (Macrogen, Seoul, Korea).

Data Analysis Quality Control
Quality control (QC) of raw paired-end reads was done by FastQC v0.11.7. The Q20% and Q30% are the proportion of nucleotides with a quality score > 20 and 30, respectively. In order to reduce biases in the analysis, artifacts such as low-quality reads, adaptor sequence, contaminant DNA or polymerase chain reaction (PCR) duplicates were trimmed by using Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic, accessed on 29 January 2021). Trimmed reads were mapped to reference genome (Lsat_Salinas_v7) with HISAT2 version 2.1.0, Bowtie2 2.3.4.1 (https://ccb.jhu.edu/software/hisat2/index.shtml, accessed on 29 January 2021). After mapping reads, StringTie version 2.1.3b (https:// ccb.jhu.edu/software/stringtie/, accessed on 29 January 2021) was used for transcript assembly and the expression profiles of the assembled transcript for each sample were calculated. The FPKM (fragment per kilobase of transcript per million mapped reads) value is used as the expression profile in the case of paired-end sequencing.

Differentially Expressed Genes (DEG)
DEG analysis was performed to compare samples grown under high temperature (37/35 • C day/night) and optimal temperature (20/18 • C day/night, control) conditions. Statistical analysis was conducted using the fold change (FC) per comparison pair. The significant results were selected based on conditions of |FC| ≥ 2 and p-value < 0.05.

Quantitative RT-PCR (qRT-PCR) Analysis
The synthesized cDNA was used for qRT-PCR analysis by using AccuPower 2X GreenStar qPCR Master Mix (Bioneer, Daejeon, Korea) and a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Relative mRNA levels were determined by normalizing the PCR threshold cycle number of each target gene with that of the reference gene, UBQ21 [22]. In the qRT-PCR analysis, three technical repeats were measured for each biological replicate analyzed. The primers used for qRT-PCR analyses are presented in Supplementary Materials, Table S1.

Phylogenetic Analysis
Protein sequences of lettuce HSPs were aligned using ClustalW and the phylogenetic tree was generated based on the neighbor-joining method implemented in MEGA X [23] under the parameters of the Jones-Taylor-Thornton (JTT) model with uniform rates among sites, and pairwise deletion of gaps.

Effects of Heat Stress on the Early Vegetative Growth of Lettuce
To examine the effect of heat stress on the early vegetative growth of lettuce, lettuce plants were exposed to high temperature conditions (37/35 • C day/night) for 24 h. The morphological traits of the lettuce were investigated after 3, 6, 10 and 13 days of the heat treatment compared to those grown under the optimal temperature condition (20/18 • C day/night, control). Heat-induced inhibition of leaf growth was observed in the lettuce with the heat stress compared to control ( Figure 1A,B). Comparison of leaf size between plants with and without heat treatment exhibited a significant difference at the 13-day time point after heat treatment: the length and width of the leaves in the plants with heat treatment decreased by 15% and 14%, respectively, compared to control. Furthermore, both fresh and dry weights of lettuce with heat treatment were also significantly reduced by 19% and 27%, respectively, compared to control ( Figure 1C). There were no signs of physiological disorders associated with high temperature stress such as tip burn, premature bolting, or rib discoloration observed in our experiment with a 24 h heat treatment. physiological disorders associated with high temperature stress such as tip burn, premature bolting, or rib discoloration observed in our experiment with a 24 h heat treatment. Asterisks indicate significant differences compared with the control (* p < 0.05, ** p < 0.01, Student's t test). Each bar represents mean ± SE of three independent experiments.

Heat Treatment Induced the Expression of a Heat Shock Protein Gene, LsHSP70-2711 in Lettuce
In order to investigate the expression levels of a known lettuce heat shock protein gene, LsHSP70-2711 (COGE: Lsat_1_v5_gn_9_22920.1) during heat stress [24], the abundance of the transcript was determined by qRT-PCR analyses in lettuce at different time points of heat treatment (1, 4, 8, 14 and 24 h). A significant increase in LsHSP70-2711 expression was observed in lettuce under heat stress condition for all time points examined, compared to the control (Supplementary Materials, Figure S2). Based on the highest expression level of LsHSP70-2711 at the 1 h time point of heat stress treatment, RNA samples prepared at the time point were applied for transcriptome sequencing to identify early heat-responsive HSP genes in lettuce. Asterisks indicate significant differences compared with the control (* p < 0.05, ** p < 0.01, Student's t test). Each bar represents mean ± SE of three independent experiments.

Heat Treatment Induced the Expression of a Heat Shock Protein Gene, LsHSP70-2711 in Lettuce
In order to investigate the expression levels of a known lettuce heat shock protein gene, LsHSP70-2711 (COGE: Lsat_1_v5_gn_9_22920.1) during heat stress [24], the abundance of the transcript was determined by qRT-PCR analyses in lettuce at different time points of heat treatment (1, 4, 8, 14 and 24 h). A significant increase in LsHSP70-2711 expression was observed in lettuce under heat stress condition for all time points examined, compared to the control (Supplementary Materials, Figure S2). Based on the highest expression level of LsHSP70-2711 at the 1 h time point of heat stress treatment, RNA samples prepared at the time point were applied for transcriptome sequencing to identify early heat-responsive HSP genes in lettuce.

Transcriptome Profiling of Lettuce with Heat Treatment
To compare the expression profiles of lettuce grown in two different conditions-with and without heat treatment-whole transcriptome sequencing was conducted using RNA isolated from each group at the 1 h time point after heat treatment. The FASTQC software was utilized for quality control of the raw sequence data. After removing low-quality reads and adaptor sequences, we generated 67,144,596 clean reads, including 6,763,559,164 total read bases from the sample under heat stress (37/35 • C day/night) for 1 h, and 75,171,018 clean reads, including 7,571,262,233 total read bases from the sample under optimal temperature (20/18 • C day/night, control). Trimmed data were mapped against the reference genome (Lsat_Salinas_v7) using the HISAT2 program. The obtained mapped reads were then assembled successfully using StringTie-e version 2.1.3b, based on the reference genome (Table 1).

Differentially Expressed Gene Analysis
Significantly up-and downregulated genes with |FC| ≥ 2 and raw p-value < 0.05 were identified (Figure 2A). A total of 7986 unigenes were discovered to be differentially expressed between samples with and without heat treatment. Compared to the sample with optimal temperature, 4402 unigenes were upregulated and 3584 were downregulated in the sample with heat treatment (Figure 2A).
The smear plot in Figure 2B represents the distribution of DEGs in the up-or downregulated section and confirms the genes are differentially expressed compared to the control according to expression volume, which is defined as the geometric mean of two groups' expression level. Twenty-four LsHSP genes upregulated more than a 100-fold under the 1 h heat treatment were selected for further studies (Table 2). Among them, the five HSP genes with the highest differential expression were LsHSP70A, LsHSP70B, LsHSP17.3A, LsHSP17.9A and LsHSP17.9B.    The two dashed blue lines represent two-fold changes. This plot is colored such that those points having a fold-change less than |FC|< 2 are shown in gray, |FC|≥ 2 in pale blue, and |FC|≥ 2, raw. p-value < 0.05 in dark blue.

Experimental Validation of Gene Expression by qRT-PCR
To further validate the comparative transcriptome results for the 24 HSP genes, transcript level variances between the heat-stressed and optimal conditions were examined using qRT-PCR analyses. The transcript levels of 24 genes were highly upregulated in the lettuce treated with 1 h heat stress, which coincided with the gene expression detected by the transcriptome analyses ( Figure 3).  . The x-axis represents temperature and the y-axis represents relative transcript level (×Log 2 fold change). Each bar represents mean ± SE of three independent experiments.

Time-Course Gene Expression Patterns of the Lettuce HSPs and Their Phylogenetic Relationships
To investigate changes in the gene expression under heat stress, the expression levels of the 24 LsHSP genes in plants exposed to heat stress for 0 h, 1 h, 4 h, 8 h, 14 h and 24 h were determined by qRT-PCR. The highest expression level of the 24 LsHSP genes was observed at the 1 h time point after heat treatment, suggesting that all the LsHSP genes selected are likely to respond early to heat stress (Figure 4). Almost all the LsHSP genes examined had a distinct peak expression in response to heat at the 1 h time point, which drastically declined at the 4 h time point, with the exceptions of LsSHSPA, LsHSP17.5 and LsHSP18.1, which exhibited more than 50% of their peak expression at the 4 h time point.

Time-Course Gene Expression Patterns of the Lettuce HSPs and Their Phylogenetic Relationships
To investigate changes in the gene expression under heat stress, the expression levels of the 24 LsHSP genes in plants exposed to heat stress for 0 h, 1 h, 4 h, 8 h, 14 h and 24 h were determined by qRT-PCR. The highest expression level of the 24 LsHSP genes was observed at the 1 h time point after heat treatment, suggesting that all the LsHSP genes selected are likely to respond early to heat stress (Figure 4). Almost all the LsHSP genes examined had a distinct peak expression in response to heat at the 1 h time point, which drastically declined at the 4 h time point, with the exceptions of LsSHSPA, LsHSP17.5 and LsHSP18.1, which exhibited more than 50% of their peak expression at the 4 h time point.  Using the 24 protein sequences of lettuce HSPs with mRNA levels greater than 100fold increase compared to those of the control, a phylogenetic tree was constructed using MEGA X ( Figure 5A). Three major clusters were identified among the lettuce HSPs based on the available domains ( Figure 5B). The first cluster of genes contained the alpha-crystallin domain, the next cluster contained the nucleotide-binding domain (NBD) of the HSP70 superfamily, and the last cluster contained the Histidine kinase/HSP90-like ATPase domain. To visualize the expression levels of the LsHSP genes under heat stress, a heat map was also generated displaying the changes in gene expression at each time point examined under heat stress together with the available domains in the HSPs ( Figure 5B). It is noteworthy that similar expression patterns were detected among members belonging to the same group in the phylogenetic tree, which raises the possibility of a conserved functional activity.
Horticulturae 2021, 7, x FOR PEER REVIEW 10 of 15 Using the 24 protein sequences of lettuce HSPs with mRNA levels greater than 100fold increase compared to those of the control, a phylogenetic tree was constructed using MEGA X ( Figure 5A). Three major clusters were identified among the lettuce HSPs based on the available domains ( Figure 5B). The first cluster of genes contained the alpha-crystallin domain, the next cluster contained the nucleotide-binding domain (NBD) of the HSP70 superfamily, and the last cluster contained the Histidine kinase/HSP90-like ATPase domain. To visualize the expression levels of the LsHSP genes under heat stress, a heat map was also generated displaying the changes in gene expression at each time point examined under heat stress together with the available domains in the HSPs ( Figure  5B). It is noteworthy that similar expression patterns were detected among members belonging to the same group in the phylogenetic tree, which raises the possibility of a conserved functional activity.

Effects of Heat Stress Treatment on the Early Vegetative Growth of Lettuce
Heat stress is one of the major problems that limit plant growth and metabolism, leading to significant loss of yield potential for various agricultural crops around the globe [25]. Due to drastic and rapid changes in the global climate, the negative effects of these stresses on the crop have been aggravated [26,27]. Each plant species has a specific temperature range for minimum, maximum, and optimum growth and development [28]. Lettuce is a cool-season crop with an optimal growing temperature ranging from 17 to 28 °C [29]. We demonstrated that the growth of lettuce was significantly suppressed by heat

Effects of Heat Stress Treatment on the Early Vegetative Growth of Lettuce
Heat stress is one of the major problems that limit plant growth and metabolism, leading to significant loss of yield potential for various agricultural crops around the globe [25]. Due to drastic and rapid changes in the global climate, the negative effects of these stresses on the crop have been aggravated [26,27]. Each plant species has a specific temperature range for minimum, maximum, and optimum growth and development [28]. Lettuce is a cool-season crop with an optimal growing temperature ranging from 17 to 28 • C [29]. We demonstrated that the growth of lettuce was significantly suppressed by heat treatment for 24 h, indicating that a short period of abnormal high temperatures can have serious implications for lettuce yields [6].

Heat Stress Treatment Increased the Expression of LsHSP70-2711 in Lettuce
To identify lettuce HSP genes that respond early to heat stress, lettuce samples were harvested at different time points under heat stress and the expression level of LsHSP70-2711, a known lettuce HSP gene was examined [24]. The expression level of the gene was the highest at the 1 h time point after heat treatment, and then decreased gradually until the 16 h time point. At the 24 h time point, when ambient temperature increased by 2 • C, a slight increase in the gene expression was observed. A previous work supports our experimental data, where LsHSP70-2711 exhibited the highest expression levels after heat stress for 1 h, although it was not among the 24 greatest differentially expressed HSP genes we identified [24].
Heat stress can trigger several mechanisms of defense such as induced gene expression to produce stress-associated molecular chaperones called HSPs. In general, major HSPs play important roles in preventing mis-folded, denatured and aggregated proteins caused by stress conditions [30,31], which contribute to abiotic and biotic stress tolerance and facilitate recovery and survival after a return to normal growth conditions [32]. Thus, the peak of LsHSP70-2711 expression at the 1 h time point in response to heat stress supports the notion that LsHSP70 may play an important role in early responses to stress, which are associated with the accumulation of unfolded or denatured proteins.

Transcriptome Profiling of Heat-Treated Lettuce and DEG Analysis
Transcriptome analyses have recently been reported in vegetables such as tomato [33], pepper [34], radish [35], and bottle gourd [36] under heat stress conditions. However, transcriptome analyses regarding heat stress in lettuce have not yet been reported. Thus, based on the temporal gene expression pattern of LsHSP70-2711 under heat stress, we compared the expression profile of lettuce treated with and without heat treatment for 1 h through whole transcriptome sequencing. We identified the most upregulated five genes, LsHSP70A, LsHSP70B, LsHSP17.3A, LsHSP17.9A and LsHSP17.9B, which are all HSPs that were revealed by transcriptomic analyses. We also selected 24 lettuce HSP genes that showed more than a 100-fold upregulation as candidate genes for further analyses. Those HSP genes with the highest expression at the 1 h time point after heat stress may contribute to the recognition of stress and production of HSPs to prevent the irreversible aggregation of other proteins and may participate in refolding proteins during the heat stress condition [37].

Temporal Expression Patterns of Lettuce HSP Genes under Heat Stress
We performed qRT-PCR to examine the expression patterns of the identified 24 HSP genes during heat stress. The expression of the 24 genes was dramatically elevated in lettuce with heat treatment for 1 h compared to control. However, the expression of 24 HSP genes tended to be reduced as the heat stress condition continued. Thus, the level of difference in the expression of 24 LsHSP genes in lettuce plants with and without heat treatment was greatest at the 1 h time point after heat stress, but was reduced as plants adapted to the stress over time. These observations are consistent with the previous result showing that maximum expression levels of HSP18.1 and HSP17.9 occurred at the beginning of a 4 h heat stress period and declined by the end of the 4 h stress despite the continued exposure to high leaf temperature [38]. Some endogenous factors other than heat stress may also regulate the expression of HSP genes when plants are consistently exposed to heat stress.
In the current study, most genes had a dramatically decreased level of gene expression at the 4 h time point, while LsSHSPA, LsHSP17.5 and LsHSP18.1 had relatively gradual decreases in expression. These three genes retained a level of expression of over 50% at the 4 h time point compared to the peak expression at the 1 h time point under the heat stress condition. Various HSPs have been identified in many crop plants that are related to heat stress, with different expression patterns under heat stress and with overlapping and distinct functions [39]. The gene expression of 9 HSP genes in rice (Oryza sativa) had different maximum expressions at various time points depending on the gene [40], supporting our finding that some HSP genes are quickly and sharply induced in response to heat stress, while others can be induced slowly to induce adaptation to higher temperatures.
The 24 LsHSP genes were grouped based on the phylogenetic analysis and/or major domains available in each protein. Small HSPs (sHSPs) with low molecular weight of 15-42 kDa were the most dominant proteins produced in plants during heat stress [38,[41][42][43] and these proteins have a conserved alpha-crystallin domain (ACD) [44]. In our study, we found that the most dominant genes among top 24 HSP genes code for sHSPs, and particularly, the expression of LsHSP17.9 genes was highly upregulated under the heat stress conditions. The major function of sHSPs is involvement in the degradation of the misfolded proteins, mostly through the ubiquitin-proteasome pathway [45] and/or preventing irreversible unfolding or wrong protein aggregation by binding to partially folded or denatured proteins [43]. One of the characteristics that distinguishes sHSPs from other HSP class members is that they exhibit an ATP-independent chaperone activity [46]. Thus, the high expression level of sHSP genes may be primarily responsible for the regulation of the early response when plants are exposed to heat stress.
In addition to sHSP genes, HSP70 genes were also significantly induced in response to heat stress. Previous research has shown that HSP70s are a class of molecular chaperones that play a crucial role in protecting plant cells from the detrimental effects of heat stress [47]. Structurally, HSP70s contain three major functional domains: a N-terminal nucleotidebinding domain (NBD) with ATPase activity, a substrate-binding domain (SBD), and a C-terminal lid region that covers the SBD [48][49][50]. The binding of ATP at the NBD, and its subsequent hydrolysis influences the substrate binding affinity of the SBD through allostery, suppressing protein aggregation and refolding misfolded proteins [51]. HSP83 is homologous to the chaperone HSP90, which has a histidine kinase/HSP90-like ATPase domain and also plays an important role in protein folding from the early stage of protein synthesis together with HSP70 and other associated co-chaperones [52]. The expression of HSP83 has been shown to increase during heat shock in Arabidopsis thaliana [53]. We also found that the expression of LsHSP83 was enhanced under heat stress conditions. Thus, it is suggested that LsHSP70 and LsHSP83 play a critical role in protecting lettuce cells from heat stress.

Conclusions
In the present study, we identified 24 lettuce HSP genes through transcriptome sequencing that were highly upregulated within 1 h of heat stress treatment. They were divided into three clusters based on phylogenetic relationships together with the conserved domains available in each protein. Furthermore, we demonstrated potential relationships between the temporal expression patterns of LsHSP genes at different time points of heat treatment and the phylogenetic relationships. Functional validation of each LsHSP gene remains to be examined in our future work.