Transcriptome sequencing of garlic reveals key genes related to the heat stress response

With global warming, heat stress has become an important factor that seriously affects crop yield and quality. Therefore, understanding plant responses to heat stress is important for agricultural practice, but the molecular mechanism of high-temperature tolerance in garlic remains unclear. In this study, ‘Xusuan No. 6’ was used as the experimental material. After heat stress for 0 (CK), 2 and 24 h, transcriptome sequencing was used to screen metabolic pathways and differentially expressed genes (DEGs) closely related to heat stress and was further verified by quantitative real-time polymerase chain reaction (qRT-PCR). A total of 86,110 unigenes obtained from the raw transcriptome sequencing data were spliced. After 2 h of heat treatment, the expression levels of 8898 genes increased, and 3829 genes were decreased in leaves. After 24 h, the expression levels of 7167 genes were upregulated, and 3176 genes were downregulated. Gene Ontology enrichment analysis showed that DEGs were mainly enriched in seven categories: cellular processes, metabolic processes, binging, catalytic activity, cellular anatomical entity and protein-containing complex response to stimulus. Kyoto Encyclopedia of Genes and Genomes pathway enrichment showed that DEGs are involved in protein processing in the endoplasmic reticulum, plant hormone signal transduction, phenylpropanoid biosynthesis, and photosynthetic antenna proteins. Six genes were selected and further verified by qRT-PCR. In this study, the full-length transcriptome of garlic was constructed, and the regulatory genes related to the heat resistance of garlic were studied. Taken together, these findings can provide a theoretical basis for the cloning of heat resistance genes in garlic and for the analysis of heat resistance mechanisms.

Garlic (Allium sativum L.) is an important cash crop belonging to the Allium family.Moreover, garlic is a coldloving crop with poor adaptability to high temperatures, and this seriously affects the growth and development of garlic, leading to a decline in its yield and quality 1,2 .
With changes in the global climate, heat stress has become one of the major environmental stressors limiting plant metabolism and crop productivity 3,4 .Studies have shown that heat stress can negatively affect plant growth, physiological processes, and metabolism 5 .Therefore, to cope with the damage caused by high temperatures, plants have formed complex physiological and biochemical adaptation mechanisms in the long-term evolution process to clear or repair denatured proteins and other biological macro-molecules and maintain cell homeostasis 6,7 .High temperatures usually damage proteins in the endoplasmic reticulum (ER), chloroplasts and cytoplasm 8 .Endoplasmic reticulum (ER) stress is a stress response caused by the accumulation of unfolded proteins in the ER lumen at high temperatures 9 .Some important regulatory factors have been found to help translate and fold proteins in the ER and cytoplasm, thereby regulating protein homeostasis and alleviating heat damage 10 .Liu et al. identified two types of unfolded protein response (UPR) signaling pathways in plants: one uses two transmembrane alkaline bright basic-leucine zipper transcription factors (bZIP17 and bZIP28), and the other involves a double protein kinase (RNA splicing factor IRE1) and its target RNA (bZIP60) 11 .In addition, plant hormones, as signaling molecules, are widely involved in plant responses to abiotic stress 12,13 .Abscisic acid (ABA), salicylic acid (SA), auxin/indoleacetic acid (IAA), and other plant hormones are also involved in the response of plants to high-temperature stress 14 .Salicylic acid can remove hydrogen peroxide produced by cucumber under high-temperature stress, thereby improving the heat resistance of cucumber 15 .The jasmonate www.nature.com/scientificreports/signaling pathway is involved in the molecular regulation mechanism of plant response to high-temperature stress and provides an important gene resource for wheat heat tolerance molecular breeding 16 .
RNA-sequencing (RNA-seq) technology has been widely used in the study of the plant response to abiotic stress, including heat stress in maize (Zea mays L.) seedlings 17 , heat stress in eggplant leaves 18 , and drought stress in peanut seedlings 19 .Wang et al. performed transcriptome sequencing in garlic variety 'Cangshan Siliuban' and found that genes involved in hormone signaling and cell wall remodeling play an important role in garlic response to high salt stress 20 .These RNA-seq technologies have been widely used in the study of abiotic plant reactions, and many stress-response genes related to different metabolic processes have been identified, enriching information on the plant stress regulation network.However, there are few reports on the response genes and heat resistance mechanism of garlic to high-temperature stress.In this study, to analyze the molecular mechanism of response and adaptation of garlic seedlings to high-temperature stress, we exposed garlic seedlings to a high-temperature environment and conducted transcriptome sequencing to study the response genes at different time points of garlic high-temperature stress, providing a theoretical basis for the analysis of the high-temperature tolerance mechanism in garlic.

