Transcriptomic Changes Induced by Drought Stress in Hardneck Garlic during the Bolting/Bulbing Stage

Induced by Drought Stress in Hardneck Garlic during the Bolting/Bulbing Abstract: Garlic ( Allium sativum L.) is an economically important, monocotyledonous plant with a strong taste and odor. Drought stress adversely affects its growth, development, and yield, particularly during the bolting/bulbing stage. Herein we performed RNA-seq to assess transcriptomic changes induced by drought stress in bolting/bulbing hardneck garlic plants (Purple Glazer). We observed that drought stress signiﬁcantly reduced photosynthesis rate, fresh weight, and leaf water content. Transcriptomic analysis of garlic leaves under normal conditions and drought stress led to the identiﬁcation of 5215 differentially expressed genes (2748 up- and 2467 downregulated). The upregulated DEGs were primarily involved in “biological process”, “metabolic process”, “oxidation-reduction process”, carbohydrate and lipid metabolism, and “proteolysis”, whereas the downregulated DEGs were mainly involved in “biological process” and metabolism of various molecules. In addition, genes encoding abscisic acid biosynthetic and catabolic enzymes, heat shock proteins, and E3 ubiquitin ligases were signiﬁcantly altered by drought stress, indicating involvement in drought tolerance. A further comparison with the DEGs related to salinity stress-treated garlic revealed 867 and 305 DEGs with a similar and reverse expression alteration tendency, respectively.


Introduction
Due to climate change and global warming, drought has become a significant adverse environmental factor, posing a major threat to agriculture. The development of droughttolerant crops is crucial to world food security in the long term. Therefore, a comprehensive understanding of plant responses to drought stress is a prerequisite to develop droughttolerant plants, particularly crops.
Plants display an array of physiological, biochemical, and molecular responses to survive drought stress. To reduce water loss, the closing of stomata and reduction in photosynthesis in leaves occur under water-limiting conditions [1]. This process is mainly mediated by abscisic acid (ABA), a plant hormone, which is derived from carotenoids, a group of natural pigments in plastids [2]. ABA plays a central role in drought tolerance by regulating stomatal closure, root growth and architecture, and transcription of stressresponsive genes [3]. Genes involved in ABA biosynthesis and signal transduction are evidently regulated by drought stress [4][5][6].
The generation of reactive oxygen species (ROS) in plants is another typical response to drought stress. Excessive ROS subsequently causes damage to lipids, proteins, and DNA molecules. On the other hand, plants synthesize antioxidant enzymes and other antioxidants to scavenge ROS [7]. Besides, abiotic stresses, including drought, lead to the accumulation of denatured or dysfunctional proteins, which are harmful to plant growth and development [8]. Heat shock proteins (HSPs) are a group of chaperones associated

Drought Stress Treatment and Physiological Measurements
Approximately 5-month-old garlic plants were transplanted into a greenhouse and allowed to grow for another 2 weeks at 25 °C and 14-h light/10-h dark cycle. The plants that started to bolt were exposed to drought stress (1 June 2020-8 June 2020) by withholding water, while well-watered plants served as control. Subsequently, we measured the photosynthesis rate of control and drought-treated plants using LI-COR 6400XT (LI-COR, Lincoln, NE), according to manufacturer instructions, followed by shoot height and shoot fresh weight (FW) measurement. Shoot dry weight (DW) was measured after shoots were dried at 80 °C for 3 days. Leaf water content was calculated using the formula [(FW − DW)/FW] × 100. All measurements were obtained from at least three plants. Statistical analysis was performed using Student's t-test. Some leaves were collected in liquid nitrogen and stored at −80 °C for RNA extraction.

RNA Extraction and Library Construction
Total RNA was extracted from garlic leaf tissues using an E.Z.N.A. ® Plant RNA Kit (Omega Bio-Tek Company, Norcross, GA, USA), according to manufacturer instructions. RNA quality and quantity were first estimated using a NanoDrop™ One/OneC Microvolume UV-Vis Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and

