RNA-Seq Analysis Demonstrates Different Strategies Employed by Tiger Nuts (Cyperus esculentus L.) in Response to Drought Stress

Drought stress, an important abiotic stress, has affected global agricultural production by limiting the yield and the quality of crops. Tiger nuts (Cyperus esculentus L.) are C4 crops in the Cyperaceae family, which have high-quality wholesome ingredients. However, data on mechanisms underlying the response of tiger nuts to drought stress are few. Here, the variety of Jisha 1 and 15% polyethylene glycol (PEG; a drought stress simulator) were used to study the mechanisms of stress response in tiger nuts. Our evaluation of the changes in physiological indicators such as electrolyte leakage (El), malondialdehyde (MDA), hydrogen peroxide (H2O2), superoxide anion (O2−) and activities of reactive oxygen species (ROS) showed that 12 h was the most suitable time point to harvest and analyze the response to drought stress. Thereafter, we performed transcriptome (RNA-Seq) analysis in the control (CK) and stress treatment groups and showed that there was a total of 1760 differentially expressed genes (DEGs). Gene Ontology (GO) analysis showed that the DEGs were enriched in abscisic acid (ABA) terms, and pathways such as starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis (ko00940) and plant hormone signal transduction (ko04075) were significantly enriched in the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. In addition, quantitative real-time PCR (qRT-PCR) analysis of the DEGs demonstrated an upregulation of ABA and lignin content, as well as enzyme activities in enriched pathways, which validated the RNA-Seq data. These results revealed the pathways and mechanisms adopted by the tiger nuts in response to drought stress.


Introduction
Tiger nut (Cyperus esculentus) is a C4 plant of the Cyperaceae family, which has a high yield and contains underground tubers with high nutrient value [1]. The tiger nut crop is widespread in tropical, temperate, and cooler zones [2] and originated from the Mediterranean area. It was domesticated and utilized as an important food source in Ancient Egypt [3,4]. The underground tubers of the tiger nuts do not only have a huge amount of nutrients, including starch, oil, sugars, protein, dietary fibers, vitamins C and E and minerals [4,5], but bioactive substances, such as phytosterols, alkaloids, saponins, tannins, flavonoids and terpenoids, as well [6,7]. These ingredients have beneficial effects in patients with diabetes, cardiovascular disease and obesity [7,8]. Since the tiger nut is

Plant Materials and Treatments
The "Jisha 1" plant material was provided by the Economic Botany Institute of Jilin Academy of Agricultural Sciences. Plump seeds were sterilized with 10% NaClO for 5 min and rinsed with distilled water five times. Then, the seeds were placed in a 30 • C incubator without light for five days for cultivation. The consistent growth sprouts were selected as plant materials for further experiments. A 15% PEG concentration (the solvent was distilled water) was used to simulate drought stress (D), while distilled water treatment was set as control treatment (CK) [29]. The samples were collected after 0 h, 3 h, 6 h, 12 h, 24 h and 48 h of treatment, respectively. Ten sprouts in the same treatment were deemed as an experimental unit. The experiments were replicated five times.

Analysis of Physiological Indicators
Samples of 0.5 g of fresh roots were used to test the trends of the physiological indicators, which were associated with abiotic stress in plants. Electrolyte leakage (El) of samples was tested using the methods of Sutinen [30], while the malondialdehyde (MDA) of each sample was measured using the ELISA Kit (M0106, Michy, Suzhou, China) following the manufacturer's instructions; Hydrogen peroxide (H 2 O 2 ) concentration was tested using ELISA Kit (M0107, Michy, Suzhou, China) following the manufacturer's instructions [31]; and superoxide anion (O 2 − ) was measured using ELISA Kit (M0114, Michy, Suzhou, China) following the manufacturer's instructions. All of these physiological indicators had five replicates.
In addition, 0.5 g fresh roots samples were used to test the trends of the activities in the ROS system. The activities of POD, SOD, CAT and APX were extracted by ELISA Kits. The activities of POD were determined using ELISA Kit (M0105, Michy, Suzhou, China) following the manufacturer's instructions; the activities of SOD were tested using ELISA Kit (M0101, Michy, Suzhou, China) following the manufacturer's instructions; the activities of CAT were measured using ELISA Kit (M0103, Michy, Suzhou, China) following the manufacturer's instructions; and the activities of APX were determined using ELISA Kit (M0403, Michy, Suzhou, China) following the manufacturer's instructions. The activities of these enzymes were analyzed in five replicates while the data of the activities were read by a microplate reader (SpectraMax ® 190, Molecular Devices, Los Angeles, CA, USA).