Physiological response of garlic under high temperatures
Under high-temperature stress, the leaves of garlic wilted after 24 h of heat stress (Fig. 1A).As shown in Fig. 1B, membrane peroxidation produced malondialdehyde (MDA) in large quantities, and the activity of MDA decreased gradually after 24 h but was maintained at a high level.The activity of peroxide (POD) increased rapidly after 2 h of high-temperature stress, indicating that the reactive oxygen scavenging system responded rapidly.With continuous exposure to a high temperature, the POD reached its peak at 24 h.High-temperature treatment significantly affects the physiological process of photosynthesis and ultimately leads to an increase or decrease in chlorophyll content.After high-temperature treatment, the chlorophyll content of garlic leaves gradually decreased.

Assembly and functional annotation
As shown in Fig. 2A, 86,110 unigenes were obtained, with an average length of 849 bp and 38.58% GC content.All Unigenes were annotated in the Non-Redundant Protein Database (Nr), KOG, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Swiss-Prot databases, and the annotation rates were 32.6%, 19.10%, 32.20%, and 22.49%, respectively (Fig. 2B).Transcriptome sequencing quality analysis showed that multiple high-quality base sequences were obtained in the control group (CK) and after high-temperature treatment for 2 and 24 h, with Q20 and Q30 reaching more than 97% and 93%, respectively, and GC% reaching more than 40%, indicating that the transcriptome sequencing data in this study were of high quality and could be used for further research (Table 1).

Functional classification
A total of 42,483 unigenes were grouped into 26 functional categories (Fig. 3).Among them, cluster "general function prediction only" was the largest group, followed by "posttranslational modification, protein turnover, chaperones, " "signal transduction mechanisms, " and "translation ribosomal structure and biogenesis." Through the analysis of DEGs in garlic exposed to high-temperature treatment, it was found that they were primarily concentrated in stress-related Gene Ontology (GO) annotations and KEGG pathways (Fig. 4).GO analysis of DEGs showed that the differences were mainly concentrated in biological processes, cell components and molecular functions.In the molecular function category, "binding" and "catalytic activity" were the most dominant classes, followed by "transporter activity" and "structural molecule activity.In the cellular process category, "cellular anatomical entity" and "protein-containing complex" were the most highly represented groups.After high-temperature stress treatment, metabolic processes and cellular processes play an important role in biological processes.
KEGG enrichment demonstrated the metabolic pathways and bioinformatics functions of DEGs, as shown in Fig. 5.In this study, 2777 unigenes were distributed to five groups.The largest subgroup was "Global and overview maps, " followed by "carbohydrate metabolism, " which both belonged to the metabolism category.

Transcriptional response of garlic seedlings to heat stress
Principal component analysis (PCA) was used to determine the similarity of gene expression in CK and garlic leaf samples after 2 and 24 h under high-temperature treatment (Fig. 6).Gene expression was significantly divided into three principal components in garlic seedlings under control and heat stress conditions.The first principal component accounted for 72.6% and the second for 21.6%.
The transcription data of garlic treated with high-temperature stress for 2 and 24 h were compared with those of the control, and the expression levels of each gene in the four samples were compared and filtered with |log 2 FC|≥ 1 and FDR < 0.05.There were 12,727, 10,343 and 15,401 up-or downregulated unigenes detected in the "CK-versus -T2, " "CK-versus -T24, " and "T2 versus T24" comparisons, respectively (Fig. 7A).The comparisons revealed that the three groups shared 1584 different genes (Fig. 7B).www.nature.com/scientificreports/