RNA Extraction and Library Construction
Total RNA was extracted from garlic leaf tissues using an E.Z.N.A. ® Plant RNA Kit (Omega Bio-Tek Company, Norcross, GA, USA), according to manufacturer instructions. RNA quality and quantity were first estimated using a NanoDrop™ One/OneC Microvolume UV-Vis Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and then on an Agilent 2100 Bioanalyzer System (Santa Clara, CA, USA). Only RNA samples with RNA integrity number of >8.0 were used for library construction. RNA-seq libraries of drought-treated and control garlic leaves were constructed using 4 µg total RNA with the TruSeq Stranded mRNA Sample Preparation Kit (Illumina Company, San Diego, CA, USA). Adaptor sequences are listed in Supplementary Table S1. The libraries were sequenced on an Illumina NovaSeq 6000 sequencing platform to generate 2 × 150 bp paired-end reads.

RNA-Seq Data Analysis and Functional Annotation
Raw RNA-seq reads were first processed by trimming the adaptor and low quality sequences using Trimmomatic [16]. The reads were then aligned to the rRNA database using Bowtie [17] allowing up to three mismatches [18]. The resultant rRNA-mapped reads were removed, and clean RNA-seq reads were mapped to the garlic reference genome using STAR [19,20]. Following alignments, raw counts for each gene were calculated and then normalized to fragments per kilobase of exon model per million mapped fragments (FPKM). We herein used FPKM values to measure gene expression levels; therefore, the reported LFCs are fold changes of mean FPKM values of drought-treated samples to those of control samples. The raw p-values were adjusted for multiple tests using the false discovery rate [21]. DEGs were identified using DESeq2 [22], with the criteria being log2(fold change) ≥1 or ≤−1 and adjusted p ≤ 0.05. We annotated the predicted protein-coding genes in the garlic genome by comparing their protein sequences to the GenBank non-redundant (NR) database using BLASTp (E-value cutoff = 1 × 10 −5 ). We also compared protein sequences to the InterPro database using InterProScan [23]. GO annotations were assigned using Blast2GO [24] based on BLASTp results from the NR database and the output from InterProScan. Enriched GO terms were identified using the Blast2GO suite. Pathway annotation was performed with BlastKOALA [25], and significantly enriched pathways were identified using the hypergeometric test.

Drought Stress Significantly Affected the Growth of Garlic Plants
To examine the impact of drought stress on garlic growth, garlic plants that started to bolt were transplanted into a greenhouse and subjected to drought stress by withholding water. On exposure to drought stress, the plants began to wilt; senescence occurred in the tips of leaves. Seven days post-treatment, the plants were clearly affected by drought conditions, showing wilted brown leaves. Under control, well-watered conditions, the plants showed normal growth ( Figure 1B). At the end of the treatment, the photosynthesis rate, shoot height, shoot fresh weight, and leaf water content were 0.43 µmol m −2 s −1 , 36.2 cm, 18.4 g, and 81.2%, respectively, whereas under normal conditions, the values were 6.72 µmol m −2 s −1 , 52.6 cm, 35.4 g, and 86.1% (Figure 2), respectively, indicating that drought stress detrimentally influenced plant metabolism and growth.

Functional Categorization of DEGs
GO annotation analysis of DEGs revealed that the upregulated DEGs were primarily involved in "biological process", "metabolic process", and "oxidation-reduction process" ( Figure 4A). Besides, pathways associated with carbohydrate and lipid metabolism were upregulated, such as "carbohydrate metabolic process", "cellular carbohydrate catabolic process", "oligosaccharide metabolic process", "monosaccharide metabolic process", "lipid metabolic process", and "cellular lipid metabolic process", indicating that drought stress altered carbohydrate and lipid metabolism in garlic leaves ( Figure 4A, Supplementary Table S5). The removal of denatured or oxidized proteins is important for plants to survive drought stress, and drought has been found to increase proteases and peptidases involved in this process [26]. In garlic leaves, "proteolysis" (75 genes) was enriched among the upregulated DEGs ( Figure 4A).  All measurements were obtained from at least three plants. Significance (*) was determined (t-test); p < 0.05 indicated statistical significance.

