The Hsf genes sequences of Arabidopsis and Utricularia gibba were retrieved from the Plant Transcription Factor DataBase (http://planttfdb.cbi.pku.edu.cn/) (Jin et al., 2014). Based on the genome and proteome sequences of sesame downloaded from the Sinbase (http://ocri-genomics.org/Sinbase/) (Wang et al., 2014), an extensive search was performed to identify all members of the Hsf gene family. Using HMMER3 tool of the software Unipro UGENE (Okonechnikov et al., 2012), the Hsf domain, denominated as PF00447 in the Pfam HMM library (http://Pfam.sanger.ac.uk/) (Finn et al., 2014), was searched against the sesame proteome.

The SMART program (Letunic et al., 2015) and Pfam were used to confirm the signature DBD domain of Hsfs (Bateman et al., 2004). All candidate Hsf genes with no DBD were removed. Additionally, the MARCOIL program (Delorenzi and Speed, 2002) was used to detect the coiled-coil structures and eliminate any sequence without a coiled-coil structure. The physical and chemical parameters of the retained Hsf genes were analyzed using ProtParam tool (http://web.expasy.org/protparam/). The classification of candidate Hsf proteins into three different classes A, B, and C was performed by Heatster (http://www.cibiv.at/services/hsf/) (Scharf et al., 2012) based on the information of oligomerization domains (Nover et al., 2001).

To understand the evolutionary relationships between Hsf genes in sesame, Arabidopsis and U. gibba, multiple alignments of the N-proximal regions of Hsf proteins from the three species were conducted using ClustalW program built in the MEGA 6.0 software (Tamura et al., 2013) with a gap open penalty of 10 and gap extension penalty of 0.2. The result of alignment was used for the construction of an un-rooted Maximum-Likelihood (ML) tree with a 1000 bootstrap value.

Based on Sinbase information, exon-intron substructure map was constructed by the Gene Structure Display Server (GSDS 2.0), web-based bioinformatics tool (http://gsds.cbi.pku.edu.cn/) (Hu B. et al., 2015). Conserved motifs of Hsf proteins were analyzed using MEME Suite version 4.10.2 (Bailey et al., 2009) with the parameters set as follows: minimum width of 6, maximum width of 50, the maximum number of motifs was set to 20 and the number of repetitions set to “any number”. All other parameters were set at default.

The Hsf genes were mapped onto the 16 Linkage Groups (LGs) of sesame genome on the basis of the information from Sinbase using the software MapChart 2.3 (Voorrips, 2002). The duplicated genes were identified and connected by lines according to Lin et al. (2011). As regard to estimating the synonymous (Ks) and non-synonymous (Ka) substitution rates, the cDNA and the protein sequences of the duplicated genes pairs were submitted to PAL2NAL (http://www.bork.embl.de/pal2nal/) (Suyama et al., 2006). The approximate date of the duplication event was calculated as described by Wang et al. (2015).

In addition, the amino acid sequences of predicted Hsf proteins were BLASTp searched against protein sequences of Arabidopsis in NCBI. Hits with E-value ≥ 1e-40 and at least 70% similarity were considered significant. The relationships between the orthologous genes of sesame and Arabidopsis were plotted using Circos (http://circos.ca/) (Krzywinski et al., 2009).

The presence of simple sequence repeats (SSRs) were searched against the DNA sequences of the Hsf genes using the web based software Websat (http://wsmartins.net/websat/) (Martins et al., 2009) with the following parameters: two to six nucleotide motifs were considered and the minimum repeat unit was defined as five reiterations for di-nucleotides and four reiterations for other repeat units (Dossa et al., 2016b). Furthermore, specific protein interaction networks based on Arabidopsis interactome were constructed applying STRING software (http://string-db.org/) (Szklarczyk et al., 2011).

To analyze the expression patterns of sesame Hsf genes, we used RNA-seq data from root, stem, seeds at different stages and leaf downloaded from SesameFG (http://202.127.18.220/hg/). Reads per kilobase per million mapped reads (RPKM) values were log10 transformed and, finally, the heat map of hierarchical clustering was constructed using the MEV Software (Saeed et al., 2006).

A total of 30 Hsf genes were confirmed in sesame genome with apparently complete Hsf-type DNA-binding domains and oligomerization domains ranging from 255 to 2646 bp in length (Table 1, Supplementary Table 2). The deduced proteins sizes vary from 84 amino acids (HSF9) to 889 amino acids (HSF21) and the isoelectric point ranging from 4.73 (HSF6) to 9.45 (HSF25). Table 2 provides more detailed information about each protein physicochemical characteristics.

Except for 3 genes (HSF28, HSF29, and HSF30) which were located on the unanchored scaffolds that could not be mapped to a specific LG, all the Hsf genes were unevenly distributed among 12 LGs out of the 16 LGs of sesame genome, suggesting that Hsf genes may have been widely distributed in the genome of the common ancestor (Figure 1). The numbers of sesame Hsf genes in each LG vary widely with LG3 and LG6, harboring the highest number of Hsf genes (3; 10%) while LG5, LG10, LG14, and LG15 have only one gene (3.33%).

Analyzing of the duplication events can help to elucidate the evolution mechanisms of the sesame Hsf gene family. In total, 7 duplicated gene pairs across 11 LGs were identified including six segmental duplication events between different LGs (HSF1 and HSF4, HSF12 and HSF13, HSF11 and HSF18, HSF10 and HSF21, HSF7 and HSF20, HSF17 and HSF26) and one duplication event within the same LG6 (HSF14 and HSF15). The duplicated gene pairs are linked with black lines in Figure 1.

To get more insights into the evolution of these duplicated genes, synonymous (Ks), non-synonymous (Ka) and their ratio (Ka/Ks) values were used to examine the selective pressure and duplication time of the 7 sesame Hsf segmentally duplicated genes. In general, a Ka/Ks > 1 means positive selection, Ka/Ks < 1 indicates purifying selection, and Ka/Ks = 1 stands for neutral selection (Nekrutenko et al., 2002). The Ka/Ks ratio for all sesame Hsf duplicated genes ranged from 0.1367 to 0.263, thus Ka/Ks < 1 for all the 7 duplicated gene pairs (Supplementary Table 3). These results suggested that the duplicated Hsf genes were under strong purifying selection pressure. In addition, the duplicated events were estimated to have occurred approximately ~67 million year ago (51.5–93.7 MYA; Figure 2).

All the 30 Hsf genes were submitted to Heaster program to uncover the conserved functional domains. The details of various functional domains such as DBD, HR-A/B, NLS, NES, and AHA motifs are presented in Table 3. As the most conserved domain in the Hsf genes, DBD composed of approximately 90 amino acids was found in all the 30 genes. In addition to DBD, HR-A/B, another core conserved domain, was also present in all sesame Hsf proteins, while the other conserved domains were specific to each class. According to the distinction between the HR-A and HR-B motifs, Hsfs were divided into A, B, and C classes. The class A was the most representative one with 16 genes out of 30 genes. Twelve genes belonged to the class B and the remaining 2 sesame Hsf genes were classified into the class C.

MEME web server was used to verify the results of domain prediction (Figure 3A; Supplementary Table 4). Twenty different motifs were found distributed throughout the Hsf proteins sequences, lengths ranging from 8 to 50 amino acids. Most of sesame Hsf genes displayed motifs 1, 2, and 5 which corresponded to the DBD domain whereas some motifs were characteristic of each class. Moreover, each class exhibited different motif organizations and two closely clustered genes on the tree generally displayed similar motif patterns. Overall, despite variability in size and sequence, the predicted Hsf DBD, HR-A/B region and other conserved domains were cross-confirmed in each sesame Hsf gene by the two combined methods.

Analysis of the exon/intron boundaries of sesame Hsf genes coupled with a phylogenetic tree exhibited a highly conserved exon/intron organization in 24 genes out of the 30 Hsf genes possessing strictly two exons. In addition, we identified 2 intronless genes (HSF9 and HSF5), 3 genes possessing 3 exons (HSF3, HSF7, and HSF23) and the intriguing gene HSF24 displaying up to 12 exons (Figure 3B). HSF24 might contain the original genes (Nuruzzaman et al., 2010). In general, although most closely clustered genes in the phylogenetic tree displayed similar numbers of exon, there is a wide variation among the organizations and lengths of exon/intron.

To evaluate the evolutionary relationships of sesame Hsf genes, phylogenetic analysis was performed based on the full-length amino acid sequences of Hsf proteins from sesame, U. gibba and Arabidopsis. The tree clearly distinguished the three classes, namely A, B, and C (Figure 4). Based on the classification of the Hsf genes of the well-studied crop Arabidopsis, the different functional groups belonging to each class have been identified. The class A members were further divided into 8 subclasses (A1–A8) with the subclass A1 more represented (4 members). The class B members were also classified into 4 subclasses namely B1, B2, B3, and B4, while all of the class C members were found to belong to the same subclass C1.

Synteny analysis of sesame Hsf genes compared to those of Arabidopsis could provide more functional insight. Hence, we further performed a genome-wide comparative analysis to identify orthologous Hsf genes between sesame and Arabidopsis genomes. In total, 18 orthologous gene pairs were found across 9 LGs in sesame (Figure 5). The LG3, LG6, and LG7 exhibited the largest orthology (3 gene pairs) while the least orthologous genes (1 gene pair) were detected on LG2, LG4, LG5, LG10, and LG15. The orthologous gene pairs and their localization in each genome are presented in Supplementary Table 5. Out of the three classes of sesame Hsf genes, only two classes (A and B) were represented in the orthologous gene pairs. Moreover, 5 Hsf genes in Arabidopsis genome retained 2 genes in the sesame genome while conversely, only one gene in sesame genome (HSF19) conserved a tripled copy in Arabidopsis genome (AT5G16820, AT3G02990, and AT1G32330).

The protein–protein interaction network is a helpful preface to explore the biological functions of unknown proteins. Based on the interactome of Arabidopsis, the protein–protein interactions, including functional and physical interactions among Hsf genes of sesame were predicated (Figure 6). Thirteen sesame Hsf genes were involved in different interaction patterns. The Arabidopsis genes found in this network are mostly described to be related to defense, plant development, proteins folding and abiotic stress related genes. Hence, we hypothesized that the sesame Hsf genes might also be involved in similar function pathways. The first group including 11 sesame Hsfs (HSF1, HSF10, HSF12, HSF22, HSF8, HSF20, HSF5, HSF27, HSF16, HSF3, and HSF18) interact directly with HSP genes, hence they might be some positive or negative regulators of HSP genes in responses to abiotic stresses and plant defenses. The 2 genes HSF19 and HSF9 seemed to be involved in the transcription regulation of other Hsf genes.

Based on the importance of Hsf genes in responses to various environmental stresses, tagging these useful genes could help in improving marker-aided breeding and establishing better crop-improvement strategies. Out of the 30 genes, only 7 SSR markers covering 4 LGs were identified (Supplementary Table 6). Two SSR types were observed: di-nucleotide and the most abundant one, tri-nucleotide motifs.

The expression patterns of sesame Hsf genes in different organs and tissues viz. root, leaf, stem tip, seed with different coat colors harvested at various development stages, were investigated. Heat maps using both phylogenetic tree and hierarchical clustering methods were generated based on the RPKM values for each gene in all samples (Figure 7, Supplementary Figure 1). The results showed that almost all sesame Hsf genes were expressed in the different tissues and organs. The hierarchical cluster analysis divided the 30 genes into five groups which were a mix of genes from different Hsf classes. Interestingly, expression levels differed from organ to organ and from Hsf gene to another even in the same class. For instance, HSF4 were expressed highly in all tissues and organs while HSF1 belonging to the same class as HSF4 displayed a strong expression level in root, leaves and stem but low expression levels in seeds. The result suggests that just because the phylogenetic relationships of these genes are close doesn't mean they may develop similar biological functions. It is worth noting that most of genes displayed higher expression in root compared to the other tissues. This is understandable since the root is a fast-growing tissue whose metabolism and development processes are all more active than in other tissues. Moreover, Hsf genes displayed mostly their lowest expression levels in sesame seeds. Finally, the Hsf class B members appeared to be the most expressed genes in organs and tissues of sesame compared to the other two classes.

To understand the potential role of Hsf genes in sesame responses to drought, we monitored all sesame Hsf genes expression levels during progressive drought and after re-watering in two contrasted accessions (drought-tolerant and drought-sensitive). Since most Hsf genes displayed higher expression in root compared to the other tissues, qRT-PCR analysis was performed in sesame root samples. The results showed that some genes were highly expressed while others had low expression levels and revealed 3 types of expression patterns: (1) Up-regulated genes, (2) down regulated genes, and (3) no regulated genes (Figure 8). In the first category (1), these genes showed increase of their expression levels during drought stress. This increase of expression levels could be observed early at the 3rd day (HSF1, HSF2, HSF4, HSF11, HSF16, HSF22, HSF24, HSF25, HSF28, HSF29, and HSF30), or tardily at the 6th day (HSF7, HSF13, HSF20, and HSF26). For the second category of genes (2), they displayed the decrease of their expression levels under drought stress (HSF3, HSF10, HSF12, HSF14, HSF15, HSF17, HSF18, HSF19, HSF21, HSF23, and HSF27). Hence, in most cases, while the expression levels of the genes in category (1) decreased after re-watering period, the genes in category (2) increased their expression levels. Finally, the genes in category (3), showed minor fluctuations during drought stress and after re-watering (HSF5, HSF6, and HSF8).

Comparing expression patterns of Hsf genes in the drought-tolerant accession and the drought-sensitive one, results showed that most of time, the expression levels of Hsf genes were broadly higher in the drought tolerant accession. Overall, it appeared that regulation by drought is not specific to some Hsf classes but depends on the genes.

Benefiting from genome availability, Hsf gene family has been extensively studied in many crops. Twenty five genes were reported in Arabidopsis (Jin et al., 2014), 35 in Chinese cabbage (Song et al., 2014), 25 in rice (Guo et al., 2008), 26 in tomato (Yang et al., 2016), 22 in U. gibba (Jin et al., 2014), 27 in potato (Tang et al., 2016), 25 in maize (Lin et al., 2011), 29 in Chinese white pear (Qiao et al., 2015), 25 in Chinese pepper (Guo et al., 2015) etc., showing the multiplicity of this gene family in plants. In this study, we reported 30 Hsf genes in sesame genome, classified into 3 major classes. As reported in other crops, the class A genes are the most represented, followed by the class B members (Guo et al., 2008, 2015; Lin et al., 2011). However, no Hsf member belonging to the class A9 (At5g54070) was found in sesame similarly as reported in Chinese cabbage (Huang et al., 2015).

Compared to its immediate relatives U. gibba (~82 Mb), potato (~840 Mb), tomato (~950 Mb) and the model plant Arabidopsis (~135 Mb), sesame (~370 Mb) bears more Hsf genes in its genome, suggesting that the number of Hsf genes is not necessarily linked to genome size but might be explained by the different evolutionary histories occurred within each species. Whole genome duplications have had a strong impact on species diversification and may have triggered evolutionary novelties (Jaillon et al., 2009). WGDs within the angiosperms presumably are the cause of varying numbers of Hsfs between different plant species (Scharf et al., 2012). Sesame is estimated to have diverged from the tomato-potato lineage approximately 125 MYA (million years ago) and from U. gibba approximately 98 MYA (Wang et al., 2015). Tomato-potato lineage has recently experienced a whole genome triplication (WGT) and it was expected to have higher copies of Hsf genes in their genomes while sesame and U. gibba underwent independently a WGD at different points of evolution. Hence, we hypothesized that there might be some extensive Hsf genes loss after WGT events principally in the tomato-potato lineage. Conversely, known as arid areas crop, sesame might retain for adaptation purposes a large part of its Hsf genes during evolution. These observations are consistent with Lin et al. (2014) who demonstrated that the vast majority of Hsf gene duplications resulted from whole genome duplication events rather than tandem duplication and significant differences in gene retention exist from species to species.

To identify the evolutionary origins of Hsf genes in sesame, we further analyzed the distribution of the duplicated genes. Seven pairs of duplicated genes were found in sesame Hsf family, which was less than potato (8 pairs) but higher than tomato (4 pairs). Six gene pairs are involved in segmental duplication while the pair of genes HSF6 and HSF14 is involved in regional duplication within the LG6 and has probably evolved from their common ancient gene through the gene duplication event within this LG. These results suggest that segmental duplications were the primary force underlying the expansion of Hsf genes in sesame. Our finding corroborates well the previous reports in other plants about expansion of Hsf genes thought to result primary from segmental duplication (Lin et al., 2011, 2014; Song et al., 2014; Guo et al., 2015). Because most plants with diploidized polyploids retained numerous duplicated chromosomal blocks within their genomes, segmental duplication occurred more frequently than tandem duplication and transposition events in plants (Cannon et al., 2004). Based on synonymous substitution rates (Ks) which is a common procedure to determine the evolutionary age and divergence level of gene copies (Wang et al., 2013), these duplication events date back to around 67 MYA, similarly as predicted by Wang et al. (2015) who proposed a period between 54 and 72 MYA.

The close relationship between sesame and the well-studied plant Arabidopsis helped in identifying some orthologous genes and allowed us to predict their functions in sesame. Indeed, gene orthology analysis can be used as a preliminary method to explore the function of candidate genes (Wang M. et al., 2014). Only class A and class B genes were found to have orthologs Hsf in Arabidopsis genome and in-depth screening of these orthologs revealed two main groups of Hsfs: the most represented are those involved in regulation of transcription of others genes and few are involved in responses to heat and other abiotic stresses. Hanada et al. (2008) reported that transcription factors encoding nucleic acid binding proteins originated mostly through segmental duplication. Hence, we deduced that segmental duplication events in the expansion of Hsf genes in sesame may induce transcriptional regulation functions to these genes. For instance the gene HSF2, HSF8, HSF9, HSF16, and HSF27 are orthologs of Arabidopsis genes (AT2G26150, AT4G36990, and AT132330) reported to act as activator or repressor of heat inducible Hsf genes in Arabidopsis (Ikeda et al., 2011; Weng et al., 2014). Furthermore, based on the interalog of Arabidopsis, the protein-protein interactions of some sesame Hsf genes were constructed, and results showed that most of genes interact with Hsps. Hsf genes have been demonstrated to play key roles in the tolerance to various stresses by reacting with different genes especially Hsps (Swindell et al., 2007; Hu et al., 2009). Unfortunately, the regulation of the Hsps by the Hsf genes in sesame is still not documented and future works should focus on this issue to unravel the mechanisms of regulation as well as the possible interaction with other transcription factors such as DREB2A as demonstrated by Scharf et al. (2012).

Some authors described Hsf genes as tissue- and stage-specific in several plants (Swindell et al., 2007; Giorno et al., 2010). This was not the case in our study and, as previously reported in L. japonicus (Lin et al., 2014), all sesame Hsf genes were expressed in all organs and tissues investigated, suggesting the functional conservation of this gene family and their regulatory roles at multiple developmental stages. However, they displayed different expression levels, pointing out their functional differences. Data also showed that some duplicated genes exhibited different expression patterns in organs and tissues of sesame. This is, for instance, the case of the duplicated genes HSF1 and HSF4, indicating that as a major feature of the long-term evolution, duplicated genes have diversified their functions (Blanc and Wolfe, 2004). Class B-Hsfs and class A-Hsfs genes were the most abundant in the different tissues and organs in sesame, showing their active roles in vegetative growth as well as seed maturation. Similar expression patterns of these genes were reported in barley, potato and legumes (Lin et al., 2014; Reddy et al., 2014; Tang et al., 2016).

Sesame is primarily grown for its seeds which bear one of the highest oil content. Thus, genes involved in sesame seed maturation are of great importance. Surprisingly, we detected that sesame genome lacks HsfA9 group which has been described to play a specific role in seed development in many species like Arabidopsis, tomato and tobacco (Kotak et al., 2007; von Koskull-Döring et al., 2007). Such observation implies that sesame might hold other functional genes which perform this parallel function. For instance, HSF4, HSF9, HSF11, and HSF15 exhibited abundant transcripts in sesame seeds whichever the seed color or the development stage. Hence these genes are likely to play key roles in sesame seed development. In rice, OsHsfA7 has also been proposed to play the seed-specific function (Chauhan et al., 2011).

The major objective for agronomic research is to enhance crop productivity under various abiotic stresses (Puranik et al., 2012). Heat stress is often compounded by additional abiotic stresses such as drought and salt stress (Bita and Gerats, 2013). Although it is well known that Hsf genes are involved in the acclimation of plant to heat, other stresses such as drought are less documented.

In this study, two contrasting accessions of sesame based on their response to drought were used for time-course Hsf genes expression and 90% of these genes were drought responsive. This result indicated the involvement of Hsfs in drought stress tolerance in sesame, which is in agreement with previous studies in different plant species (Chauhan et al., 2011; Bechtold et al., 2013; Li et al., 2014). In addition, we observed that most of these drought responsive Hsf genes in this study are highly expressed in the drought-tolerant accession compared to the sensitive one, confirming their key roles in drought tolerance. About 50% of Hsf genes were up-regulated while 36% were down-regulated, showing that these genes might play specific roles and sesame adapts to drought stress through compensation and interaction of its Hsf genes expression. Most of down-regulated genes were class A-Hsfs, while many up-regulated genes belonged to the class B and C. This result is intriguing, knowing the importance of some class A-Hsfs such as HsfA2 and HsfA1 members which have been identified to be responsive to heat stress in Arabidopsis and tomato (Mishra et al., 2002; Ogawa et al., 2007) and may probably be involved in other abiotic stresses such as drought (Ogawa et al., 2007; Bechtold et al., 2013; Zhang et al., 2015). However, under drought stress, most of these genes were down-regulated in sesame. In the works of Guo et al. (2015), they also highlighted the involvement of some class B-Hsfs in osmotic stress tolerance in pepper. More notably, they observed that class A-Hsfs expressed differentially under various stresses. In soybean, the Hsf type-A1 gene GmHsf-17, was strongly induced under heat stress but not influenced under drought (Li et al., 2014). These results indicate that there are species-specific and stress-specific features in the functions of Hsf members in regulating genes involved in plant stress responses and we propose class B-Hsfs as major transcriptional repressors or co-activators cooperating with some class A-Hsfs (Bharti et al., 2004; Scharf et al., 2012), to confer drought resistance in sesame.

Although HSF genes are well characterized in the heat stress-related pathway, they are poorly understood in other stress responses in Arabidopsis (Huang et al., 2016). For instance, the few studies focused on drought stress revealed the role of HsfA6b, HsfA6a, and HsfA1b in drought stress responses in Arabidopsis (Bechtold et al., 2013; Hwang et al., 2014; Huang et al., 2016). Through syntheny analysis between sesame and Arabidopsis, we found the gene HSF10 as orthologs of HsfA6b in Arabidopsis, HSF19 as orthologs of HsfA1b in Arabidopsis while no orthology was found for the gene HSFA6a in sesame. Gene expression analysis of these 2 genes (HSF10 and HSF19) in sesame revealed that they were significantly down regulated under progressive drought stress in the resistant accession as well as in the sensitive one. This suggests that these genes might also be involved in drought response pathways in sesame.

To cope with prolonged drought stress, sesame uses different set of Hsf genes belonging to class B4, class B3 and class A1 to maintain housekeeping gene expression. In Arabidopsis, similar pattern has been reported under heat stress with AtHsfA1a and AtHsfA1b involved in the early response to heat stress (Lohmann et al., 2004) while AtHsfA2 enhanced and maintained the response when plants are subjected to long-term stress (Charng et al., 2007). The two members of Class C-Hsfs were induced under drought stress in sesame and may play important function during drought stress. In contrast to class A-Hsfs and class B-Hsfs, members of class C have not so far been fully characterized. Further functional works are needed to excavate the role of these genes in drought tolerance in sesame. The unaltered genes (all belonging to the class A-Hsfs) might operate on the regulation pathways of other abiotic stresses (Victor and Benecke, 1998). During plant recovery from drought damage, the up-regulated genes expression levels slowly decreased while the down-regulated genes increased. These results imply that the functions of down-regulated genes are more potent under normal conditions in sesame plants. However, gene expression is a complex biological process and more thorough studies are needed to decipher the regulatory mechanisms between Hsfs and their co-expressed genes not only under drought stress but in combination with other stresses.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.