GO analysis of DEGs
Based on GO annotation analysis, DEGs were divided into three categories-biological process, cell component, and molecular function-with 50 different classification groups (Fig. 8).As shown in Table 2, regarding biological processes, DEGs were mainly concentrated in cellular processes, metabolic processes, biological regulation, and responses to stimuli.With respect to molecular functions, DEGs were mainly concentrated in binging, catalytic activity, and transporter activity.As cell components, DEGs were mainly concentrated in cellular anatomical entities and protein-containing complexes in response to stimuli.In addition, with the increase in high-temperature  www.nature.com/scientificreports/stress, the number of upregulated DEGS increased significantly, and the overall number of DEGS also increased.Therefore, we hypothesize that garlic may respond to high-temperature stress by affecting the expression levels of relevant DEGs at the transcriptional level.

KEGG analysis of DEGs
KEGG enrichment is used to demonstrate the metabolic pathways and bioinformatics functions of differential genes.KEGG enrichment of DEGs in this study was divided into 5 KEGG pathway branches according to their metabolic and signaling pathways, including metabolism, genetic information processing, environmental  www.nature.com/scientificreports/information processing, cellular process and organismal systems (Fig. 9).In the "CK -versus -T2" comparison, 1961 DEGs were enriched into 130 pathways, among which 1236, 475, 101, 69, and 80 DEGs were involved in metabolism, genetic information processing, environmental information processing, cell process and organic system, respectively.In the "CK-versus -T24" comparison, 2561 DEGs were enriched into 130 pathways, among which 1771, 518, 107, 101, and 64 DEGs were involved in metabolism, genetic information processing, environmental information processing, cell process and organic system pathway branches, respectively.In the "T2versus -T24" comparison, 3259 DEGs were enriched into 130 pathways, among which 2252, 660, 130, 119, and 98 DEGs were involved in metabolism, genetic information processing, environmental information processing, cell process and organic system pathway branches, respectively.In the three comparison weights, the DEGs involved in metabolic pathways were the most involved, and the DEGs involved in cellular process pathway branches were the least involved.In addition, with the extension of high-temperature stress exposure, DEGs in the organic system of the "CK-versus -T24" group decreased, and DEGs in the other four KEGG pathway branches gradually increased, especially in metabolic pathways.
In the CK versus T2 and CK versus T24 comparison groups, eight pathways were found to be the same, including protein processing in the ER (ko04141), plant hormone signal transduction (ko04075), photosynthesisantenna proteins (ko00196), riboflavin metabolism (ko00740), diterpenoid biosynthesis (ko00904), taurine and hypotaurine metabolism (ko00430), fatty acid metabolism (ko00062) and glycerolipid metabolism (ko00561) (Fig. 10).Among these, the protein processing in the ER (ko04141) and plant hormone signal transduction (ko04075) pathways were significantly enriched by DEGs.KEGG pathway enrichment analysis showed that protein processing in the ER was the most significantly enriched pathway, suggesting that it plays a central role in garlic under high-temperature stress.

Protein processing in ER-related genes under heat stress: expression change
Heat stress can cause protein folding error, ER stress and cell death.Therefore, the ER plays an important role in plant response to high-temperature stress.In this study, after high-temperature stress, 155 and 125 DEGs were found after 2 and 24 h of high-temperature treatment in the protein processing in the ER (ko04141).Most of them were associated with heat stress, including HSP40, HSP70 and HS90 (Fig. 11, Table S4).In the CK-versus -T2 and CK -versus -T24 comparison groups, 90 and 71 heat shock protein (HSP) genes, respectively, were identified among the DEGs.Moreover, the expression levels of most HSP genes in CK and at 2 and 48 h showed a trend of first increasing and then decreasing (Fig. 12).These results suggest that garlic may rapidly respond to heat stress by regulating the expression of HSP and ER-related genes at the transcriptional level to clear misfolded proteins induced by heat stress.