RNA-Seq Analysis
The tender roots at the bud stage under two treatments for 12 h were used as samples for RNA-Seq analysis, and each treatment had three biological replicates. These samples were rapidly frozen in liquid nitrogen. And then RNA was extracted by RNApure Plant Kit (CW0559, CWBIO, Beijing, China) following the manufacturer's instructions. The RNA was measured in 1% agarose gel electrophoresis and NanoDrop instrumentation (OneC, Thermo, Waltham, MA, USA) using the methods of Zhang [20], which determined whether the quality of the RNA was qualified. The qualified RNA of these samples underwent RNA-Seq by Bio-maker (Beijing, China), while those with no reference genome were analyzed by Biocloud (https://international.biocloud.net/, accessed on 1 January 2022). UniGene was analyzed by Trinity software [32], while the differentially expressed genes (DEGs) were assessed by DESeq2 software [33]. In addition, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used for gene and pathway enrichment analyses [34].

QRT-PCR Analysis
The extracted RNA in RNA-Seq was also used to synthesize single-strand cDNA using Evo M-MLV RT Premix for qPCR (AG11706, Accurate Biology, Hunan, China). The qRT-PCR was performed in a Light Cycler 480II system (Roche, Roche Diagnostics, Basel, Switzerland) using Hieff UNICON ® Universal Blue qPCR SYBR Green Master Mix (11184E, Yeasen, Shanghai, China). The reversed transcription products were diluted 10-fold and then used for qRT-PCR in a 20 µL reaction volume. While the qRT-PCR conditions were used according to the manufacturer's instructions (Table S1). Primers for DEGs were designed by the NCBI database, and CeUCE2 was used as reference for the qRT-PCR [34] analysis (Table S2). The relative expression levels of the selected DEGs were calculated with biological and technological replicates, using the 2 −∆Ct method [17].

Validation of the Physiological Indicators in Enriched Pathways
The ABA content was determined by HPLC-MS/MS (AB SCIEX, shimadzuLc-20AD and AB5500, Framingham, MA, USA) at the Quality Inspection Center (Dalian, China) [35]. The physiological indicators (such as starch, sucrose and lignin) [36] and the enzyme activities [37] in the enriched pathways were compared between the control and stress treatment groups and were determined by ELISA Kit (Michy, Suzhou, China) with five replicates and a microplate reader. Data were analyzed using SPSS19.0 (Armonk, NY, USA), and a p < 0.05 was used to demonstrate statistical significance.

The Physiological Trends under Stress
The roots of tiger nuts at the bud stage under control (CK) and drought stress were analyzed for the 0 h, 3 h, 6 h, 12 h, 24 h and 48 h time points. The membrane lipid peroxidation indicators such as El, MDA, H 2 O 2 and O 2 − were tested for stress trends. The results showed that the four indicators had increased significantly under drought stress compared with CK between 3 h and 12 h, but growth slowed after 12 h ( Figure 1A-D). These findings illustrated that 12 h was the optimum time point to harvest and analyze the mechanisms of stress response by the plants.

The Physiological Trends under Stress
The roots of tiger nuts at the bud stage under control (CK) and drought stress were analyzed for the 0 h, 3 h, 6 h, 12 h, 24 h and 48 h time points. The membrane lipid peroxidation indicators such as El, MDA, H2O2 and O2 − were tested for stress trends. The results showed that the four indicators had increased significantly under drought stress compared with CK between 3 h and 12 h, but growth slowed after 12 h ( Figure 1A-D). These findings illustrated that 12 h was the optimum time point to harvest and analyze the mechanisms of stress response by the plants. In addition, we determined the four enzymatic activities in the ROS (including SOD, POD, CAT and APX) at the defined sampling time points. The assessment of the activities of SOD, POD, CAT and APX showed that drought stress significantly increased the enzymes' activities compared with CK, and the increment slowed after 12 h (Figure 2A-D). These findings were in sync with the indicators in membrane lipid peroxidation, thus demonstrating that 12 h is an optimum sampling time point for studying the mechanisms of response to drought stress. In addition, we determined the four enzymatic activities in the ROS (including SOD, POD, CAT and APX) at the defined sampling time points. The assessment of the activities of SOD, POD, CAT and APX showed that drought stress significantly increased the enzymes' activities compared with CK, and the increment slowed after 12 h (Figure 2A-D). These findings were in sync with the indicators in membrane lipid peroxidation, thus demonstrating that 12 h is an optimum sampling time point for studying the mechanisms of response to drought stress.

