Transcription Factors and microRNA Genes Regulatory Network Construction Under Drought Stress in Sesame (Sesamum indicum L.)

Background: Drought is one of the most common environmental stresses affecting crops yield and quality. Sesame is an important oilseed crop that most likely faces drought during its growth due to growing in semi-arid and arid areas. Plants responses to drought controlled by regulatory mechanisms. Despite this importance, there is little information about Sesame regulatory mechanisms against drought stress. Results: 458 drought-related genes were identied using comprehensive RNA-seq data analysis of two susceptible and tolerant sesame genotypes under drought stress. These drought-responsive genes were included secondary metabolites biosynthesis-related Like F3H, sucrose biosynthesis-related like SUS2, transporters like SUC2, and protectives like LEA and HSP families. Interactions between identied genes and regulators including TFs and miRNAs were predicted using bioinformatics tools and related regulatory gene networks were constructed. Key regulators and relations of Sesame under drought stress were detected by network analysis. TFs belonged to DREB (DREB2D), MYB (MYB63), ZFP (TFIIIA), bZIP (bZIP16), bHLH (PIF1), WRKY (WRKY30) and NAC (NAC29) families were found among key regulators. mRNAs like miR399, miR169, miR156, miR5685, miR529, miR395, miR396, and miR172 also found as key drought regulators. Furthermore, a total of 117 TFs and 133 miRNAs that might be involved in drought stress were identied with this approach. Conclusions: Most of the identied TFs and almost all of the miRNAs are introduced for the rst time as potential regulators of drought response in Sesame. These regulators accompany with identied drought-related genes could be valuable candidates for future studies and breeding programs on Sesame under drought stress.

potential of Sesame production and adversely affects its growth, development, metabolism, and yield [7]. Therefore, improving drought-tolerant 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 identi ed and studied. The plant responds 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 binding 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 identi ed and proved that play key roles in drought tolerance [8,9]. In Arabidopsis, it has been shown that the expression of maize ZmMYB3R improves drought and salinity tolerance [10]. The expression of Soybean GmNAC085 transcription factor in Arabidopsis signi cantly led to improved drought tolerance of transgenic plants by positively regulating growth rate and reducing transpiration and cell membrane damage [11]. It was reported in rice that overexpression of OsNAC5, 6, 7, and 9 genes result in enlarged root and enhanced drought tolerance [12]. The SlWRKY58 and SlWRKY72 genes have been shown to be associated with improved drought tolerance in tomatoes [13]. Dossa et al. [14] 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 stress. At this level we can refer to microRNAs (miRNAs), non-coding single-stranded RNAs that regulate gene expression using sequence-speci c cleavage or translation inhibition of the target transcripts [15,16]. Plant miRNAs derive from their precursors that are transcripts of exon and intron of coding regions or intergenic regions [17,18,19]. After multiple processing steps, the mature miRNA forms from the miRNA gene. Mature miRNA then loaded onto the argonaute protein and along with some other proteins forms the RISC complex [17]. After all, miRNA by base pairing with mRNA would guide RISC to cleave the target or repress its translation [15,16]. 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 [20,21]. A study on maize showed that miRNAs such as miR474, miR168, miR528, and miR167 might be involved in regulating tolerance mechanisms against drought stress [22].
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 [23].
The development of high-throughput technologies in recent years and their combination with bioinformatics analysis enabled the researchers to have a clearer vision of drought response mechanisms in plants and helped them to identify the drought tolerance candidate genes more e ciently for use in the breeding program. Recently, studies on the role of transcription factors in Sesame responses to drought stress were increased and have been growing [14,24,25,26]. However, our knowledge about many of Sesame TFs and their roles under drought stress is limited and insu cient. Meanwhile, most of these researches considered TFs response separately, and collective and networkbased studies that focused on the interactions between TFs and stress-related genes have 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 Sesame drought-related genes by performing a comparative RNA-seq-based analysis and then construct gene regulatory networks (GRNs) using the predicted relations of identi ed genes with Sesame TFs and miRNAs. Finally, hub genes were identi ed by network analysis as potentially key regulators in drought stress to be veri ed and used in future breeding programs.

Enrichment Analysis
To determine the functional properties of each set, genes were annotated by Gene Ontology (GO) terms and Kyoto Encyclopedia Genes and Genomes (KEGG) pathways. Cellular parts such as nuclear chromatin, apoplast, nuclear chromosome part, cell wall, and some others were signi cantly enriched for DRC genes (Fig. 3-a). They were also signi cantly enriched for the response to light intensity, hydrogen peroxide, heat, toxic substance, and some others from the biological process category (Fig. 3-b).