RNA-Seq Analysis under Normal Conditions and Drought Stress
Six libraries for RNA-seq were prepared from three replicates of leaves collected under normal conditions and drought stress (Supplementary Table S1). A total of 424,788,158 raw reads and 388,865,998 clean reads were accordingly generated. Approximately 87.05-94.83% reads from each sample were mapped to the assembled garlic genome [19] ( Supplementary  Table S2), of which around 83.47% were aligned to the regions of the predicted proteincoding genes, while the remaining 16.53% were aligned to the potential non-coding regions (Supplementary Table S3).
Using the criteria adjusted p ≤ 0.05 and fold change ≥ 2.0, 5215 DEGs were identified, among which, 2748 genes were upregulated and 2467 were downregulated ( Figure 3, Supplementary Table S4). The downregulated DEGs were mainly involved in "biological process" and metabolism of various molecules, such as "organic substance metabolic process", "primary metabolic process", "nitrogen compound metabolic process", "macromolecule metabolic process", and "organonitrogen compound metabolic process", suggesting a decrease in these pathways within plant cells ( Figure 4B). Moreover, the "photosynthesis", "photosynthesis, light harvesting", and "chlorophyll metabolic process" pathways were downregulated ( Figure 4B), which was consistent with the decrease in the photosynthesis rate in garlic leaves upon exposure to drought stress ( Figure 2A).

Functional Categorization of DEGs
GO annotation analysis of DEGs revealed that the upregulated DEGs were primarily involved in "biological process", "metabolic process", and "oxidation-reduction process" ( Figure 4A). Besides, pathways associated with carbohydrate and lipid metabolism were upregulated, such as "carbohydrate metabolic process", "cellular carbohydrate catabolic process", "oligosaccharide metabolic process", "monosaccharide metabolic process", "lipid metabolic process", and "cellular lipid metabolic process", indicating that drought stress altered carbohydrate and lipid metabolism in garlic leaves ( Figure 4A, Supplementary Table S5). The removal of denatured or oxidized proteins is important for plants to survive drought stress, and drought has been found to increase proteases and peptidases involved in this process [26]. In garlic leaves, "proteolysis" (75 genes) was enriched among the upregulated DEGs ( Figure 4A).

Altered Pathways in Garlic Leaves by Drought
KEGG pathway enrichment analysis showed that the pathways associated with metabolism (secondary metabolites, starch and sucrose, phenylpropanoid, glycerolipid, and so on) and signal transduction (plant hormone, estrogen, MAPK signaling pathways, and plant-pathogen interaction) were up-regulated in the leaves after drought treatment (Figure 5A, Supplemental Table S6). On the other hand, the "carbon fixation in photosynthetic organisms", "photosynthesis", "photosynthesis-antenna proteins", "porphyrin and chlorophyll metabolism" were downregulated, in line with the reduced photosynthesis rate in leaves under drought conditions ( Figure 5B). The other downregulated pathways included "Biosynthesis of amino acids", "Plant hormone signal transduction", "Glyoxylate and dicarboxylate metabolism", "Ribosome biogenesis in eukaryotes", and others ( Figure  5B, Supplemental Table S6). The downregulated DEGs were mainly involved in "biological process" and metabolism of various molecules, such as "organic substance metabolic process", "primary metabolic process", "nitrogen compound metabolic process", "macromolecule metabolic process", and "organonitrogen compound metabolic process", suggesting a decrease in these pathways within plant cells ( Figure 4B). Moreover, the "photosynthesis", "photosynthesis, light harvesting", and "chlorophyll metabolic process" pathways were downregulated ( Figure 4B), which was consistent with the decrease in the photosynthesis rate in garlic leaves upon exposure to drought stress (Figure 2A).

