Two black sesame cultivars, Jinhuangma (JHM) and Poyanghei (PYH), widely cultivated and used in China were assessed in this study. They were provided by the Key Laboratory of Crop Physiology, Ecology, and Genetic Breeding, Ministry of Education, College of Agronomy, Jiangxi Agricultural University (Nanchang, Jiangxi province, China). These varieties were selected based on their cultivation history and performance in fields under various environmental conditions. In fact, our group screened hundreds of sesame plant materials via field experiments and selected these two native varieties to clarify the mechanisms underlying black sesame tolerance to drought stress. The cultivar JHM was relatively tolerant to drought, while PYH was sensitive and possessed a high per plant yield in the north of Jiangxi province.

The experimentation was performed in a greenhouse at Jiangxi Agricultural University in 2020. Seeds of the two cultivars were sown float tray until the two true-leave stages. Then, they were transferred into plastic pots (Diameter×Height: 19×27 cm) containing 7.5 kg of soil. To better control water status, the soil water content was measured using a soil moisture probe, 20 cm long (probes inserted vertically into the pots). All the sesame seedlings were watered normally (soil watered daily to 25 ± 5%) until they reached the flowering period (85 days after germination). Thereafter, the drought stress was imposed on other treatments by keeping the soil moisture content at the level of 10 ± 5% for three (T1), five (T2), and seven days (T3), which stand for mild, moderate, and severe drought stress, respectively. Each treatment was composed of twenty individual plants. Middle leaf and root samples from drought stress were sampled at the end of each treatment for physiological analyses. Prior to the induction of the drought stress (the starting day), middle leaves and roots were sampled to constitute control samples (CK) in order to investigate changes in a dynamic manner. Other root samples were prepared for RNA sequence. All the samples were immediately frozen in liquid nitrogen and stored at -80 °C until use. Each sample was analyzed in triplicate.

The chlorophyll content of leaves was determined using the SPAD meter as described previously (Naus et al., 2010; Ling et al., 2011). MDA content was measured via the thiobarbituric acid method (Castrejón and Yatsimirsky, 1997). The activity of SOD and CAT was assayed following the method described by García-Triana (García-Triana et al., 2010) and Zhao and Shi (Zhao and Shi, 2009), respectively. POD activity was measured as described in a previous report (Wang et al., 2014). Free proline content was determined based on the spectrophotometric method described by Vieira (Vieira et al., 2010). Soluble sugar content was determined using the anthranone reagent (Bodelón et al., 2010). Finally, soluble protein content was measured by the coomassie brilliant blue G-250 (Lowry et al., 1951).

Total RNA from root samples was extracted using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer and checked using RNase-free agarose gel electrophoresis. Next, each sample RNA was PCR amplified and sequenced using the Illumina Novaseq6000 Sequencing System by Gene Denovo Biotechnology Co. (Guangzhou, China). To get high-quality clean reads, reads were further filtered by fastp (Chen et al., 2018). The mapped reads of each sample were assembled with StringTie v1.3.1 in a reference-based approach (Pertea et al., 2015; Pertea et al., 2016). For each transcript an FPKM value was calculated to quantify its expression abundance and variations using RSEM software (Dewey and Bo, 2011). Principal component analysis (PCA) was performed in R using the “prcomp” package.