Gene Regulatory Networks
Three gene regulatory networks (GRNs) were constructed for each set using predicted relations between studied genes, TFs, and miRNAs of Sesame (Additional le 1: Table S3-8). The degree of nodes in all networks followed a power-law distribution which indicated that they were scale-free (0.7<R 2 <0.87).
General properties of networks are shown in Table 1. To identify key regulators, the degree and betweenness centrality were calculated for each network. Degree centrality is the simplest measure of gene connectivity in the GRNs that 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 [27]. Thus, nodes with high degree and betweenness centrality considered as hub genes (Additional le 1: S9-12).

Discussion
The present study conducted to explore regulatory networks that control Sesame respond to drought stress. To this end, the drought-responsive genes rstly were identi ed using a comparative analysis of transcriptome data of two susceptible and tolerant Sesame genotypes. To investigate the general and speci c tolerance response of Sesame in the facing with drought stress, two DRC and DTR gene sets were selected based on similarity and difference between susceptible and tolerant genotypes DEGs, respectively. As shown in results, genes from families like LEA, Dehydrin, and HSP were found upregulated in the DRC set. Many studies con rmed the positive role of these gene families against oxidative stresses including drought [28, 29,30]. Other genes like alcohol dehydrogenase 1 (ADH1) and oleosin S1-2 were also detected in the DRC set. Shi et al. [31] in a study on 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 (such as sugars and sucrose) [31]. Little is known about the role of oleosin during drought stress. Oleosins provide stability to oil bodies and prevent their fusion to an enlarged abnormal oil-structure. The irregularly enlarged oil bodies caused the deformation of the cell nucleus and disrupt its function [32,33]. Oleosin conferred freezing tolerance to Arabidopsis seeds through this mechanism [32]. Oleosins may likely stabilize lipid bodies during the dehydration under drought conditions with a similar function. Genes such as pectate lyase 5 and auxin-binding protein 19a (ABP19a) were found down-regulated among DRC genes. Pectate lyase was down-regulated in Arabidopsis and Tomato during drought stress, similar to our results [34,35]. Pectat lyase is involved in cell wall modi cation and formation of root structures [34,36]. Auxin-binding proteins have also been found down-regulated in maize and potato during drought stress [37,38]. Similarly, ABP was reported to have a role in cell wall modi cation and cell expansion [39]. It seems the genes related to cell wall modi cation and degradation were suppressed by drought stress in both Sesame genotypes. Cytokinin dehydrogenases (CKXs) are avoenzymes that catalyze the oxidation of cytokinins. Overexpression of CKX resulted in high drought tolerance by slowed growth, enhanced root system, and reduced shoot structures. However, these features such as slow growth and decreased leaves and shoots were present even under normal condition which are undesirable [40,41]. CKX1-like was down-regulated in all samples (presented in 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 avonoid 3'-monooxygenase-like and 2-oxoglutarate-dependent dioxygenase AOP1 were present. Flavonoid 3'-monooxygenase also known as avonoid 3'-hydroxylase (F3H), is a member of 2-oxoglutarate-dependent dioxygenase family that has a key role in avonoids biosynthesis. Flavonoids are plant secondary metabolites which involved in various biological functions including response to abiotic stresses. Many reports in various plants have been showed that F3H improved drought tolerance [42,43]. The AOPs are a subfamily of 2-oxoglutarate-dependent dioxygenases which 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 loses [44]. Secondary metabolites and oxoglutarate-dependent dioxygenase family seem to play important roles in drought tolerance, as both AOP1 and F3H are members of them. The casparian strip membrane protein (CASP) and extracellular ribonuclease LE-like (RNase LE-like) also found among top up-regulated genes in the DTR set. CASPs are transmembrane proteins that mediate casparian strips formation by localizing lignin deposition sites [45]. It has been suggested that ligni cation and casparian strip formation might lead to drought tolerance by preventing toxic compounds uptake and water loss [46, 47,48]. RNase LE is an extracellular class I RNases [49]. The class I RNases has been reported that induced by Pi de ciency, senescence, wounding and osmotic stresses including drought [50]. RNase LE has a phosphate-scavenging function in response to P i de ciency by hydrolyzing extracellular RNAs [51,52]. One of the consequences of drought stress in plants is the decrease in concentrations of nutrients especially P which caused by the decrease in P uptake of roots [53]. Thus, RNase LE up-regulation might contribute to the drought stress tolerance using this function.
The sucrose transport protein (SUC2-like) was also up-regulated in the DTR set. It was shown in Arabidopsis that SUC2 up-regulated during water de cit and increased C export to the roots by sucrose phloem loading [54]. Genes such as NIM1-interacting 2-like, glutamate receptor 2.9-like, receptor-like protein, and pathogenesis-related protein 5-like were found down-regulated in the DTR set. Although there have been reports that overexpression some of these genes have a positive effect on tolerance to abiotic stresses [55,56,57], but better-known function of all of them is resistance to disease and biotic stresses [58, 59, 60, 61, 62, 63]. Based on the results, it seems likely that tolerant genotype at high drought levels down-regulate the expression of some of the disease-related genes (not all) and thus spend their stored energy on producing other vital compounds to survive under severe drought stress. Overall, previous studies support most of our selected genes as drought-related.
Interactions between identi ed drought-related genes and gene regulators including TFs and miRNAs were predicted. Based on results, the transcription factor IIIA (TFIIIA) and ZAT5 from the Cys2/His2 zinc nger family were hub genes in the DRC-TF network. Overall, C2H2 zinc nger proteins improve plant tolerance during drought stress by increasing level of ABA, proline, soluble sugars and ROS scavenging enzymes like superoxide dismutase (SOD) and peroxidase (POD) through ABA-dependent and ABAindependent pathways, as described by Wang  On the other hand, the PIF1 and PIF3 were among DTR-TF network hub genes. Phytochrome-interacting factors (PIFs) are members of basic helix-loop-helix (bHLH) transcription factors family. ZmPIF1 was shown that enhances drought tolerance and prevents water loss by reducing the stomatal opening in transgenic rice plants [71]. ZmPIF3 could also improve drought tolerance of transgenic rice plants by positive regulation of relative water content, chlorophyll content, uorescence, and drought-responsive genes expression [72]. Knotted1-like-homeobox 1 (KNOX1) was another hub TF that found in the DTR-TF network. It was suggested that KNOX1 could be used to improve crops drought resistance [73]. 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 that positively regulate tolerance to abiotic stresses including drought [74,75]. Yang et al. [76] showed that overexpression of OsbZIP62 improved drought tolerance by regulating the expression of stress-related genes [76].
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 signi cantly increased in Einkorn wheat under drought stress [77]. 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 [78]. Studies on miR156 also demonstrated its regulatory roles in response to drought and salinity stresses in alfalfa (Medicago sativa) [79,80]. The miR396b and miR172d were hub genes in the DRC-miR network. miR396 identi ed as a key regulator for drought tolerance in many plants including rice [81]. The gma-miR396 family was found as positive (leaves) and negative (roots) regulator under low water availability stress in Arabidopsis plants [82]. miR172 family has been reported that involved in drought stress response in potato (Solanum tuberosum) and Echinacea purpurea L. [83,84]. In the DTR-miR network, miR169 and miR395 were identi ed as hub miRNAs. Several studies demonstrated that miR169 plays a signi cant role in drought stress responses in various plants [85,86,87]. miR169 family has been shown down-regulated in leaves but upregulated in roots of wheat under drought stress [88]. miR395 also has an important role in response to drought stress and was found upregulated in several plants under drought stress [88,89,90].
Important regulatory interactions between DRGs and regulators (TFs and miRNAs) in the miR-TF-Gene networks were identi ed using the MCODE clustering algorithm. 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 identi ed as drought-responsive genes [91,92]. MYB63, TF found in DRC modules is a transcriptional activator of lignin biosynthetic genes [93]. Lignin has a positive role in drought tolerance [94] 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 another transcription factor found in DRC modules was reported to regulate drought resistance [95]. As shown earlier, miR169 and miR156 have a role in regulating drought responses. miR399 was also shown that regulates plant responses to salt, ABA, and drought [96,97]. On the other hand, SUC2-like and CKX7 genes were presented in DTR modules which their 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 that enhances tolerance to salinity and drought stresses [98]. Among TFs presented in DTR modules, members of MYB, WRKY and NAC family can be seen. MYB108 has been found to regulate responses to biotic and abiotic stresses including drought [99,100,101]. WRKY30 has been shown that positively regulate drought tolerance [102,103]. NAC29 also enhanced salt and drought tolerance in transgenic Arabidopsis [104]. miRNAs like miR171 and miR529 were found in DTR modules. miR171 family was identi ed as droughtresponsive microRNAs which regulate plant response during the stress [105,106,107]. miR529 also has been shown that involved in drought stress responses [105,107,108]. The relation of miR169 and miR395 with drought responses have been previously discussed. Given that these genes are related to drought stress, identi ed relations and other genes of modules might be involved in Sesame drought stress responses.

Conclusion
In this study, we identi ed a total of 458 drought-related genes (290 DRC and 168 DTR) using RNA-seq data analysis of two susceptible and tolerant sesame genotypes under drought stress. Drought-related genes like LEA protein Dc3, DHN1-like, ADH1, HSP, CKX7, F3H, SUC2, SUS2 were identi ed. Then, their relation with gene regulators (TFs and miRNAs) were predicted using bioinformatics tools. A total of 117 TFs and 133 miRNAs that might be related to drought stress were identi ed. Using predicted relations between DRGs and TFs or miRNAs the regulatory gene networks were constructed. Possible key regulators of Sesame drought response were identi ed using centrality measures. The TF and miRNA regulatory networks were merged and miR-TF-DRG networks were constructed. Key regulatory relations between miRNAs, TFs and DRGs were identi ed using modules detecting. The TFs like TFIIIA, MYBrelated protein 306, DREB2D-like, PIF1, PIF3, KNOX1, and CPRF2 found as hub TFs. miRNAs like miR9773, miR5658, miR156, miR396, miR395, miR172, and miRNA169 were among hub genes. Moreover, TFs including MYB63, MYB108, bZIP16, WRKY30, and NAC29 were detected in modules of DRG-TF-miR networks. miR171, miR399, miR171, and miR529 also were identi ed in modules. Most of the identi ed TFs and almost all miRNAs of this study have not yet been reported for Sesame drought responses. These ndings could be used in the next studies on Sesame to validate their function in drought stress. Finally, the validated genes and regulators could be used in Sesame future breeding programs to improve drought tolerance.

Data Processing
The quality of the raw data was checked using FASTQC [110]. Then, excluding low-quality nucleotides and reads and also adaptor removing were done using Trimmomatic [111]. After ltering 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 11560 using HISAT2 [112]. Count matrix was calculated based on matched reads using HTSeq [113]. Count data were normalized using the TMM method included in the edgeR package [114] under the R platform [115]. Differentially expressed genes (DEGs) analysis was performed on normalized data using the edgeR package [114]. DEGs with FDR<0.01 and |logFC|>3 were retained for each sample. Genes that differently expressed in all samples were selected as drought-responsive core (DRC) gene set. On the other hand, genes that exclusively differently expressed in tolerant samples under severe drought (T4 and T4T3) compared to sensitive samples at the same condition, as well as other tolerant samples were also selected as drought tolerance-related (DTR) gene set. The analysis was then performed separately on these two sets of genes.

Enrichment Analysis
Each of the sets was functionally annotated using Gene Ontology (GO) main categories including cellular components (CC), molecular functions (MF) and biological processes (BP). GO terms for each gene obtained from their Arabidopsis homologues that were identi ed by blastp of BLAST+ software [116]. Target genes were also assigned to biological pathways using the Kyoto Encyclopedia Genes and Genomes (KEGG) pathway database [117]. Enrichment analysis was done using clusterPro ler under the R platform [118].

TF-Gene Network Construction
Sesame TFs that targeted gene sets were identi ed using TF Enrichment tool of PlantRegMap database [119]. TFs that possess a signi cantly over-presented number of targets in gene sets based on Fisher's exact tests were selected (adjusted p-value cut-off= 0.05). TF-Gene networks were constructed using regulatory relations between selected TFs and genes of each set.

miR-Gene Network Construction
A list of mature Sesame miRNA sequences was obtained from the Sinbase database [120] and previous studies [121,122]. miRNA that targeted gene sets were predicted by psRNATarget web-based tool (expectation-value cut-off= 5) [123,124]. Using predicted miRNA-target gene interactions miR-Gene regulatory network was constructed for each set.

miR-TF-Gene Network Construction
The bidirectional interactions of identi ed miRNAs and TFs with each other were predicted using the Regulation Prediction tool of PlantRegMap database (TFs target miRNAs) and psRNATarget website (miRNAs target TFs). The integrated miR-TF-Gene GRNs were constructed by merging TF-Gene, miR-Gene and miR-TF networks of each set.

Network visualization and Analysis
Networks were visualized by Cytoscape software [125,126]. Centrality measures including degree and betweenness centrality were calculated for GRNs using CytoNCA, a Cytoscape plug-in [127]. Availability of data and materials All RNA sequencing data used and analysed during the current study are available in the NCBI BioProject under accession numbers PRJNA478474 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA478474).

Competing interests
The authors declare that they have no competing interests.

Funding
There was no funding for this research.  Table S1 Relative expression of top 15 Up and Down-regulated genes of DRC set under drought stress in all samples (selected based T4 samples). Samples collected at 16%, 13%, 10%, and 8% soil moisture named T1, T2, T3, and T4 for tolerant genotype and S1, S2, S3, and S4 for sensitive genotype, respectively.           Table S13 General details of detected modules in the miR-TF-DRC network.