Differentially Expression Analysis
To identify the general mechanism of drought response in Sesame, 290 differentially expressed genes (DEGs) that were shared between all genotypes and samples were selected as drought-responsive core (DRC) gene-set (Fig. 1 and Additional file 1: Table S1). Protein coding genes like ‘LEA protein Dc3’ (105157537), ‘dehydrin DHN1-like’ (105179059), ‘heat shock protein’ (105164479), ‘alcohol dehydrogenase 1’ (105159420), ‘oleosin S1-2’ (105175017) and some ncRNA genes (110012740, 105162944) were up-regulated, from this set. On the contrary, expressions of protein-coding genes such as ‘auxin-binding protein 19a’ (105157802), ‘pectate lyase 5’ (105168345), ‘cytokinin dehydrogenase 1-like’ (105162564) and ncRNAs (105173052, 105173721) were recorded as down-regulated. Most of the DRC genes had similar expression levels in both genotypes under low drought stress, but respond to higher drought intensities was different between tolerant and susceptible samples. Notably ‘SPX domain-containing protein 3-like’ (105173318), ‘monogalactosyldiacylglycerol synthase 2’ (105166432), ‘pathogenesis-related protein 1-like’ (105162859), ‘expansin-A10’ (105159984) and some uncharacterized ncRNA and protein-coding genes (110012421, 105163110, 105162944, 105168585) showed higher expression compared to sensitive samples.
On the other hand, another 168 DEGs which only differentially expressed in tolerant samples at high drought levels were selected as drought tolerance-related (DTR) gene-set to understand the specific mechanism of Sesame drought tolerance (Fig. 2 and Additional file 1: Table S2). Up-regulated genes including ‘flavonoid 3'-monooxygenase-like’ (105175177), ‘casparian strip membrane protein’ (105177354), ‘2-oxoglutarate-dependent dioxygenase AOP1’ (105157562), ‘extracellular ribonuclease LE-like’ (105156533) and ‘sucrose transport protein SUC2-like’ (105176204) were found in this set. Genes like ‘NIM1-interacting 2-like’ (105156738), ‘glutamate receptor 2.9-like’ (105155884), ‘receptor-like protein 12’ (105163973) and ‘pathogenesis-related protein 5-like’ (105166178) were also down-regulated.
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 significantly enriched for DRC genes (Fig. 3-a). They were also significantly enriched for the response to light intensity, hydrogen peroxide, heat, toxic substance, and some others from the biological process category (Fig. 3-b). Molecular functions such as DNA-dependent ATPase and galactosyltransferase activities were significantly enriched for the DRC set (Fig. 3-c). KEGG pathways including ‘protein processing’, ‘DNA replication’, ‘pentose and glucuronate interconversions’, and ‘galactose metabolism’ were significantly annotated for them (Fig. 5-a). On DTR genes biological processes including inorganic anion, anion, and nitrate transport, and response to nitrogen compound were significantly enriched (Fig. 4-a). From molecular function aspect, various transmembrane transporter (inorganic anion, anion, and secondary active), and symporter (carbohydrate:proton, carbohydrate:proton and solute:proton) activities were significant, on DTR set (Fig. 4-b). Moreover, ‘Phenylpropanoid biosynthesis’, ‘starch and sucrose’, ‘pyrimidine metabolism’, and ‘cyanoamino acid metabolism’ pathways were also significantly enriched based on the KEGG database (Fig. 5-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 file 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<R2<0.87). General properties of networks are shown in Table 1.
Table 1 General properties of gene regulatory networks.
|
Network
|
Nodes
|
Edges
|
DRGs
|
TFs
|
miRs
|
TF-DRC
|
289
|
1121
|
228
|
61
|
-
|
miR-DRC
|
375
|
947
|
247
|
-
|
128
|
TF-DTR
|
116
|
235
|
76
|
40
|
-
|
miR-DTR
|
276
|
580
|
147
|
-
|
129
|
miR-TF-DRC
|
489
|
2659
|
281
|
80
|
128
|
miR-TF-DRC
|
349
|
1280
|
156
|
64
|
129
|
DRG: Drought-Related Gene, TF: Transcription Factor Gene, miR: miRNA Gene.
|
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 file 1: S9-12).
TFs including ‘transcription factor IIIA’ (105155794), ‘MYB-related protein 306’ (105157040), ‘dehydration-responsive element-binding protein 2D-like’ (105166097), ‘zinc finger protein ZAT5’ (105161312), and some others were hub genes in TF-DRC network as shown in Table 2. Notably, most of the top hub TFs in this network were the members of MYB family. In the DTR-TF network, ‘transcription factor PIF1’ (105165703), ‘homeobox protein knotted-1-like 1’ (105168958), ‘light-inducible protein CPRF2’ (105172269) and ‘basic leucine zipper 43’ (105177566) were among top hub TFs (Table 3). TFs including, ‘light-inducible protein CPRF2’ (105172269), ‘transcription factor PIF3’ (105169008) and ‘transcription factor PIF1’ (105165703) were common between the two networks. However, they found as hub genes only in TF-DTR network. Commonly miR9773, miR5658 and miR156f were top hub miRNAs in both miR networks (Table 4). The majority of miRNAs were shared between two miR networks but there were some different hub miRNAs between two networks, like miR396b, miR172d and miR8140 in miR-DRC and miR169i and miR169a and miR395 in miR-DTR.
Table 2 Top 10 hub TFs of TF-DRC network based centrality measures.
|
Gene ID
|
Description
|
Betweenness
|
Degree
|
Family
|
|
|
Value
|
Rank
|
Value
|
Rank
|
|
105155794
|
transcription factor IIIA
|
27389.361
|
1
|
95
|
1
|
zf-C2H2
|
105157040
|
myb-related protein 306
|
6199.0103
|
2
|
46
|
2
|
MYB
|
105166097
|
dehydration-responsive element-binding protein 2D-like
|
5865.903
|
3
|
39
|
4
|
AP2
|
105161312
|
zinc finger protein ZAT5
|
5856.9683
|
4
|
25
|
13
|
zf-C2H2_6
|
105171783
|
myb-related protein Myb4
|
5674.037
|
5
|
44
|
3
|
MYB
|
105168111
|
transcription factor MYB63
|
4036.9114
|
7
|
38
|
5
|
MYB
|
105155728
|
myb-related protein 306
|
3645.386
|
8
|
33
|
6
|
MYB
|
105162220
|
myb-related protein Myb4-like
|
3306.7559
|
9
|
30
|
10
|
MYB
|
105165703
|
transcription factor PIF1
|
2107.3755
|
19
|
30
|
9
|
HLH
|
105166028
|
transcription factor E2FC
|
2938.6282
|
11
|
9
|
95
|
E2F_TDP
|
MCODE algorithm was applied to the integrated miR-TF-DRG network of each set to detect key regulatory relations of Sesame under drought stress (Additional file 1: Table S13-14). 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 DTR network (Fig. 6). Module c1 with 25 nodes and 38 edges was the largest detected cluster in the DRC network. Notably, genes such as DHN1 (105165756), ABP19a (105173675) (module a1), pectate lyase 18 (105157315) (module a2), DHN1-like (105179059) (module a3), pectate lyase 8 (105171845) and sucrose synthase 2 (105163347) (module a7), TFs like MYB63 (105168111), DREB 2D-like (105166097) (module c1), MYB-related 306 (105157040), PIF1 (105165703) (module c2), bZIP16 (105159270), CPRF2 (105172269) (module c3) and transcription factor IIIA (105155794) (module c4), and miRNAs like miR169a, miR399g (module c1), miR156b (module c2) were among detected modules in DRC network. Module d1 was also the largest cluster in the DTR network with 27 nodes and 41 edges. Genes like SUC2-like (105176204), CYP82D47-like (105178083) (module d1), cytokinin dehydrogenase 7 (105180056) (module d2) and TFs like MYB108 (105178858), WRKY30 (105160383), NAC29 (105175179) (module d1), and miRNAs including miR171b, miR529d (module d1), miR169a, miR395c (module d2), were found in modules of DTR network.
Table 3 Top 10 hub TFs of TF-DTR network based centrality measures.
|
Gene ID
|
Description
|
Betweenness
|
Degree
|
Family
|
|
|
Value
|
Rank
|
Value
|
Rank
|
|
105165703
|
transcription factor PIF1
|
2581.837
|
1
|
11
|
7
|
HLH
|
105168958
|
homeobox protein knotted-1-like 1
|
2401.3303
|
2
|
5
|
36
|
KNOX1,2
|
105172269
|
light-inducible protein CPRF2
|
1746.6625
|
5
|
11
|
8
|
bZIP_C
|
105177566
|
basic leucine zipper 43
|
1578.3323
|
6
|
8
|
13
|
bZIP_1
|
105178858
|
transcription factor MYB108
|
1504.9016
|
7
|
5
|
39
|
MYB
|
105178907
|
calmodulin-binding transcription activator 5
|
1357.1499
|
8
|
8
|
14
|
CG-1
|
105174207
|
dof zinc finger protein DOF3.5
|
1300.5356
|
9
|
8
|
12
|
zf-Dof
|
105169008
|
transcription factor PIF3
|
984.0747
|
16
|
9
|
10
|
HLH
|
105179013
|
probable WRKY transcription factor 50
|
1070.4827
|
14
|
6
|
31
|
WRKY
|
105178586
|
protein FAR-RED ELONGATED HYPOCOTYL 3
|
989.2037
|
15
|
8
|
15
|
FAR1
|
Table 4 Top 10 hub miRNAs of miRNA networks based centrality measures.
|
miR-DRC
|
|
miR-DTR
|
Gene
|
Betweenness
|
Degree
|
|
Gene
|
Betweenness
|
Degree
|
|
Value
|
Rank
|
Value
|
Rank
|
|
|
Value
|
Rank
|
Value
|
Rank
|
miR9773
|
21457.668
|
1
|
41
|
1
|
|
miR9773
|
13463.416
|
1
|
28
|
1
|
miR156f
|
14251.105
|
2
|
32
|
2
|
|
miR5658
|
4960.732
|
3
|
13
|
6
|
miR156e
|
13407.189
|
3
|
31
|
3
|
|
miR156f
|
4769.096
|
4
|
13
|
5
|
miR5658
|
10172.956
|
4
|
24
|
5
|
|
miR156b
|
4678.1597
|
5
|
13
|
4
|
miR156b
|
8392.074
|
5
|
27
|
4
|
|
miR169i
|
3451.3743
|
8
|
7
|
43
|
miR529d
|
6910.4907
|
7
|
18
|
9
|
|
miR156a
|
2916.2463
|
17
|
12
|
8
|
miR396a
|
4287.992
|
20
|
19
|
7
|
|
miR169a
|
3415.93
|
9
|
9
|
17
|
miR172d
|
6028.376
|
10
|
19
|
8
|
|
miR395k
|
3137.5925
|
12
|
9
|
16
|
miR396b
|
3994.1416
|
24
|
17
|
10
|
|
miR529d
|
2701.2659
|
18
|
8
|
33
|
miR156a
|
5670.0586
|
11
|
22
|
6
|
|
miR828b
|
2529.5598
|
21
|
9
|
19
|