The sequenced mungbean genotype VC1973A (named Zhonglu in China) was used for all experiments (Kang et al., 2014). Mungbean seeds were germinated and grown in pots in a growth room, with 16 h 28°C light/8 h 22°C dark cycles. Humidity was maintained at approximately 30%.

Two weeks old mungbean seedlings were used for stress treatments. Cold stress treatment was performed by transferring plants to a 4°C chamber for 6 h. Heat stress was carried out by putting plants in a 40°C chamber for 6 h. For the drought stress treatment, plants were grown without watering for 6 days. For the salt stress treatment, plants were watered with 200 mmol salt solution, and tissues were collected 48 h after treatment. Shoots and roots were collected separately after stress treatments and then stored at -80°C before RNA extraction. Moreover, various tissues, such as roots, stems, leaves, flowers, pods and seeds, were also sampled for gene expression analysis. Roots, stems and leaves were collected from 3 weeks old seedlings, and flowers were sampled from 6 weeks old plants. The pods were sampled for analysis at the beginning of the pod stage. The seeds, at the full seed stage were used for gene expression analysis. Each sample was analyzed using three biological replicates.

The conserved amino acid sequence of DNA-binding domains (DBD, PF00447) was download from the protein family database Pfam¹, and full-length amino acid sequences of Hsf proteins from Arabidopsis, soybean and potato were used as BLAST queries against the mungbean database² and National Center for Biotechnology Information (NCBI). All output genes with default (Limit Expect Value 1e-5) were analyzed using the Pfam database (E = 1.0) and SMART³ to remove genes without conserved domain sequences. The classification of VrHsf proteins was performed using Heatster⁴. The theoretical iso-electric points, grand average of hydropathicity and protein molecular weight analyses were performed using the ProtParam tool⁵.

In total, 124 Hsf amino acid sequences from mungbean, soybean, Arabidopsis and potato were used for phylogenetic analysis (Nover et al., 2001; Chung et al., 2013; Tang et al., 2016). The N-proximal regions of Hsf proteins, from the start of the conserved DBD domain to the end of the HR-A/B region, were aligned using Clustal-X2. The alignment result was used to construct a phylogenetic tree using MEGA 6.0’s Neighbor-Joining method with pairwise deletion, 1000 bootstraps and a Poisson model (Tamura et al., 2013).

The CDS and genomic DNA sequences of mungbean Hsf genes were aligned with the Gene Structure Display Server program (GSDS)⁶ to illustrate exon/intron organization (Hu et al., 2015). VrHsf conserved domains were obtained using Heatster⁷, Pfam and SMART. The conserved motifs of VrHsf proteins were assessed via MEME tools⁸, and the parameters of the maximum number of motifs and the optimum motif widths were 16 and 6–50 amino acid residues, respectively. Promoter cis-elements were analyzed using the Plant cis-acting regulatory DNA elements (PLACE) database (Higo et al., 1999) and the distribution maps were constructed using iBS⁹.

The MapInspect software¹⁰ was used for mapping VrHsf genes to mungbean chromosomes. Duplications of Hsf genes in mungbean were analyzed and marked as previously described (Zhang et al., 2015; Tang et al., 2016). The divergence time (T) of the duplicated genes was calculated as described (Lynch and Conery, 2000; Chen et al., 2014).

All mungbean RNA was extracted using Qiagen RNeasy mini kit following the instructions (Qiagen, United States)¹¹. cDNA synthesis was conducted as previously described (Li et al., 2017). SuperScript II reverse transcriptase first-strand synthesis kit (Invitrogen) was used for the synthesis of the first strand cDNAs with 2 μg total RNA according to the instructions. The LightCycler 480 SYBR Green I Master Kit (Roche Diagnostics) was used for quantitative real-time PCR (qRT-PCR) using a LightCycler480 machine (Roche Diagnostics), according to the manufacturer’s instructions. The amplification program for qRT-PCR was performed as previously described (Li et al., 2017). For gene expression analysis, three biological replicates were used for each sample and gene expression was normalized to an Actin-expressing gene in mungbean (Vradi03g00210). All primers used for qRT-PCR analyses are listed in Supplementary Table 1.

Statistical significance was analyzed by t-test using the SAS program (SAS Institute Inc.)¹².

We used the conserved DBD amino acid sequence (PF00447) and the full-length amino acid sequences of Hsf proteins from Arabidopsis, soybean and potato as BLAST queries against the mungbean genome database and NCBI GenBank resources (Kang et al., 2014). We ultimately identified 24 VrHsf candidate genes in the mungbean genome with complete DBD and HR-A/B regions (Supplementary Figures 1, 2) and analyzed their genomic length, CDS length, number of amino acids, theoretical molecular weight, grand average of hydropathicity, and isoelectric point (Table 1). The genomic lengths of the VrHsf genes ranged from 1,435 bp (Vradi06g15090) to 12,757 bp (Vradi08g09150), the coding sequence sizes varied from 561 bp (Vradi08g02500) to 1,521 bp (Vradi07g03150), and the deduced number of amino acids ranged from 186 to 506. The grand averages of hydropathicity fell between -0.484 and -0.891, indicating that all the VrHsfs were predicted to be hydrophilic proteins. In addition, the predicted molecular weights ranged from 22.16 to 56.26 kDa, and the isoelectric points ranged from 4.69 to 9.30, indicating the structural diversity and functional variation of VrHsfs (Table 1).