Plant hormone signal transduction-related genes under changes in heat stress expression
In the plant hormone signal transduction (ko04075) pathway, 63 and 66 DEGs were found after 2 and 24 h of high-temperature treatment, respectively.There were 20 genes affecting plant growth, including ABA, AUX, IAA, SAUR, and ARF, and most of these genes were downregulated in the late stage of heat stress.Among them, four ABF genes (Unigene0090235, Unigene0030557, Unigene0060159, and Unigene0087890) were upregulated after exposure to high temperatures for 2 and 24 h.Furthermore, there were nine genes affecting stomatal closure.The expression of PYL genes was downregulated after 24 h of treatment, and the expression of PP2C genes was upregulated after 2 and 24 h of treatment (Table 3).At the same time, under high-temperature stress, species DEGs related to plant signal transduction pathways were abundant, indicating that plant hormones play an important role in regulation under high-temperature stress (Fig. 13).

Quantitative real-time PCR (qRT-PCR) validation of differential gene expression under heat stress
We used RT-qPCR to further verify the accuracy of the transcriptome data.Six unigenes were selected from the list of DEGs based on their potential functions: Unigene0083975 (HSP90-3), Unigene0007832 (NDPK1), Uni-gene0046359 (WRKY24), Unigene0055905 (CAT1), Unigene0001589 (AUX22D), and Unigene0049563 (MYB5).The analysis showed that although there were certain differences in the multiples of zero upregulation or downregulation of expression detected by RNA-seq and qRT-PCR, as shown in Fig. 14, the expression trend reflected by the results of qRT-PCR was consistent with the results of transcriptome sequencing, which may be caused by the different detection ranges and expression calculation methods of the two methods.The transcriptome sequencing data of garlic high-temperature stress were accurate and reliable.

Models of the molecular mechanisms of garlic in response to heat stress
Based on analysis of the transcriptome sequencing results, we developed a schematic model (Fig. 15).In the protein processing in the ER pathway, the expression of most HSPs is rapidly induced by high temperatures and may play an important role in preventing the formation of misfolded protein structures in garlic cells, as well as in rescuing aggregated or denatured proteins.In the plant hormone signal transduction pathway, ABA is an important hormone for regulating stress responses.After ABA bound to receptor protein PYR/PYL, the inhibition of kinase SnRK2 activity by phosphatase PP2C was relieved, and SnRK2 activity was activated to induce the plant stress response.