RNA-Seq Analysis
There were 5.74 Gb of RNA-Seq data in each sample, while Q20 was more than 97.91% and Q30 was more than 94.04% (Table S3). The quality of the RNA-Seq data was confirmed and then used for further analysis. The data were uploaded in National Centre for Biotechnology Information (NCBI) database, with the accession number PRJNA821655

RNA-Seq Analysis
There were 5.74 Gb of RNA-Seq data in each sample, while Q20 was more than 97.91% and Q30 was more than 94.04% (Table S3). The quality of the RNA-Seq data was confirmed and then used for further analysis. The data were uploaded in National Centre for Biotechnology Information (NCBI) database, with the accession number PRJNA821655 (Table S4).
The UniGene dataset was used as a reference for further analyses while the DEGs were assessed by DESeq2 ( Table 1). Out of the total DEGs, 854 DEGs were upregulated, while 906 were downregulated (Figure 3).

Enrichment Analysis
These DEGs underwent GO and KEGG enrichment analysis. In the GO enrichment analysis, the DEGs were enriched in three classes ( Figure 4A). The q-value of 55 GO terms was less than 0.05, which were the significantly enriched terms, while four terms were associated with their response to ABA. These included the cellular response to abscisic acid stimulus (GO:0071215), the response to abscisic acid (GO:0009737), the cellular response to hormone stimulus (GO:0032870) and the abscisic-acid-activated signaling pathway (GO:0009738). The GO analysis demonstrated that ABA participates in responding to stress. On the other hand, the KEGG enrichment analysis showed that there was significant enrichment of DEGs in four KEGG pathways (q-value < 0.05), which included ribosome (ko03010), starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis

Enrichment Analysis
These DEGs underwent GO and KEGG enrichment analysis. In the GO enrichment analysis, the DEGs were enriched in three classes ( Figure 4A). The q-value of 55 GO terms was less than 0.05, which were the significantly enriched terms, while four terms were associated with their response to ABA. These included the cellular response to abscisic acid stimulus (GO:0071215), the response to abscisic acid (GO:0009737), the cellular response to hormone stimulus (GO:0032870) and the abscisic-acid-activated signaling pathway (GO:0009738). The GO analysis demonstrated that ABA participates in responding to stress. On the other hand, the KEGG enrichment analysis showed that there was significant enrichment of DEGs in four KEGG pathways (q-value < 0.05), which included ribosome (ko03010), starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis (ko00940) as well as plant hormone signal transduction (ko04075) ( Figure 4B). These pathways might be mediating the responses to drought stress.

QRT-PCR Analysis
Twelve DEGs were selected for qRT-PCR analysis. The data showed that three parts of four DEGs were involved in three enriched pathways, including ABA signal transduction, starch and sucrose metabolism and phenylpropanoid biosynthesis. The qRT-PCR re-

QRT-PCR Analysis
Twelve DEGs were selected for qRT-PCR analysis. The data showed that three parts of four DEGs were involved in three enriched pathways, including ABA signal transduction, starch and sucrose metabolism and phenylpropanoid biosynthesis. The qRT-PCR results showed differences in the expression compared with the RNA-Seq findings, while the trends in CK and drought stress were similar ( Figure 5).

Enriched Pathways Analysis
The plant hormone signal transduction pathway (ko04075) was enriched in both the GO and KEGG analyses. Similarly, the ABA channel was also enriched ( Figure 6A), together with the DEGs as shown in Figure 6B. In addition, the data showed that the ABA content had increased significantly (p < 0.05) under drought stress ( Figure 6C), which showed that the ABA-specific pathway was enriched in response to stress.