Altered Pathways in Garlic Leaves by Drought
KEGG pathway enrichment analysis showed that the pathways associated with metabolism (secondary metabolites, starch and sucrose, phenylpropanoid, glycerolipid, and so on) and signal transduction (plant hormone, estrogen, MAPK signaling pathways, and plant-pathogen interaction) were up-regulated in the leaves after drought treatment ( Figure 5A, Supplemental Table S6). On the other hand, the "carbon fixation in photosynthetic organisms", "photosynthesis", "photosynthesis-antenna proteins", "porphyrin and chlorophyll metabolism" were downregulated, in line with the reduced photosynthesis rate in leaves under drought conditions ( Figure 5B). The other downregulated pathways included "Biosynthesis of amino acids", "Plant hormone signal transduction", "Glyoxylate and dicarboxylate metabolism", "Ribosome biogenesis in eukaryotes", and others ( Figure 5B, Supplemental Table S6).

HSPs and E3 Ubiquitin Ligases were Significantly Upregulated
HSPs function as molecular chaperones, mediating protein folding, assembly, trans location, and degradation under normal and stress conditions [32]. HSPs have been re ported to be commonly induced by drought and other stresses, probably enhancing plan tolerance by re-establishing normal protein conformation and protein homeostasis [33] In drought-treated garlic leaves, 124 genes encoding HSPs were significantly upregulated including 62 class I HSPs, 22 70-kDA HSPs, 14 class II HSP-like, and other HSPs (Figure 7 and Supplementary Table S4). E3 ubiquitin ligase is the major component of the ubiquitin mediated protein degradation pathway and determines substrate specificity [34]. The ex pression of 24 E3 ubiquitin ligase genes was induced by drought in garlic leaves (Figure 7 and Supplementary Table S4), indicating that ubiquitin-mediated protein degradation was significantly enhanced in response to drought.

Commonly and Reversely Regulated-DEGs in Garlic Leaves on Exposure to Drought and Salinity Stresses
Drought and salinity are two major abiotic stresses that impact agriculture. In a sep arate study, 3616 DEGs (1571 up-and 2045 downregulated) were identified by RNA-seq in garlic leaves exposed to salinity stress (GenBank accession no. PRJNA682570). We herein simply compared DEGs to identify common and unique genes altered by drough and salinity stresses. Only 9.3% (256 out of 2748 genes) genes upregulated under drough conditions were also upregulated in response to salinity stress, indicating that drough

HSPs and E3 Ubiquitin Ligases Were Significantly Upregulated
HSPs function as molecular chaperones, mediating protein folding, assembly, translocation, and degradation under normal and stress conditions [32]. HSPs have been reported to be commonly induced by drought and other stresses, probably enhancing plant tolerance by re-establishing normal protein conformation and protein homeostasis [33]. In drought-treated garlic leaves, 124 genes encoding HSPs were significantly upregulated, including 62 class I HSPs, 22 70-kDA HSPs, 14 class II HSP-like, and other HSPs (Figure 7 and Supplementary Table S4). E3 ubiquitin ligase is the major component of the ubiquitinmediated protein degradation pathway and determines substrate specificity [34]. The expression of 24 E3 ubiquitin ligase genes was induced by drought in garlic leaves (Figure 7 and Supplementary Table S4), indicating that ubiquitin-mediated protein degradation was significantly enhanced in response to drought. secondary metabolites", "starch and sucrose metabolism", "RNA degradation", "galactose metabolism", and other pathways were upregulated under salinity stress but downregulated under drought conditions (Figure 8). On the contrary, 169 genes involved in "metabolic pathways", "biosynthesis of secondary metabolites", "tyrosine metabolism", "glycerolipid metabolism", "glycine, serine and threonine metabolism", and other pathways were downregulated by salinity stress but upregulated under drought conditions (Figure 8).