Discussion
As an important environmental factor, temperature is involved in regulating many aspects of plant growth and development and has a very important impact on agricultural production 21,22 .Heat stress can lead to stalk elongation, early flowering and a serious decrease in crop yield 23 .Garlic leaves gradually withered and turned yellow after 24 h of exposure to high-temperature stress.POD activity increases rapidly to clear active oxygen species, and MDA accumulates continuously to destroy the cell membrane structure.Therefore, it is of great significance to study the mechanism of garlic coping with high-temperature stress for the cultivation of hightemperature-resistant varieties.
Studies have shown that plant responses to heat stress are carried out through a complex gene regulatory network 24 .In this study, RNA-seq was used to analyze the changes in the gene expression profiles of garlic treated at different times under high temperatures.In total, 12,727 and 10,343 DEGs were detected in response to heat stress at 2 and 24 h, respectively.Through the GO function of differential genes and the enrichment of the KEGG metabolic pathway, the primary metabolic processes, such as amino acid metabolism, carbohydrate metabolism and some secondary metabolic processes, began to be affected after 2 h of high-temperature stress.At the same time, the high-temperature stress response genes involved in the metabolic pathway of plant hormone signal transduction began to accumulate in the early 2 h of treatment.The results indicated that the stress signal reception and signal transduction processes of garlic were activated in the early stage of high-temperature stress and also began to regulate some primary and secondary metabolic processes.Under abiotic stress conditions, the accumulation of unfolded or misfolded proteins in the ER can lead to ER stres 25.Li et al. also found that the upregulation of ER protein processing in maize was most significant after hightemperature stress, which was consistent with our research results 26 .In this study, KEGG enrichment analysis showed that the ER protein processing pathway was most enriched, indicating that ER stress occurs and triggers the UPR, which may play a central role in garlic leaves under heat stress.The ER luminal binding protein (BIP) is the core molecular chaperone that assists protein folding in the ER; upregulation of BIP gene expression can promote the folding of ER protein 27,28 .In protein processing in the ER pathway, ER stress mediated by inositoldependent enzyme 1 (IRE1) is the driver of the early (1-3 h) heat stress response 29 .Under high-temperature stress, the IRE1A gene of garlic also significantly increased after 2 h of treatment, indicating that the IRE1A www.nature.com/scientificreports/gene played a certain role in the early heat shock response of garlic.High-temperature stress can cause protein degeneration in organisms 30 .Heat shock protein (HSP) acts as a molecular chaperone to assist protein refolding, stabilization, intracellular transport, and degradation, prevent the accumulation of damaged proteins, and maintain the stability of the intracellular environment 31,32 .High-temperature stress transcripome analysis of quinoa showed that HSF transcription factors such as LOC110702486, LOC110697083, LOC110729384, LOC110709409, and LOC11073445 were the core regulatory factors of transcription under high temperatures 33 .When plants are subjected to high-temperature stress, HSP accumulates a large amount of expression in a short time and participates in the stress response, but with the increase in treatment time, the stress response is weakened and HSP is degraded 5 .Transcriptome analysis showed that HSP genes were mainly concentrated in the protein processing pathways in the ER.Moreover, the expression levels of most HSP genes increased first and then decreased after 2 and 24 h of high-temperature treatment, which may be due to the short-term nature of HSP.
Plant endogenous hormones are very important regulatory factors in plant growth and development 34 .When plants are exposed to environmental stress, the synthesis, distribution and transport of endogenous hormones change significantly to activate the plant stress resistance mechanism 35 .In response to high-temperature stress, plant hormones, including ABA, IAA, AUX, SAUR, ARF, and PYL, are signaling compounds that regulate important aspects of growth, development, and environmental stress response 5 .Plants subjected to high-temperature stress can rapidly accumulate a large amount of ABA, which can bind to PYR/PYL (Pyrabactin resistance 1/ Pyrabactin resistance 1-like) to inhibit PP2C (Type 2C protein phosphatase) activity, but snPK2 (Sucrosenonfermenting 1-related protein 2) remains active and phosphorylates the downstream transcription factor ABF 36 .For example, Kim et al. found that overexpression of ABA transcription factor ABF3 can make plants resistant to high-temperature stress 37 .In this study, garlic may also form PYR-PP2C-SnRK2 signal transduction complex under high-temperature stress and then participate in the activation of ABA to initiate the mechanism of resistance to high temperature.The signal transduction pathway of auxin in response to high-temperature stress mainly includes transport inhibitor response1/auxin signaling F-box (TIR1/AFB), auxin/indoleacetic acids proteins (Aux/IAA), auxin response factors (ARFs) and other protein components 38 .Chen et al. found that the expression of many SbARFs in sorghum was upregulated by high-temperature stress, and SbARF17/24 accumulated in vascular tissue under high-temperature stress, indicating that SbARF may participate in the high-temperature response 39 .In this study, the expression levels of four ABF genes (Unigene0090235, Unigene0030557, Uni-gene0060159, and Unigene0087890) were upregulated at 2 and 24 h under heat stress.These DEGs may be key candidate genes related to auxin in garlic.At the same time, DEGs in plant signal transduction pathways were abundant under high-temperature stress.It shows that plant hormones play an important role in the regulation of high-temperature stress.
In this study, we screened and obtained important DEGs related to protein processing in the ER.At present, there are few transcriptome studies on abiotic stress in garlic seedlings.Therefore, a comparative analysis of transcriptomes based on different garlic varieties and growth stages may be of great significance in elucidating some common molecular mechanisms of garlic abiotic stress response.