Enriched Pathways Analysis
The plant hormone signal transduction pathway (ko04075) was enriched in both the GO and KEGG analyses. Similarly, the ABA channel was also enriched ( Figure 6A), together with the DEGs as shown in Figure 6B. In addition, the data showed that the ABA content had increased significantly (p < 0.05) under drought stress ( Figure 6C), which showed that the ABA-specific pathway was enriched in response to stress. The plant hormone signal transduction pathway (ko04075) was enriched in both t GO and KEGG analyses. Similarly, the ABA channel was also enriched ( Figure 6A), t gether with the DEGs as shown in Figure 6B. In addition, the data showed that the AB content had increased significantly (p < 0.05) under drought stress ( Figure 6C), whi showed that the ABA-specific pathway was enriched in response to stress.  KEGG analysis showed that the starch and sucrose metabolism pathway (ko00500) was also enriched ( Figure 7A). The expression of DEGs in the ko00500 pathway was significantly changed ( Figure 7B). In addition, there was an increase in the sucrose content while the starch content was suppressed with the increasing activities significantly (p < 0.05) of alpha-amylase and beta-amylase ( Figure 7C-F). Therefore, these results revealed that the starch and sucrose metabolism pathway was involved in responses to drought stress. was also enriched ( Figure 7A). The expression of DEGs in the ko00500 pathway was significantly changed ( Figure 7B). In addition, there was an increase in the sucrose content while the starch content was suppressed with the increasing activities significantly (p < 0.05) of alpha-amylase and beta-amylase ( Figure 7C-F). Therefore, these results revealed that the starch and sucrose metabolism pathway was involved in responses to drought stress.  In addition, there was enrichment in phenylpropanoid biosynthesis pathway (ko00940) (Figure 8A), and the expression of DEGs in the ko00500 pathway had significant changes ( Figure 8B). Similarly, lignin content had increased significantly (p < 0.05) and there were significant shifts (p < 0.05) in the activities of cinnamyl-alcohol dehydrogenase (CAD) (Figure 8C,D). These results revealed that phenylpropanoid biosynthesis was one of the pathways that mediates responses to stress. represent upregulated and downregulated DEGs contain; (B) The expression profile of DEGs in modules enriched in pathways, the color from blue to red represent the expression from low to high; (C) The sucrose content analysis under CK and drought treatments; (D) The starch content analysis under CK and drought treatments; (E) The activities of alpha-amylase analysis under CK and drought treatments; (F) The activities of beta-amylase analysis under CK and drought treatments.
In addition, there was enrichment in phenylpropanoid biosynthesis pathway (ko00940) (Figure 8A), and the expression of DEGs in the ko00500 pathway had significant changes ( Figure 8B). Similarly, lignin content had increased significantly (p < 0.05) and there were significant shifts (p < 0.05) in the activities of cinnamyl-alcohol dehydrogenase (CAD) (Figure 8C,D). These results revealed that phenylpropanoid biosynthesis was one of the pathways that mediates responses to stress.