Commonly and Reversely Regulated-DEGs in Garlic Leaves on Exposure to Drought and Salinity Stresses
Drought and salinity are two major abiotic stresses that impact agriculture. In a separate study, 3616 DEGs (1571 up-and 2045 downregulated) were identified by RNA-seq in garlic leaves exposed to salinity stress (GenBank accession no. PRJNA682570). We herein simply compared DEGs to identify common and unique genes altered by drought and salinity stresses. Only 9.3% (256 out of 2748 genes) genes upregulated under drought conditions were also upregulated in response to salinity stress, indicating that drought and salinity treatment activated distinct responses in garlic leaves. These genes were associated with, for example, "metabolic pathways", "biosynthesis of secondary metabolites", "protein processing in endoplasmic reticulum", and "MAPK signaling pathway" (Figure 8). On the other hand, 24.8% (611 out of 2467) genes downregulated under drought conditions were also suppressed by salinity. These genes were associated with, for example, "metabolic pathways", " biosynthesis of secondary metabolites", carbon metabolism, "photosynthesis", "glyoxylate and dicarboxylate metabolism", and "plant hormone signal transduction" (Figure 8), suggesting that both drought and salinity negatively affected the metabolism of various molecules, photosynthesis, as well as signal transduction. Interestingly, 136 genes involved in "ribosome biogenesis in eukaryotes", "biosynthesis of secondary metabolites", "starch and sucrose metabolism", "RNA degradation", "galactose metabolism", and other pathways were upregulated under salinity stress but downregulated under drought conditions (Figure 8). On the contrary, 169 genes involved in "metabolic pathways", "biosynthesis of secondary metabolites", "tyrosine metabolism", "glycerolipid metabolism", "glycine, serine and threonine metabolism", and other pathways were downregulated by salinity stress but upregulated under drought conditions (Figure 8).

Discussion
In the case of hardneck garlic accessions, scapes and bulbs initially form around the same time in May or June. This developmental stage has been reported to be affected by various environmental factors, such as temperature and photoperiod [35], and in garlic, it is the most sensitive stage to drought stress [15]. In this study, we performed RNA-seq to examine transcriptomic changes in garlic plants exposed to drought stress during the bolting/bulbing stage. Mapping the sequencing reads to the garlic reference genome and statistical analysis led to the identification of 5215 DEGs, providing insights into the molecular mechanisms involved in drought tolerance. Differentially expressed genes (DEGs) identified upon exposure to salinity and drought stresses were analyzed using a Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/). The Kyoto Encyclopedia of Genes and Genomes pathways that DEGs were assigned to are shown in tables.