The VrHsf proteins are classified into three groups (A, B, and C) based on differences in the numbers of amino acids inserted between the HR-A and HR-B cores (Guo et al., 2016). Class A has a 21 amino acid insertion between the HR-A and HR-B regions, and Class C has a 7 amino acid insertion, whereas Class B Hsfs are more compact (Scharf et al., 2012). We aligned HR-A/B regions and used Heatster to classify VrHsfs into three groups. Among the 24 VrHsf members, 13 VrHsf proteins were grouped into class A, 10 VrHsf members belonged to class B and only 1 member was classified as class C (Table 1 and Supplementary Figure 2). We constructed a phylogenetic tree based on 124 Hsf amino acid sequences, including 24 mungbean Hsfs, 52 soybean Hsfs (Scharf et al., 2012), 21 Arabidopsis Hsfs (Scharf et al., 2012), and 27 potato Hsfs (Tang et al., 2016), using the conserved regions from the start of the conserved DBD domain to the end of the HR-A/B domain (Figure 1 and Supplementary Table 2) to investigate the evolutionary relationships among Hsf families and to gain insight into the potential function of VrHsfs from the well-studied Hsf families in other species. VrHsfs in class A were grouped into 6 distinct sub-classes (A1, A3–A7), and no VrHsf members were classified as A2, A8, or A9. VrHsfs in class B were sub-classified into five groups (B1-B5) (Table 1, Figure 1, and Supplementary Table 2). We also created a phylogenetic tree using only VrHsf N-proximal regions (Figure 2A). The class A members formed a single group in comparison with class B, which was divided into 3 sub-groups. VrHsfB1 sub-subclass had a distinct relationship with other VrHsfB sub-families, and VrHsfB5 had a much closer relationship with class A and class C families (Figure 2A).

We constructed the exon-intron organization of the 24 VrHsf genes using the genomic and coding sequences. The gene structures showed that all the VrHsf genes exhibited a highly conserved exon-intron organization, with two exons and one intron (Figure 2B). The intron lengths varied, similar to Arabidopsis (Nover et al., 2001). We used MEME to predict putative motifs to further reveal the conservation and diversity of VrHsf genes (Figures 2C, 3 and Table 2), and we identified 16 distinct motifs among mungbean Hsf proteins (Figure 3). All the VrHsf members contain motifs 1, 2, 3, and 4. Motif 1 is closely connected with motifs 2 and 3. Combined, motifs 1, 2, and 3 represent the most conservative domain, the DBD domain. We considered motifs 4 and 5 the HR-A/B regions (Figure 2C). Some motifs were only found in specific VrHsf proteins, for example, motif 5 is specific to the class A family, and most class B members have motif 8, except for VrHsfB5 and VrHsfB3a (Figure 2C). Moreover, we identified NLS motifs in 21 VrHsf proteins and NES motifs in 9 class A and 2 class B VrHsfs. The class A Hsfs have short activator peptide motifs (AHA) and class B Hsfs contains the repressor domain tetrapeptide (Scharf et al., 2012). We detected AHA motifs, which may work as transcriptional activators, in 10 class A VrHsfs, while VrHsfA1c, VrHsfA1d, and VrHsfA3b have no AHA domains (Table 2). Eight of the 10 class B VrHsf genes contain the tetrapeptide LFGV, the RD domain (Table 2), and therefore may function as repressors in stress response pathways (Scharf et al., 2012). In addition, VrHsfC1 has no RD or AHA domains. In sum, these observations indicate the functional divergence among VrHsf genes.

