Insight into gene regulatory networks involved in sesame (Sesamum indicum L.) drought response

Drought is one of the most common environmental stresses affecting crops yield and quality. Plants’ responses to drought are controlled by regulatory mechanisms. Sesame is an important oilseed crop that most likely faces drought during its growth due to growing in semi-arid and arid areas. Despite this importance, there is little information about sesame regulatory mechanisms against drought stress. In this study, highly differentially expressed genes were identified using in-silico transcriptome analysis of two sesame genotypes (one sensitive and one tolerant) under drought stress. Interactions between identified genes and regulators including Transcription Factors (TF) and microRNAs (miRNA) were predicted using bioinformatics tools and then related regulatory networks were constructed. A total of 117 TFs and 133 miRNAs that might be involved in drought stress were identified with this approach. Key regulators of sesame under drought stress were detected by network analysis. Regulatory modules involved in drought response were extracted from the integrated networks using the MCODE algorithm, to explore important relationships. Finally, some of the identified drought-tolerance-related genes were examined using qRT-PCR analysis in two contrasting sesame genotypes under drought conditions. Interestingly, the studied genes showed a different expression pattern in the tolerant genotype compared to the sensitive. Taken together, our results suggest that miR530/TGA1/CER1, miR319/NAC29L, and miR171/SUC2L modules might play key roles in sesame drought tolerance by regulating wax biosynthesis, leaf senescence, and sugar transport during drought, respectively. Overall, the identified drought-related genes including TFs and miRNAs along with their relationships could be valuable candidates for future studies and breeding programs on enhancing sesame tolerance under drought stress.