The DESseq2 software (Love et al., 2014) was used to detect DEGs between two different groups with the criteria of false discovery rate (FDR) below 0.05 and absolute fold change ≥ 2. The HISAT2 program (Kim et al., 2015) was used to align the clean reads to the sesame reference genome (“S_indicum_v1.0”, https://www.ncbi.nlm.nih.gov/data-hub/taxonomy/4182/) and to obtain information regarding genomic loci and characteristics unique to the sequenced samples. GO (Genes Ontology, http://geneontology.org/), and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/kaas) enrichment analysis for the DEGs were performed using GO seq and KOBAS (2.0) software, respectively. According to the GO annotation result, the DEGs were mapped to GO terms in the Gene Ontology database, and significant enrichment terms were detected at the threshold P-value <0.05. Similarly, KEGG pathways were assigned to the assembled sequences using the online KEGG Automatic annotation server and enrichment analysis.

We isolated total RNA from each root sample and synthesized first-strand cDNAs following the reported methods by Wei et al. (Wei et al., 2019). Real-time quantitative PCR (RT-qPCR) was carried out in CFX96 (BioRad) with the SYBR Green Perfect mix (TaKaRa, Dalian, China). All samples were analyzed in triplicates. Relative expression levels of each gene were computed using the 2^–ΔΔCT method (Livak and Schmittgen, 2001). The sesame gene β-actin (ncbi_105159390) was used to normalize the genes’ expression levels (Li et al., 2017; Su et al., 2022). The primers were designed with Primer5 and are listed in Table S1 .

Statistical analyses of all traits were conducted using SPSS 17.0 software, and the data are presented as the mean ± SD of three replicates. The standard error is shown as an estimate of variability, and Duncan’s multiple test was used to determine statistical differences at P < 0.05.

To access the drought-responsive mechanisms of JHM and PYH during anthesis, we investigated morphological and physiological changes after three (T1), five (T2), and seven (T3) days of stress induction. Morphological observations showed that the drought stress caused symptoms of yellowing, drooping, and wilting of leaves of plants of both cultivars ( Figure S1 ). However, compared to PYH plants, JHM plants were less affected at T3 ( Figures S1C–F ), indicating they have suffered less damage from drought stress. The yellowing symptom is generally caused by a decrease in chlorophyll content of plants exposed to drought stress (Bhargava and Sawant, 2012; SeyedYahya and Hamideh, 2016; Gurumurthy et al., 2019). Supportively, leaves chlorophyll content analysis revealed a significant decrease in chlorophyll content of the two cultivars ( Figure S2 ). For instance, after five days of stress exposure, the leaf chlorophyll content of JHM and PYH plants exhibited a decrease of 23.75% (from 39.87 to 30.40 SPAD) and 30.78% (from 41.80 to 28.93 SPAD), respectively.

As drought stress promotes the synthesis of oxidants, we analyzed the malondialdehyde (MDA) content and the activity of antioxidative enzymes superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) in leaves and roots of the two cultivars at the different time points ( Figures 1 , S3 ). The MDA content of the leaves and roots of the two cultivars was significantly increased along with the drought stress duration ( Figures 1A , S3A ). Compared to JHM, the MDA content of PYH was significantly higher, indicating that the degree of membrane lipid peroxidation was more severe in PYH plants. The activity of SOD in roots and leaves of the two sesame cultivars significantly increased up to T2 and then decreased, while POD and CAT activities were increased along with the stress duration ( Figures 1B-D , Figures S3 B–D ). It is worth noting that the activity of antioxidant enzymes in JHM under drought stress conditions was significantly higher than in PYH, implying that JHM had a stronger enzymatic defense system than PYH. For instance, the activities of SOD, POD, and CAT in the leaf of JHM reached a maximum value of 166.10, 44.4, and 257.4 u/g FW (fresh weight), respectively, under the drought stress conditions compared to 119.5, 37.7, and 182.9 u/g FW, respectively, in PYH ( Figures 1B–D ).

We further assayed the content of osmolytes, including soluble sugar (SS), soluble protein (SP), free proline (PRO), and glutathione (GSH) in the leaf and root of the two sesame cultivars under the drought stress conditions at the different time points. Except for the roots’ SS content, both the osmolytes showed a significant increase of content in leaves and roots at T1 and T2 and then decreased at T3 ( Figure 2 ), suggesting that prolonged drought stress of more than a week might severely affect sesame productivity. As expected, the increase in osmolytes contents in JHM was more considerable than in PYH ( Figure 2 ).

To get more insights into the drought-responsive mechanisms in JHM and PYH plants, roots samples of CK, T1, T2, and T3 from the two sesame cultivars in three biological replicates, were subjected to RNA sequencing via the Illumina sequencing platform. The summary of the transcriptomics data is presented in Table S2 . The total clean reads generated varied from 38.01 to 55.45 million. The unique mapping reads matching the sesame reference genome ranged from 77.71 to 93.48% and 78.76 to 93.94% for JHM and PYH, respectively. Correlation analysis showed strong positive correlations between samples within the same group, indicating high reproducibility between the biological replicates ( Figures S4A ). We performed the principal component analysis (PCA) to differentiate between the groups. The results showed that the transcriptome of JHM and PYH roots was very different and changed according to the drought stress severity ( Figure 3 ). There were 22674 genes (including 1142 novel genes) in all samples ( Table S3 ). The total sequenced genes from all samples accounted for 94.45% of the sesame reference genome ( Table S3 ). Reads alignment analysis showed that over 75% of the genes are located in exonic regions ( Figure S5 ).

To uncover drought-induced changes in transcriptional levels in the two cultivars during anthesis, we carried out differentially expressed genes (DEGs) analysis along with the drought duration. In total, we identified 24,037 DEGs, including 13,951 and 10,086 up- and down-regulated genes in JHM, respectively. Meanwhile, 23,604 DEGs, including 12,586 and 11,018 up- and down-regulated genes, respectively, were identified in PYH. In both cultivars, the number of up-regulated DEGs showed similar patterns of increasing and then decreasing from T2 ( Figure 4A ). In contrast, the number of down-regulated DEGs increased along with the drought stress duration. We searched for genes that were significantly affected at all time points in JHM and PYH. We detected 3,881 and 3,409 overlapped DEGs at the different time points in JHM and PYH, respectively ( Figures 4B, C ).