Plant materials and growth conditions
The garlic cultivar 'Xusuan No. 6' was used as the experimental material, and it was conserved at the Xuzhou Institute of Agricultural Sciences in Jiangsu Xuhuai Area.It was selected and the bulbs that had been released from dormancy full and undamaged were selected and then grown in a mixture of organic soil and vermiculite in a greenhouse (34° 27′ N, 117° 29′ E).The temperature of the culture in the artificial climate chamber was set to 25 °C, with a relative humidity of 30% and a photoperiod of 16 h of light and 8 h of dark.To ensure that the plants are not exposed to drought at high temperatures, after two weeks, the garlic seedlings were transplanted to a hydroponic tank with 1/2 Hoagland nutrient solution and maintained in this solution for six days.In addition, plants with consistent and robust growth were selected for high-temperature treatment (38 °C).The relative humidity of the air was set to 30%, and each treatment was repeated 3 times.At 0 (CK), 2 and 24 h of treatment, the upper leaves of the garlic plant were removed, frozen in liquid nitrogen and stored at − 80 °C for later use.

Determination index
The activities of POD and MDA were determined according to the manufacturer's instructions (Beijing Solarbio Technology Co., Ltd., Beijing, China).The chlorophyll content was determined using the ethanol extraction colorimetric method.In the presence of a small amount of quartz sand and calcium carbonate powder, leaf chlorophyll was extracted with 95% ethanol.The entire extraction process was carried out under dark conditions.The absorbance values at 665, 649, and 470 nm were determined using a spectrophotometer, with 95% ethanol serving as the control.The chlorophyll contents were then calculated based on these absorbance values.Each index was determined 3 times, and the average value was calculated using Excel 2016 software.The method was analyzed and compared with Duncan's method using SPSS software.

RNA extraction, library construction, and RNA-seq
Total RNA was extracted using a plant RNA extraction kit (Tiangen, Beijing, China) in accordance with the manufacturer's instructions.The integrity of nucleic acid samples was tested by agarose gel electrophoresis.The purity of nucleic acid is determined by detecting OD value of nucleic acid by NanoDrop.The Agilent 2100 assay

Figure 1 .
Figure 1.Phenotypic and physiological analysis of changes in garlic under heat stress.(A) Plant phenotype of garlic under heat temperature for 0 (CK) and 24 h; (B) Measurement of physiological indicators; different letters indicate significant difference between treatments at the 5% probability level, the same as below.

Figure 2 .
Figure 2. (A) Length distribution of unigene in the garlic transcriptome.(B) Venn diagram of the number of unigenes annotated in NR, KEGG, KOG, and Swiss-Prot database.

Figure 3 .
Figure 3. KOG function classification of unigenes in the garlic transcriptome.

Figure 4 .
Figure 4. GO analysis of unigenes identified in the garlic transcriptome.

Figure 5 .
Figure 5. KEGG classification of unigenes in the garlic transcriptome.

Figure 6 .
Figure 6.Principle component analysis plot of samples of garlic seedlings grown under normal conditions or exposed to heat stress.CK represents the control seedlings, while T2 and T24 represent heat-stressed seedlings.

Figure 9 .
Figure 9. KEGG function analysis of differentially expressed genes.

Figure 10 .
Figure 10.Top 20 KEGG pathways in garlic leaves under heat stress.The first and outer circles show the top 20 KEGG pathways that were enriched, while the scale outside the circle indicates the number of genes.Different colors represent different ontologies.Next, the KEGG pathway number in the background gene, along with the Q value, can be seen in the middle circles.Dark colors indicate genes that are upregulated, and light colors indicate genes that are downregulated.Below is a display of the specific value.The inner circle: rich factor).

Figure 11 .
Figure 11.Pathway of processing pathways in the ER.

Figure 14 .
Figure 14.qRT-PCR validation of the expression levels of DEGs.

Figure 15 .
Figure 15.Schematic model of ABA signaling pathways and ER in the heat stress response of garlic.

Table 1 .
Base information statistics.

Table 2 .
GO enrichment analysis of the DEGs.

Table 3 .
Changes in DEGs in the KEGG pathway of plant hormone signal transduction at different time points during heat stress.