Introduction
Sesame (Sesamum indicum L.) is an ancient oil crop that was cultivated around the world for many years. It is now identified as an important crop due to its numerous applications, including foods, pharmaceutical, and industrial (Morris 2002). Sesame seeds are an excellent source of protein, oil, and antioxidants. The protein and oil content of the seed was reported to be about 25% and 50% respectively (Bedigian 2010). Some reports suggested possible protective effects of sesame oil consumption against heart disease, diabetes, arthritis, and cancers (Cooney et al. 2001;Yokota et al. 2007;Lin et al. 2014). These nutritional, medicinal, and industrial benefits have increased the importance of sesame and led to an ever-growing demand for it.
Sesame traditionally is cultivated in arid and semi-arid regions of the world and drought is one of the most common abiotic stresses in these areas. So, it is very likely for sesame to be faced with moderate to severe drought stress during at least one or two stages of its growth. Although, sesame significantly has been adapted to dry environments and could regain moisture from the lower layers of the soil by developing its tap root system under drought conditions . However, drought stress mostly limits the potential of sesame production and adversely affects its growth, development, metabolism, and yield (Islam et al. 2016). Therefore, improving drought tolerance could be one of the most important sesame breeding goals.
To improve plants against environmental stresses including drought, their response patterns to desired stress should be identified and studied. The plant's response to drought is controlled by several regulatory mechanisms. Transcriptional regulation by transcription factors is one of these mechanisms that plant used to deal with stresses. Transcription factors (TFs) are proteins that bind to cis-regulatory elements located in the promoter of their target gene and regulate its expression. In recent years, several families of plant TFs including MYC, MYB, bZIP, NAC, WRKY, and AP2/ERF have been identified and proved their key roles in drought tolerance (Joshi et al. 2016). In Arabidopsis, it has been shown that the expression of maize ZmMYB3R improves drought and salinity tolerance . The expression of soybean GmNAC085 transcription factor in Arabidopsis significantly led to improved drought tolerance of transgenic plants by positively regulating growth rate and reducing transpiration and cell membrane damage (Nguyen et al. 2018). It was reported in rice that overexpression of OsNAC5, 6, 7, and 9 genes result in enlarged roots and enhanced drought tolerance (Chung et al. 2018). The SlWRKY58 and SlWRKY72 genes were associated with improved drought tolerance in tomatoes (Karkute et al. 2018;Dossa et al. 2016) introduced AP2si16 as a potential candidate gene in breeding programs for drought tolerance, in a study on the AP2/ERF genes of sesame.
The plant regulatory mechanisms are not limited to the transcriptional level. Post-transcriptional regulation is another prominent plant control process in the face of drought. At this level, we can refer to microRNAs (miR-NAs), non-coding single-stranded RNAs that regulate gene expression using sequence-specific cleavage or translation inhibition of the target transcripts (Li et al. 2013b;Park and Shin 2014). Plant miRNAs derive from their precursors that are transcripts of exon or intron of coding regions or intergenic regions (Voinnet 2009;Nozawa et al. 2012;Budak and Akpinar 2015). After multiple processing steps, the mature miRNA forms from the miRNA gene. Mature miRNA is then loaded onto the argonaute protein and along with some other proteins forms the RISC complex (Voinnet 2009). Eventually, miRNA by base pairing with mRNA would guide RISC to cleave the target or repress its translation (Li et al. 2013b;Park and Shin 2014). In recent years, studies on miRNA-mediated gene regulation of plants under abiotic stresses, including drought have been growing, and their results revealed the crucial role of miRNAs in the plant response to stress (Ferdous et al. 2015). A study on maize showed that miRNAs such as miR474, miR168, miR528, and miR167 might be involved in regulating tolerance mechanisms against drought stress (Wei et al. 2009). In barley, miR169b and miR1432 have been reported to control the expression of drought-responsive genes including NFY-A and calmodulin-related proteins, and thus regulate stress response (Ferdous et al. 2017).
The development of high-throughput technologies in recent years and their combination with bioinformatics tools enabled the researchers to have a clearer vision of drought response mechanisms. However, the large number of targets proposed by these methods is a major challenge for researchers. Network analysis and module detection methods have been proved as efficient approaches to extract key genes from large omics datasets in biological researches, especially gene regulatory studies (Saelens et al. 2018;Sanchez and Mackenzie 2020). Recently, studies on the role of transcription factors in sesame responses to drought stress were increased and have been growing (Dossa et al. 2016;Wang et al. 2018b;Zhang et al. 2018). However, our knowledge about many sesame TFs and their roles under drought stress is limited and insufficient. Meanwhile, most of these researches considered TFs response separately, and collective and network-based studies that focused on the interactions between TFs and stress-related genes have been less seen. On the other hand, studies about sesame miRNAs are rare and their regulatory roles in response to drought stress are unclear. The present study aimed to identify important regulators potentially involved in sesame drought response by integrating in-silico data using network-based approaches. Finally, some of the key relationships detected by module detection algorithms were examined in two contrasting sesame genotypes under drought stress using qRT-PCR analysis.
The quality of the raw data was checked using FASTQC (Andrews 2010). Then, excluding low-quality nucleotides and reads and also adaptor removing were done using Trimmomatic (Bolger et al. 2014). After filtering data, clean reads were mapped to the reference genome sequences of sesame available at NCBI Genome database (https:// www. ncbi. nlm. nih. gov/ genome) under genome ID 11,560 using HISAT2 (Kim et al. 2019). Count matrix was calculated based on matched reads using HTSeq (Anders et al. 2015). Count data were normalized using the TMM method included in the edgeR package (Robinson et al. 2010) under the R platform (R Development Core Team 2020). Differentially expressed genes (DEGs) analysis was performed on normalized data using the edgeR package (Robinson et al. 2010).
Highly DEGs with FDR<0.01 and |logFC|>3 were retained for each sample. Genes that were highly differentially expressed in all samples were selected as droughtresponsive core (DRC) gene set. On the other hand, genes that only highly differentially expressed in tolerant samples under severe drought (which include the highly DEGs that were present only in T4 or those were common between T3 and T4 and were not present in T1 and T2) but not in the sensitive samples at the same condition (S4) were selected as drought tolerance-related (DTR) gene set. The analysis was then performed separately on these two sets of genes.

Gene regulatory networks construction
Sesame TFs that targeted gene sets were identified using the TF Enrichment tool of the PlantRegMap database (Jin et al. 2017). TFs that possess a significantly over-presented number of targets in gene sets based on Fisher's exact tests were selected (adjusted p-value cut-off= 0.05). TF-drought related gene (DRG) networks were constructed using regulatory relations between selected TFs and genes of each set.
A list of mature sesame miRNA sequences was obtained from the Sinbase database (Wang et al. 2015a and previous studies (Joshi and Mandavia 2018;Marakli 2018). miRNA that targeted gene sets were predicted by psRNA-Target web-based tool (expectation-value cut-off= 5) (Dai and Zhao 2011;Dai et al. 2018). Using predicted miRNAtarget gene interactions miR-DRG regulatory network was constructed for each set.
The bidirectional interactions of identified miRNAs and TFs with each other were predicted using the Regulation Prediction tool of the PlantRegMap database (TFs target miRNAs) and psRNATarget website (miRNAs target TFs). The integrated miR-TF-DRG gene regulatory networks (GRNs) were constructed by merging TF-DRG, miR-DRG, and miR-TF networks of each set.

Network analysis and modules detection
Networks were visualized by Cytoscape software (Smoot et al. 2011). Centrality measures including degree and betweenness centrality were calculated for GRNs using CytoNCA, a Cytoscape plug-in (Tang et al. 2015). The Molecular Complex Detection (MCODE) algorithm (degree cutoff=2, node score cutoff=0.2, k-Core=2) (Bader and Hogue 2003) was applied to detect network modules using clusterMaker, a Cytoscape clustering plug-in (Morris et al. 2011).

Plant growth and drought treatments
Seeds of two sesame genotypes with SaSiG004 (droughtsensitive) and SaSiG006 (drought-tolerance) accession IDs were selected from a set of sesame genotypes classified based on drought tolerance in the plant breeding department of Sari Agricultural Sciences and Natural Resources University. The seeds were sown in pots filled with clay-loam soil enhanced with leaf compost and fertilizer. The experiment was performed as a factorial (two factors including genotypes and drought treatments) in a completely randomized design with three replicates in the greenhouse. The average temperature of the greenhouse was 35/26°C day/night during 1 3 the experiment. The soil moisture (VWC) of each pot was measured using the ProCheck sensor with a 10HS moisture sensor (METER Group, Inc., USA) during the experiment. The seedlings were well-watered until 40 days after sowing to maintain soil moisture at optimum level (about 30% VWC). The seedlings were thinned twice (at 14 and 21 days after sowing) and five uniform seedlings were kept. Progressive drought stress was applied to plants by withdrawing irrigation, starting from 40 days after sowing (early flowering stage). Leaf samples were collected from plants before stress (C0) and when their soil moisture reached about 15% (D1), 12% (D2), 9% (D3), and 6% (D4), after 4, 6, 8, and 12 days of no watering, respectively.

RNA extraction and qRT-PCR analysis
Total RNA of leaf samples was extracted using TRIzol Reagent (Thermo Fisher Scientific) according to the manufacturer's instructions. The Residual DNA of the extract was digested using DNase I (Thermo Fisher Scientific). Thereafter, the cDNA library was synthesized from RNA extract using the RevertAid RT Reverse Transcription Kit (Thermo Fisher Scientific). The miRNA-specific cDNAs also were generated from RNA/RT primer mixture based on the stemloop method (Chen et al. 2005). The qRT-PCR was performed using the SYBR Green Master Mix (Thermo Fisher Scientific) on a CFX96 Real-Time PCR Detection System (Bio-Rad Laboratories, Inc., USA). The Actin97 gene of sesame was used as endogenous controls for further analysis. The primers for DRG and TF genes were designed using Primer3 software (Untergasser et al. 2012) and the miRNA primers were designed by sRNAPrimerDB (Xie et al. 2019) (Supplemental Table S1). The relative expression of each gene was calculated using 2 −ΔΔCt as described by Livak and Schmittgen (2001).

Gene regulatory networks
Six gene regulatory networks (GRNs) were constructed (three GRNs for each set) using predicted relations between studied genes, TFs, and miRNAs of sesame (Supplemental Table S4-9). The degree of nodes in all networks followed a power-law distribution which indicated that they were scalefree (0.7<R 2 <0.87). The general properties of networks are shown in Table 1. To identify key regulators, the degree and betweenness centrality were calculated for each network. Nodes with a high degree and betweenness centrality were considered hub genes (Supplemental Table S10-13).
MCODE algorithm was applied to the integrated miR-TF-DRG network of each set to detect key regulatory relations of sesame under drought stress (Supplemental Table S14-15). Seven modules (c1, c2, c3, c4, c5, c6, and c7) were identified in the DRC network, and five (d1, d2, d3, d4, and d5) in the DTR network (Fig. 2). Module c1 with 25 nodes and 38 edges was the largest detected cluster in the DRC network. Notably, genes such as DHN1 (105, Fig. 3a, the CER1 expression was up-regulated in both genotypes under low drought levels (D1 and D2). However, tolerant genotype showed a higher expression level of this gene, compared to sensitive genotype. At higher drought levels (D3 and D4), the expression of CER1 became down-regulated in sensitive genotype, while it was still upregulated in tolerant genotype. The TGA1 was up-regulated at all drought levels in the tolerant genotype and had a significantly higher expression than the sensitive samples (Fig. 3b). On the contrary, TGA1 was stable or down-regulated under drought stress in the sensitive genotype, except at the D1 level, which showed a slight increase in relative expression. The miR530 was down-regulated at all drought levels in tolerant samples and had a lower expression than in the sensitive samples (Fig. 3c. But in the sensitive genotype,  1 3 miR530 was only down-regulated under severe drought (D3 and D4). The expression of miR319 was up-regulated in the tolerant genotype under all drought levels, except at the D4 level (Fig. 3d). On the contrary, it was down-regulated at all drought levels in the sensitive samples and showed a lower relative expression than in the tolerant ones. In the sensitive samples, the NAC29L expression was up-regulated at all drought levels (Fig. 3e). On the other hand, the expression of NAC29L was down-regulated in the tolerant genotypes under severe drought conditions (D3 and D4). Overall, the relative expression of NAC29L was higher in the sensitive samples compared with the tolerant ones, at all drought levels. The SUC2L was up-regulated under all drought levels in the tolerant genotype, while it was mostly down-regulated in the sensitive genotype under drought conditions, except at the D1 level (Fig. 3f). The miR171 expression was down-regulated in both genotypes at almost all drought levels (Fig. 3g). However, the tolerant genotype noticeably showed a lower relative expression of miR171 than the sensitive genotype.

Discussion
The present study was conducted to explore regulatory networks that control sesame response to drought stress. To this end, firstly the drought-responsive genes were identified using a comparative analysis of transcriptome data of two sesame genotypes (one sensitive and one tolerant). To investigate the general and specific tolerance response of sesame while facing drought, two DRC and DTR gene sets were selected based on similarity and difference between highly DEGs of sensitive and tolerant genotypes, respectively. As shown in the results, genes from families like the late embryogenesis abundant (LEA), dehydrin (DHN), and heat shock proteins (HSP) were found up-regulated in the DRC set. Many studies confirmed the positive role of these gene families against oxidative stresses including drought (Augustine 2016;Magwanga et al. 2018;Yu et al. 2018). The alcohol dehydrogenase 1 (ADH1) was also detected in the DRC set. Shi et al. (2017) in a study on the function of ADH1 in Arabidopsis under abiotic stress showed that overexpression of ADH1 enhanced tolerance to abiotic stress including drought by conferring ABA sensitivity, increasing expression of stress-related genes, and accumulations of more protective osmolyte. Pectate lyase was down-regulated in Arabidopsis and Tomato during drought stress, similar to our results (Le Gall et al. 2015;Iovieno et al. 2016). Auxin-binding proteins (ABPs) have also been found down-regulated in maize and potato during drought stress (Moon et al. 2018;Wang et al. 2019a). Cytokinin dehydrogenases (CKXs) are flavoenzymes that catalyze the oxidation of cytokinins. Overexpression of CKX resulted in high drought tolerance by slowed growth, enhanced root system, and reduced shoot structures (Prerostova et al. 2018). CKX1-like was down-regulated in all samples (presented in the DRC gene set) but CKX7 has been found up-regulated in the DTR gene set. After all, cytokinin seems to play an important role in response to drought conditions.
In DTR set up-regulated genes like flavonoid 3'-monooxygenase-like and 2-oxoglutarate-dependent dioxygenase AOP1 were present. Flavonoid 3'-monooxygenase also known as flavonoid 3'-hydroxylase (F3H), is a member of the 2-oxoglutarate-dependent dioxygenase family that has a key role in flavonoids biosynthesis. Flavonoids are plant secondary metabolites involved in various biological functions including response to abiotic stresses. Many reports in various plants have shown that F3H improved drought tolerance (Liu et al. 2013;Han et al. 2017). The AOPs are a subfamily of 2-oxoglutarate-dependent dioxygenases involved in glucosinolates (secondary metabolites) biosynthesis. It was suggested that glucosinolate accumulation has a role in drought tolerance through controlling stomatal closure and preventing water losses (Eom et al. 2018). The casparian strip membrane protein (CASP) was found among the top up-regulated genes in the DTR set. CASPs mediate casparian strip formation by localizing lignin deposition sites (Roppolo et al. 2014). It has been suggested that lignification and casparian strip formation might lead to drought tolerance by preventing toxic compounds uptake and water loss (Chen et al. 2011;Robbins et al. 2014). Overall, previous studies support most of our selected genes as drought-related.
GRNs were constructed based on the predicted interactions between drought-related genes and regulators. TFs and miRNAs with high centrality measures including degree and betweenness were considered as the key regulators. Degree centrality is the simplest measure of gene connectivity in the GRNs that is calculated by the number of direct connections of that gene with other genes. Regulator genes with the highest degrees control a large number of drought-responsive genes and therefore could play an important role in regulating plant response to stress and inducing tolerance. Betweenness was measured by the number of shortest paths that pass through each node. Genes with high betweenness are crucial regulators in GRNs that act as bridges between other genes and connect them (Koschützki and Schreiber 2008;Sircar and Parekh 2019) identified important genes based on centrality measures in a network-based study on Oryza sativa genotypes under drought stress.
Based on results, the transcription factor IIIA (TFIIIA) and ZAT5 from the Cys2/His2 zinc finger family were hub genes in the DRC-TF network. Overall, C2H2 zinc finger proteins improve plant tolerance during drought stress by increasing the level of ABA, proline, soluble sugars, and ROS scavenging enzymes like superoxide dismutase and peroxidase through ABA-dependent and ABA-independent pathways, as described by Wang et al. (2019b). MYBrelated protein 306 and some other MYB TFs were found in the DRC-TF network (Table 2) as hub genes belonging to the R2R3-MYB subfamily. The R2R3-MYB TFs regulate important plant biological processes such as primary and secondary metabolism, cell fate and identity, developmental processes, and responses to biotic and abiotic stresses (Dubos et al. 2010;Tang et al. 2019) showed that overexpression of OsMYB6, an R2R3-MYB TF has enhanced drought tolerance. Another identified TF with high centrality measures in the DRC-TF network was dehydration-responsive element-binding 2D-like (DREB2D-like), a member of the DREB family of the AP2/EREBF superfamily. Ansari et al. (2017) in a study on muskmelon genotypes (Cucumis melo L.) showed that high up-regulation of DREB2D improved tolerance of tolerant genotype under severe drought stress.
On the other hand, the PIF1 and PIF3 were among the DTR-TF network hub genes. Phytochrome-interacting factors (PIFs) are members of the basic helix-loop-helix (bHLH) transcription factors family. It was shown that ZmPIF1 enhances drought tolerance and prevents water loss by reducing the stomatal opening in transgenic rice plants (Gao et al. 2018). ZmPIF3 could also improve drought tolerance of transgenic rice plants by positive regulation of Fig. 3 Relative expression of 7 selected transcripts, including a) CER1, b) TGA1, c) miR530, d) miR319, e) NAC29L, f) SUC2L, and g) miR171, in the sesame tolerant and sensitive genotypes under drought conditions. Samples were collected before the stress (C0) and when soil moisture reached to 15% (D1), 12% (D2), 9% (D3) and 6% (D4), after 4, 6, 8, and 12 days of no watering, respectively ◂ 1 3 relative water content, chlorophyll content, fluorescence, and drought-responsive genes expression (Gao et al. 2015). Knotted1-like-homeobox 1 (KNOX1) was another hub TF found in the DTR-TF network. It was suggested that KNOX1 could be used to improve crops' drought resistance (Townsley et al. 2013). CPRF2 and basic leucine zipper 43 (bZIP43), members of plant bZIP TF family also were seen among DTR-TF network hub genes. Many members of the bZIP family were reported to positively regulate tolerance to abiotic stresses including drought (Banerjee and Roychoudhury 2017;Wang et al. 2018a;Yang et al. 2019b) showed that overexpression of OsbZIP62 improved drought tolerance by regulating the expression of stress-related genes.
As previously shown, in both miR-DRG networks miR9773, miR5658 and miR156f had the highest centrality measures. Studies on miR9773 and miR5658 and their regulatory function during abiotic stresses including drought have just begun. The expression level of ta-miR9773 was significantly increased in Einkorn wheat under drought stress (Baloglu et al. 2017). It was suggested that sly-miR5658 regulates drought stress response by controlling hormone signal transduction pathway genes in both sensitive and tolerant genotypes of tomato (Candar-Cakir et al. 2016). Studies on miR156 also demonstrated its regulatory roles in response to drought in alfalfa (Medicago sativa) (Arshad et al. 2017). The miR396b and miR172d were hub genes in the DRC-miR network. miR396 was identified as a key regulator for drought tolerance in many plants including rice (Nadarajah and Kumar 2019). The Gma-miR396 family was found as a positive (leaves) and negative (roots) regulator under low water availability stress in Arabidopsis plants . miR172 family has been reported to be involved in drought stress response in potato (Solanum tuberosum) and Echinacea purpurea L. (Hwang et al. 2011a;Moradi and Khalili 2018). In the DTR-miR network, miR169 and miR395 were identified as hub miRNAs. Several studies demonstrated that miR169 plays a significant role in drought stress responses in various plants (Zhang et al. 2011;Ni et al. 2013). miR169 family has been shown down-regulated in leaves but upregulated in roots of wheat under drought stress (Akdogan et al. 2016). miR395 also has an important role in response to drought stress and was found up-regulated in several plants under drought stress (Li et al. 2013a;Akdogan et al. 2016).
In genetic studies, modules are frequently used to identify important regulatory relationships between regulators and putative genes (Saelens et al. 2018). Therefore, the clustering algorithm was applied to the miR-TF-DRG networks to integrate and improve the understanding of the drought-regulatory mechanisms. Genes like DHNs, ABP19a, and pectate lyases presented in DRC modules were previously shown as the drought-responsive gene. Sucrose synthase (SUS2) was also found in DRC modules which were identified as drought-responsive genes (Yang et al. 2019a). MYB63, TF found in DRC modules is a transcriptional activator of lignin biosynthetic genes (Zhou et al. 2009). Lignin has a positive role in drought tolerance ) and due to this fact, MYB63 is possibly involved in plant drought responses. DREB 2D-like, MYB-related 306, PIF1, and CPRF are also drought-related TFs as previously discussed. The bZIP16 is another transcription factor found in DRC modules that was reported to regulate drought resistance (Chen et al. 2012). As shown earlier, miR169 and miR156 have a role in regulating drought responses. It was also shown that miR399 regulates plant responses to salt, ABA, and drought (Baek et al. 2016). On the other hand, CKX7 was presented in DTR modules for which its relations with drought stress previously have been described. CYP82D47-like which is a member of the CYP82 family also has been found in DTR modules. The CYP82 of Soybean has been reported to enhance tolerance to salinity and drought stresses (Yan et al. 2016). Among TFs presented in DTR modules, members of the MYB, WRKY, and NAC families can be seen. MYB108 has been found to regulate responses to biotic and abiotic stresses including drought (Baldoni et al. 2015;Yang et al. 2019c). WRKY30 has been shown to positively regulate drought tolerance (Zhu et al. 2018;El-Esawi et al. 2019). miRNAs like miR171 and miR529 were found in DTR modules. miR171 family was identified as drought-responsive microRNAs which regulate plant response during the stress (Zhou et al. 2010;Hwang et al. 2011b;Zhang et al. 2017). It was shown that miR529 is also involved in drought stress responses (Bakhshi et al. 2016;Zhang et al. 2017). The relations of miR169 and miR395 with drought responses have been previously discussed. Overall, previous researchers have reported the effectiveness of using network-based studies in identifying regulatory genes and relationships under abiotic stress conditions (Song et al. 2017;Ko and Brandizzi 2020).
In the following, the expression of seven transcripts from module d1 was examined to validate identified regulatory relationships under drought stress. As shown in the results, the expression of miR530 in the tolerant genotype was upregulated and was lower than in the sensitive genotype, under all drought levels. Downregulation of miR530 was reported in some other plants under water stress conditions (Balyan et al. 2017;Ning et al. 2017;Hu et al. 2021). Similarly, the miR530 was downregulated in the flag-leaves of tolerant rice cultivar but was upregulated in the sensitive rice cultivar under drought stress (Balyan et al. 2017). According to the results, the miR530 targets both TGA1 and CER1.
The expression of both genes in the tolerant genotype was upregulated and was higher than in the sensitive genotype during the drought stress. Induced TGA1 and CER1 expression are well compatible with the downregulation of their inhibitor, miR530, in the tolerant genotype under drought conditions. The TGA1 is a member of TGA (TGACG motif-binding factor) transcription factors, a sub-family of the basic leucine zipper (bZIP) gene family. Recent studies revealed that TGA genes are involved in regulating abiotic stress responses, especially drought tolerance Chen et al. 2021). Moreover, the result of the present study also showed that TGA1 is one of the regulators of CER1 expression. The CER1, a very-long-chain aldehyde (VLCA) decarbonylase, is a core component of the very-long-chain alkane synthesis complex that catalyzes VLCAs synthesis from very-longchain-fatty-acyl-CoA (VLCFA). VLCAs are one of the main components of cuticular wax, which makes up about 50-80 of the total wax content (Bourdenx et al. 2011). Recent studies have well demonstrated the important role of cuticular waxes in enhancing drought tolerance (Lee and Suh 2013;Xue et al. 2017). Overexpression of CER1 in various plants resulted in enhanced drought tolerance by increasing cuticular wax formation (Bourdenx et al. 2011;Wang et al. 2015bWang et al. , 2020. Therefore, plants with a higher potential for induction of CER1 under drought stress are likely to have a higher tolerance, as the sesame tolerant genotype in this study. Overall, it seems the tolerant genotype directly and indirectly (by TGA1 induction) induced upregulation of CER1 through inhibition of miR530 expression during drought stress. Regulatory activity of this module (miR530/TGA1/CER) increased wax accumulation and enhanced tolerance under drought stress.
On the other hand, miR319 expression was upregulated in the tolerant genotype while was conversely downregulated in the sensitive, under drought stress. miR319 is one of the oldest and most conserved miRNA families in plants. The overexpression of the rice miR319 in bentgrass significantly enhanced drought tolerance (Zhou et al. 2013;Zhou and Luo 2014). Based on our results, miR319 regulates NAC29L expression. Despite the increased expression at low levels (D1), the lower NAC29L expression and its decreasing trend at higher drought levels were consistent with the upregulation of miR319 (NAC29L inhibitor) in the tolerant genotype compared to the sensitive. The NAC TF family plays a key role in the leaf senescence process (Podzimska-Sroka et al. 2015). The leaf senescence is considered a type of programmed cell death that could be inappropriately induced under drought stress and disrupts plants' adaptation mechanisms (Rivero et al. 2007). It was shown in Arabidopsis that NAC29 mediates inhibition of stomatal closure, increase in water loss, and consequent initiation of the leaf senescence process by activating its target gene, SAG113 (Zhang and Gan 2012). NAC29 also induces chlorophyll degradation in Arabidopsis leaves by enhancing the AAO3 expression (Yang et al. 2014). Such a senescence-related mechanism mediated by NAC genes imposes adverse effects on the plants under drought conditions. Therefore, inhibition of NAC29L expression might be intended to prevent or delay a similar process in sesame tolerant genotype during drought stress.
Moreover, the NAC29L was identified as a potential regulator of SUC2L, in this study. Therefore, the expression of SUC2L was expected to be downregulated in the tolerant genotype under drought stress, due to the inhibition of NAC29L. However, contrary to this expectation, the expression of SUC2L was upregulated in the tolerant genotype at all drought levels. While, it was conversely downregulated in the sensitive genotype, almost at all drought levels (except at the D1 level). But, on the other hand, miR171, another predicted regulator of SUC2L and NAC29L, was downregulated in the tolerant genotype during the drought. Although miR171 was also downregulated in the sensitive samples, its relative expression was significantly higher than in the tolerant genotype at all drought levels. Thus, the upregulation of SUC2L in the tolerant genotype could be explained with higher inhibition of miR171 expression, compared with the sensitive genotype during drought stress. It seems the regulatory activity of miR171 on SUC2L expression was more efficient than NAC27L. Although, SUC2L expression might be enhanced by other regulators that were not identified in this study.
Anyway, SUC2 (also named SUT1) is a type I member of the sucrose transporter (SUT) family. It is involved in the sucrose transportation into the cell through a symport system and is required for apoplastic phloem sucrose loading (Gottwald et al. 2000;Xu et al. 2018). Besides, the sucrose transporters (SUTs) are also responsible for the soluble sugars accumulation in cells . Previous studies suggested that SUC2 activity is required for abiotic stress tolerance including drought (Gong et al. 2015;Du et al. 2020). In Arabidopsis and soybean, upregulation of the SUC2 gene enhanced sucrose transportation from leaves to roots and consequently maintained root growth under drought stress (Durand et al. 2016;Du et al. 2020). The sesame tolerant genotype might also have employed such a mechanism to tolerate drought through upregulation of SUC2L. Overall, qPCR results, in line with the previous studies, confirm the efficiency of module detection algorithms in identifying key drought regulatory relations.

Conclusions
In this study, we identified a total of 117 TFs and 133 miR-NAs that might be related to drought stress by the combination of in-silico data from transcriptome analysis and bioinformatics tools. The regulatory gene networks were constructed and possible key regulators of sesame drought response were explored using network analysis. Important regulatory relations between miRNAs, TFs, and DRGs were 1 3 identified using modules detection analysis. The TFs like TFIIIA, MYB 306, DREB2D-like, PIF1, PIF3, KNOX1, and CPRF2 were found as important drought-related TFs. miR-NAs like miR9773, miR5658, miR156, miR396, miR395, miR172, and miRNA169 also were among key drought regulators. Moreover, seven potentially drought tolerancerelated transcripts of sesame, from identified modules, were examined under drought conditions through a comparatives qRT-PCR analysis. The studied genes showed a different expression pattern in the tolerant genotype compared to the sensitive. It seems that miR530/TGA1/CER1, miR319/ NAC29L, and miR171/SUC2L modules play curial roles in sesame drought tolerance through regulating wax biosynthesis, leaf senescence, and sugar transport during drought, respectively (Fig. 4). Moreover, some detected genes that have not yet been reported to function in drought stress may also play a vital role in regulating sesame drought response. These findings could be used in the next studies on sesame and future breeding programs to improve drought tolerance.

Data availability
In-silico data used and analyzed during the current study are available in the NCBI BioProject under accession numbers PRJNA478474 (https:// www. ncbi. nlm. nih. gov/ biopr oject/ PRJNA 478474). Other data used to support the findings of this study are available in this manuscript and supporting materials.
Code availability Not applicable.

Conflicts of Interest/Competing Interests
The authors declare that they have no competing interests.
Ethics Approval Not applicable.