We performed GO and KEGG analyses to unveil molecular mechanisms involving the DEGs. The most GO terms that involve DEGs at early (T1) and moderate (T2) drought stress during anthesis in JHM included RNA modification, organic substance metabolic process, and ribonucleoprotein complex biogenesis ( Figure S6A ). At T3, the main GO terms that involve DEGs in JHM were hormone-mediated signaling pathway, regulation of meristem development, response to hormones, response to stimulus, signal transduction, and developmental growth. While in PYH, the most identified GO terms were stomatal movement, response to water deprivation, proteolysis, RNA modification, anion transport, and acetyl-CoA metabolic process at early and moderate drought stages, and response to water deprivation, ethylene metabolic process, anion transport, and steroid biosynthesis process at T3 ( Figure S6C ). Noteworthily, GO term related to the stomatal movement was identified only in PYH at the early drought stage and included four DEGs, ncbi_105159372 (PLDDELTA), ncbi_105162005 (AHK5), ncbi_105162425(CAS), and ncbi_105174436 (HSC-2).

The KEGG analysis revealed that several pathways related to plants’ drought stress tolerance mechanisms were significantly induced in JHM compared to PYH ( Figures S6B–D ). For instance, at T3, the DEGs in JHM were mainly assigned to photosynthesis - antenna proteins, biosynthesis of amino acids, fatty acid metabolism, peroxisome, lysine degradation, ascorbate and aldarate metabolism, pyruvate metabolism, plant hormone signal transduction, tryptophan metabolism, biosynthesis of secondary metabolites, and glutathione metabolism ( Figure S6B ). In contrast, in PYH, the most significant pathways at the same time were plant hormone signal transduction, MAPK signaling pathway, carbon metabolism, biosynthesis of amino acids, and arachidonic and alpha-linolenic acids metabolism ( Figure S6D ).

A total of 5,453 DEGs, including 505 common (up- or down-regulated), were identified between JHM and PYH ( Figures 4A , 3A ). GO analysis revealed that at T2 and T3, most DEGs between JHM and PYH plants were related to vacuolar transport, chemical homeostasis, response to abiotic stimulus, response to oxidative stress, and secondary metabolites biosynthesis processes ( Figure 3D ). Meanwhile, the most induced pathways between JHM and PYH plants were photosynthesis, carotenoid biosynthesis, and biosynthesis of secondary metabolites ( Figure 3E ). Flavonoid biosynthesis and starch and sucrose metabolism were specifically significantly induced at T3. Interestingly, the DEGs enriched in photosynthesis were 2.0- to 7.9-fold more highly induced in JHM than in PYH.

In order to identify potential candidate genes for targeted improvement of drought stress tolerance in sesame, we constructed Venn diagrams among the up- and down-regulated DEGs between JHM and PYH. There were 54 and 42 DEGs significantly induced in JHM and PYH, respectively, at the four-time points ( Figures 3B, C ). It is worth noting that the expression of 11 up- and 11 down-regulated genes between JHM and PYH were highly affected (|FPKM| >10 during at least one-time point) along with the drought treatments. Thus, we selected these genes as potential candidate genes for future studies aiming to enhance drought stress tolerance in sesame ( Table S4 ).

Glutathione and ethylene play critical roles in plants’ tolerance to abiotic stresses. Glutathione reductase (GR) catalyzes to maintain cellular levels of reduced glutathione, which is essential for reactive oxygen species control (Couto et al., 2016). We then examined the expression of genes involved in the glutathione and ethylene biosynthesis pathway ( Figure 5A ), including At3g24170 (GR, glutathione reductase, ncbi_105158649), S-adenosylmethionine-dependent methyltransferase (LAMT, ncbi_105155206), 1-aminocyclopropane-1-carboxylate synthase (ACS1, ncbi_105164055), and 1-aminocyclopropane-1-carboxylate oxidase 1 (ACO1, ncbi_105161839). As shown in Figure 5B , the sesame GR was up-regulated along with the drought stress duration in both JHM and PYH plants. However, it was slightly more induced in JHM plants than in PYH plants. The expression of LAMT was also induced by the drought stress in the two cultivars. ACS1 was mainly induced at T3, while ACO1 at T1 ( Figure 5B ). At T3, the expression levels of the four genes were 1.0- to 8.1-fold highly induced in JHM plants than in PYH plants.

We screened the expressed genes and found that the most expressed transcription factor families at all time points of the drought stress in JHM plants were ERF (20%), followed by NAC (15%), MYB (10%), GeBP (10%), C₃H (10%), and bHLH (10%) ( Figure 6A ). Meanwhile, NAC (20%), C₃H (15%), MYB (10%), bZIP (10%), and bHLH (10%) were the most expressed TFs in PYH plants ( Figure 6B ). We then filtered out the top expressed TFs in the two cultivars under the drought condition for future studies ( Figure 6C ). In JHM, the top expressed TFs included ERF091 (ncbi_105175644), LAF1 (ncbi_105160487), JUB1 (ncbi_105175173), NAC100 (ncbi_105157091), and MYB4 (ncbi_105176094). These genes are known to be essential for plant survival from abiotic stresses. For example, JUB1 is involved in various metabolic processes, such as trehalose, proline, hyperosmotic, and flavonoid biosynthetic processes (Alshareef et al., 2019; Chen et al., 2021).

To validate the RNA-seq data, we selected seven genes for RT-qPCR analysis. As shown in Figure 7 , the expression patterns of the selected genes via the RT-qPCR and RNA-seq were consistent (R² = 0.85), confirming the reliability of our results.