Discussion
The bud stage is the most sensitive stage in plants under abiotic stresses, which could affect growth in later stages [38]. Drought stress limits the growth, development and survival rate of plants, leading to tremendous losses in yield, thus threatening global agricultural production [39]. In our study, 15% PEG was used to simulate drought stress [29,[40][41][42]. Therefore, the analysis of the response of sprouts to drought stress at the bud stage is important for actual agricultural production [43]. Previous studies have shown that roots can be used to study responses to drought stress in common bean (Phaseolus vulgaris) [44], rice (Oryza sativa) [45], soybean (Glycine max) [46] and peanut (Arachis hypogaea) [47]. In this study, 15% PEG was used to simulate the drought stress and roots of the tiger nuts were used as sample tissues, as previously described [48,49]. Drought-exposed plants exert different mechanisms to regulate the responses to stress [39]. A 12 h duration under drought stress has been shown to be an optimal sampling time to study stress responses [50]. Exposing plants to stress for a short time does not cause changes in the phenotype, but physiological indicators might change significantly [51]. El is an indicator which reflects the state of a plant's membrane system [52]. El increases when plants are subjected to stress or other damaging stimuli [53]. The El was shown to increase under drought stress, but the rate of increase after 12 h was not significant. MDA content is a parameter which reflects antioxidant potential, and indirectly reflects the degree of tissue peroxidative damage [54]. There was significant increase in the MDA content under stress before 12 h, which was not observed after the 12 h time point. In addition, H 2 O 2 and O 2 − reflect the plant's resilience under stress [55]. Our study showed that 12 h of exposure to stress increased the levels of H 2 O 2 and O 2 − . On the other hand, the ROS system plays a dual role in the growth and development of plants [56], and enzyme activities directly reflect ROS in plants. In this study, the POD, SOD, CAT and APX activities increased more within 12 h of exposure. In line with previous data, these data revealed that 12 h was a suitable sampling time to study the mechanisms involved in responding to drought stress in tiger nuts [50,57,58].
There are many mechanisms involved in responding to abiotic stress, especially drought stress, in plants [59]. Drought stress induces ABA accumulation and triggers rapid biochemical and physiological responses [60]. The ABA pathway is a well-established pathway which has been shown to respond to drought stress [61]. In this study, DEGs were enriched in ABA-associated pathways in the GO and KEGG analyses, while the ABA content changed under drought stress, which demonstrated that ABA-associated pathways participate in responding to drought stress in tiger nuts. The starch and sucrose metabolism pathway is closely associated with ABA accumulation [62][63][64], which was also associated with responses to drought stress [36]. Starch is made in chloroplasts in leaves which the content decreased under drought stress [65], while sucrose is produced in leaves following photosynthesis, which the content increased to balance osmotic pressure within plant cells under drought stress [64]. Also in this study, the starch and sucrose metabolism pathway was an enriched pathway in the KEGG analysis, which showed an increase in the sucrose content with changes in amylase activities (α and β amylase), while the content of starch was reduced. This could be because the starch was converted to sugar, which increased osmotic pressure in responding to stress. In addition, lignin was recognized as having an important role in drought resistance [66] while the CAD enzyme was used to characterize drought tolerance [67]. Here, the phenylpropanoid biosynthesis (ko00940), as a stress-resistant pathway [68], was enriched in the KEGG analysis, while the contents of lignin and CAD activities were significantly changed (Figure 9). In addition, the three pathways have been reported to participate in responding to abiotic stress in Hibiscus cannabinus, Miscanthus and tomato (Solanum lycopersicum) [69][70][71]. Our study demonstrated the mechanisms underlying the drought stress responses by tiger nuts, which enhances drought tolerance. demonstrated the mechanisms underlying the drought stress responses by tiger nuts, which enhances drought tolerance.

Conclusions
Our analyses showed that 12 h was a suitable sampling time for tiger nuts when studying stress responses. A comparative analysis of the RNA-Seq data under drought stress showed a total of 1760 DEGs, where 854 genes were upregulated while 906 genes were suppressed. These DEGs were enriched in ABA-associated terms in the GO analysis, while the starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis (ko00940) and plant hormone signal transduction (ko04075) pathways were significantly enriched in the KEGG analysis. Taken together, our data revealed the mechanisms that underly the responses exerted by tiger nuts to drought stress.
Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Table S1: The qRT-PCR program used in this study; Table S2: QRT-PCR primers for differentially expressed genes (DEGs); Table S3: The quality analysis of the RNA-Seq data; Table S4: The RNA-Seq data in National Centre for Biotechnology Information (NCBI) database.
Author Contributions: Z.M. and Z.W.: methodology, data curation and writing-original draft; J.L. and Y.C.: data curation; H.Y., Y.G., X.Y. and J.Z.: conceptualization and methodology; Y.S. and C.L.: software; S.W. and K.L.: formal analysis and preparation of materials; J.D. and Q.Z.: conceptualization, data curation and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Figure 9. Schematic representation of the mechanisms employed by tiger nuts in response to drought stress at bud stage, where red up-arrows indicate upregulation while green arrows show downregulation.

Conclusions
Our analyses showed that 12 h was a suitable sampling time for tiger nuts when studying stress responses. A comparative analysis of the RNA-Seq data under drought stress showed a total of 1760 DEGs, where 854 genes were upregulated while 906 genes were suppressed. These DEGs were enriched in ABA-associated terms in the GO analysis, while the starch and sucrose metabolism (ko00500), phenylpropanoid biosynthesis (ko00940) and plant hormone signal transduction (ko04075) pathways were significantly enriched in the KEGG analysis. Taken together, our data revealed the mechanisms that underly the responses exerted by tiger nuts to drought stress.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/life12071051/s1, Table S1: The qRT-PCR program used in this study; Table S2: QRT-PCR primers for differentially expressed genes (DEGs); Table S3: The quality analysis of the RNA-Seq data; Table S4: The RNA-Seq data in National Centre for Biotechnology Information (NCBI) database.
Author Contributions: Z.M. and Z.W.: methodology, data curation and writing-original draft; J.L. and Y.C.: data curation; H.Y., Y.G., X.Y. and J.Z.: conceptualization and methodology; Y.S. and C.L.: software; S.W. and K.L.: formal analysis and preparation of materials; J.D. and Q.Z.: conceptualization, data curation and funding acquisition. All authors have read and agreed to the published version of the manuscript.