The VrHsf gene family members have some common functions, and some members may have evolved new functions related to the original genes after gene duplication. The chromosomal locations of VrHsf genes can represent gene distributions after duplication. To map the locations of VrHsf genes on different chromosomes, we determined the distribution of VrHsf genes based on the mungbean genome database. Two candidate genes, VrHsfA5b and VrHsfA5c, were discarded due to the lack of chromosome information (Table 1). VrHsf genes were located on 8 of the 11 mungbean chromosomes, with no VrHsfs found on chromosomes 2, 4, or 9 (Figure 4). Chromosome 8 contained the most VrHsf genes, with four class A and two class B genes (VrHsfA3a, VrHsfA6a, VrHsfA6b, VrHsfA7b, VrHsfB2c, and VrHsfB3a), followed by chromosomes 3, 7, and 11, with three genes each (Figure 4). To investigate the duplication events, which may have occurred during mungbean genome evolution, we used the MpInspect software to analyze mungbean Hsf genes. The analysis identified six interchromosomal duplications and one intrachromosomal duplication (Figure 4). Class A proteins contained four duplication events, class B contained three duplication events and class C had no duplicated genes, and these gene duplications indicated similar functions for the duplicated gene pairs. In addition, chromosome 8 contained the most duplicated gene pairs, with five duplication events (VrHsfB2a/VrHsfB2c, VrHsfB3a/VrHsfB3b, VrHsfA3a/VrHsfA3b, VrHsfA6a/VrHsfA6b, and VrHsfA7a/VrHsfA7b) and all VrHsf members located on chromosome 8 had duplicated genes, indicating that chromosome 8 contained the original genes of many duplicated VrHsf genes. In contrast, genes on chromosome 10 had no duplication events (Figure 4).

We next investigated putative cis-elements in the promoter region 2 kb upstream of the translation initiation codons of VrHsf members. We obtained many cis-elements in the promoter regions including five known abiotic stress response elements (Figure 5). The stress related elements, Heat Stress Element (HSE), Low Temperature Responsive Element (LTRE), Dehydration-Responsive Element (DRE), C-Repeat Binding Factor (CBF), and ABA Responsive Element (ABRE) were characterized in VrHsf promoter regions (Figure 5). Distribution analysis of these cis-elements showed that all the VrHsf promoters contain HSE, CBF, ABRE, and DRE elements, and 17 of the 24 VrHsf promoter regions contain LTRE elements. All genes had multiple DREs and ABREs, pointing to their key roles in response to drought stress and ABA response pathways (Figure 5). These observations of the cis-elements in promoter regions imply that the VrHsf genes exhibit functional diversity and might be responsive to many different abiotic stresses.

To address the potential functions of VrHsf genes in different tissues, we used qRT-PCR to analyze their transcription patterns in various tissues, including the root, stem, leaf, flower, pod and seed (Figure 6). The expression patterns of each VrHsf gene varied in different tissues, indicating their potential functions in these tissues. Among the class A family, VrHsfA1c, VrHsfA1d, and VrHsfA5a were expressed at relatively high levels in all tissues, while VrHsfA4 was expressed at low levels (Figure 6). Moreover, VrHsfA1c, VrHsfA1d, and VrHsfA5a were more highly expressed in the root than in other tissues, indicating their critical roles in roots. Among the VrHsfB family, VrHsfB2a, VrHsfB2b, and VrHsfB4c were also expressed at high abundances in all tissues. In contrast, VrHsfB3a, VrHsfB3b, and VrHsfB5 were expressed at low levels, and were not expressed in the pod, seed, and leaf (Figure 6). In addition, the VrHsfC1 gene showed constitutively low expression in all tissues (Figure 6).

Many duplicated gene pairs displayed similar expression patterns in some tissues (Figure 6). For example, VrHsfA1c and VrHsfA1d showed similar expression levels in the root, leaf and seed. VrHsfA6a and VrHsfA6b expression patterns closely resembled each other in the flower and stem. VrHsfB4a and VrHsfB4c displayed similar expression patterns in the stem. Moreover, expression of VrHsfB3a conformed to VrHsfB3b in the pod, seed and leaf, indicating closely related functions in these tissues for these duplicated genes. In contrast, duplicated genes showed distinct expression levels in some tissues (Figure 6), indicating that they may have evolved new functions compared to the original genes.

To investigate the potential functions of VrHsf genes in response to different stresses, we analyzed VrHsf gene expression in mungbean shoots and roots under cold, drought, heat and salt conditions. The expression of most VrHsf genes changed in roots or shoots under stress treatments (Figure 7). The expression patterns of each VrHsf gene varied under different stress treatments. Expression of VrHsfA1d, VrHsfA7b, VrHsfB2b, and VrHsfB4a increased sharply under cold treatment in the shoot. Under drought stress, VrHsfA6a and VrHsfA6b exhibited a more than 10-fold increase in expression level in the root, compared with plants grown under normal conditions. The expression of many genes increased more than 10-fold under heat stress, including VrHsfA4, VrHsfA7a, VrHsfA7b, VrHsfB2c, VrHsfB4c, VrHsfB4h, VrHsfB5, and VrHsfC1 (Figure 7). Under salt treatment, VrHsfA6b showed the greatest change in expression with more than an 11-fold increase. The expression alterations of these VrHsf genes suggest their potential roles in response to the related stresses. In addition, the duplicated gene pairs displayed similar expression patterns under certain stresses (Figure 7), such as VrHsfA3a/VrHsfA3b, VrHsfA7a/VrHsfA7b, and VrHsfB2a/VrHsfB2c in both root and shoot under heat stress, VrHsfA6a/VrHsfA6b in root under drought stress, suggesting the functional redundancy of these duplicated genes.