Discussion
In the case of hardneck garlic accessions, scapes and bulbs initially form around the same time in May or June. This developmental stage has been reported to be affected by various environmental factors, such as temperature and photoperiod [35], and in garlic, it is the most sensitive stage to drought stress [15]. In this study, we performed RNA-seq to examine transcriptomic changes in garlic plants exposed to drought stress during the bolting/bulbing stage. Mapping the sequencing reads to the garlic reference genome and statistical analysis led to the identification of 5215 DEGs, providing insights into the molecular mechanisms involved in drought tolerance.
The primary regulation of gene expression by abiotic stresses occurs at the transcription level [36]. RNA processing follows, wherein pre-mRNAs become mature mRNAs, including the modification of RNA molecules and removal of introns by splicing. Alternative splicing under abiotic stress has been reported in previous RNA-seq studies [37,38], demonstrating that RNA-seq is a powerful method to study gene expression regulation at the transcriptional and post-transcriptional levels. In addition, it should be noted that pervasive expression of the eukaryotic genome produces a great number of non-protein-coding RNAs (ncRNAs) as well as protein-coding mRNAs [39]. ncRNAs occur at the low level in plant cells, but accumulating evidence indicates their involvement in plant growth, reproduction, and abiotic stress responses [40][41][42]. In a previous study, 57,561 protein-coding genes and 20,008 non-coding RNAs were identified in garlic [19]. Although the focus of this study was to identify protein-coding genes that were differentially expressed in response to drought stress, investigating differentially expressed ncRNAs under the same conditions is an interesting topic for future studies. ABA levels are induced by drought stress and play essential roles in drought tolerance [43]. The enhanced accumulation of ABA in response to drought is mainly due to elevated transcription of vital carotenoid and ABA biosynthetic genes. For example, BCH was found to be upregulated by water stress in rice [44]. NCED was strongly induced by drought stress in Arabidopsis [45], cowpea [46], peanut [47], wheat [48], rice [4], and tomato [49]. ABA4 expression was necessary to dehydration-induced ABA biosynthesis [28]. ABA levels in response to drought stress are also determined by ABA catabolism. ABA 8 -hydroxylase (CYP707A3) was found to determine ABA threshold levels during dehydration and after rehydration [50]. In response to drought, BCH, NCED, and ABA4 levels were significantly upregulated in garlic, whereas ABA 8 -hydroxylase level was downregulated. The reverse regulation of ABA biosynthesis and catabolism possibly led to enhanced accumulation of ABA in leaves, revealing that garlic probably utilizes the same mechanism to cope with drought stress. These genes can serve as targets for improving the tolerance of garlic plants to drought stress.
Drought and other abiotic stresses often cause protein denaturation and aggregation, which is toxic to plant cells [8]. Therefore, maintaining the conformation of proteins and preventing the aggregation of misfolded ones is critical for plant survival under drought or other stresses. One strategy that plants utilize is to produce HSPs, which facilitate refolding of denatured proteins [9,10]. Although HSPs are commonly induced by heat stress [51], they have also been reported to be upregulated by other abiotic stresses [52,53]. For example, OsHSP50.2, sHSP17.7, AtHsp90.2, AtHsp90.5, and AtHsp90.7 overexpression was found to enhance the tolerance of rice and Arabidopsis to drought, indicating that HSP production is one of adaptive mechanisms used by plants to tackle drought [54,55]. Likewise, we observed that 124 genes encoding various HSPs were significantly induced by drought stress in garlic (Figure 7 and Supplementary Table S4), indicating that refolding of denatured proteins is perhaps an essential mechanisms induced by garlic to tolerate drought stress.
The ubiquitin/26S proteasome is a major protein degradation system in the nucleus and cytosol of plant cells and is involved in various processes, such as light signaling, plant hormone signaling, and abiotic stress responses [56,57]. In this system, E1, E2, and E3 ubiquitin ligases mediate ubiquitin ligation to target proteins. The genome of Arabidopsis contains two E1, 37 E2, and >1300 E3 genes [58]. A large number of E3 ubiquitin ligases indicates their involvement in diverse processes. Indeed, the expression of 24 E3 ubiquitin ligase genes was induced by drought in garlic leaves (Figure 7 and Supplementary Table S4). These E3 ubiquitin ligases probably have two functions. The first is to remove denatured or dysfunctional proteins generated in response to drought. The α2 subunit of 26S proteasome was found to play a similar role in rice, conferring thermotolerance [59]. The second is to remove suppressors of drought signaling [60], regulating drought tolerance at the post-translational level in garlic.
Drought and salinity are significant abiotic stresses that pose a major threat to agriculture. Common physiological consequences imposed by them include cell dehydration and ROS production [61]. Herein we found that 16.6% (867 out of 5215 genes) of droughtregulated genes were altered by salinity stress in the same manner. In comparison, 5.8% (305 out of 5215 genes) of drought-regulated genes were reversely influenced by salinity. These genes possibly form the molecular basis of dissimilar and similar physiological responses to drought and salinity stresses. Further functional analysis of these genes should provide a novel strategy to produce garlic plants with better tolerance to both drought and salinity or improve drought or salinity tolerance without detrimental effects on other stress defense mechanisms.

Conclusions
Transcriptomic analyses of hardneck garlic at the bolting/bulbing stage suggested that ABA biosynthesis and protein homeostasis possibly play essential roles in drought tolerance, furthering our understanding of the underlying molecular mechanisms.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-439 5/11/2/246/s1, Table S1: Sample sequencing statistics, Table S2: Summary statistics of the RNA-seq results, Table S3: Garlic genes identified in the study, Table S4: The DEGs identified in the garlic leaves after drought treatment, Table S5: GO terms analysis of the DEGs from drought-treated leaves, Table S6: KEGG pathways enriched in drought treatment.