The gene sequence of B. napus were downloaded from the genome database¹, to gather the probable candidate B. napus Hsf amino acid sequences, the Hsf-type DBD domain (Pfam: PF00447) was submitted as a query in a protein BLAST search of the B. napus genome database. Hsf gene sequences from Arabidopsis (Nover et al., 2001) were retrieved from the TAIR database (Lamesch et al., 2012) and used as queries to perform repetitive BLAST searches against the Phytozome database v9.1 (Goodstein et al., 2012). BLAST searches were also performed against the NCBI nucleic-acid sequence data repositories. All protein sequences derived from the BLAST searches were examined using domain-analysis programs. Molecular weight, iso-electric point, functional domains, and amino acid signal peptides of BnaHsfs were calculated using the ExPASy online servers².

Multiple sequence alignment of Hsf proteins from B. napus were performed using the program ClustalX 1.83 (Thompson et al., 1997). The phylogenetic tree was constructed using the neighbor-joining (NJ) method by MEGA 6 program (Tamura et al., 2013), the bootstrap value was set at 1000 replications to assess tree reliability.

The MEME program³ (Bailey et al., 2009) was used for identification of conserved motifs, with the following parameters: number of repetitions: any; maximum number of motifs:15; and the optimum motif widths: 6–200 amino acid residues. Exton/intron organization of the Hsf genes in B. napus was illustrated using Gene Structure Display Server program (GSDS⁴) (Hu et al., 2015) by alignment of the cDNAs with their corresponding genomic DNA sequences.

Prediction of putative cis-elements was performed using Signal Scan Search against the plant cis-acting regulatory DNA elements (PLACE) database (Higo et al., 1999). A 2-kb sequence upstream of ATG initiation codon of BnaHsf genes was selected as the promoter region for this analysis.

Rapeseed seeds were germinated on a filter paper, and then transplanted to soil pots growing in the greenhouse at Oil crops research institute (Wuhan, China) with conditions of a temperature of 22°C, LED sodium lamp and a humidity of about 50–70% for 6 weeks. The plants were then transferred to growth chamber programmed under specific environmental conditions for 2–3 days before stress treatment. The conditions in growth chamber were set as follows: temperature of 25°C and humidity of 50% in 16 h light; temperature of 22°C and humidity of 60% in 8 h dark.

The high CO₂ stress was performed in a growth chamber (AR-41L2, United States) in which CO₂ gas can be accurately and stably controlled in the range of 100–3000 ppm. The conditions of growth chamber were set as follows: CO₂ concentration of 1000 ppm, light intensity of 100 umol/m²/s, temperature of 25°C and 60% relative humidity. Leaf samples were collected at 1, 3, and 6 h during the treatment.

The heat and drought stress were performed in a common growth chamber. For heat stress, the chamber was set with temperature of 40°C and humidity of 60%. Leave samples under heat were collected at 1, 3, and 6 h during treatment. For drought stress, the chamber was set as follows: temperature of 25°C and humidity of 40% in 16 h light; temperature of 22°C and humidity of 50% in 8 h dark; withholding water for 7 days, leaf samples were collected at 1, 2, and 3 days during drought treatment.

All the collected leaf samples were immediately frozen in liquid nitrogen, and stored at -70°C for further analysis.

The RNA was isolated from leaf tissues using an RNA extraction kit (Takara, Dalian), according to the manufacturer’s instructions. The first-strand cDNA was synthesized by the Prime Script RT reagent Kit (Takara, Dalian). Real-time quantitative PCR was performed using 2 μl of cDNA in a 20 μl reaction volume with SYBR Premix Ex Taq (Takara) on a 7500-Fast real-time PCR System (Applied Biosystems). Gene-specific primers were designed (Supplementary Table 1). The rapeseed TMA7 gene (BnaC05g11560D) was used as an internal control to normalize the expression level of the target gene, which has highly stable expression level in different tissues and under various growth conditions. Each treatment was repeated three times independently. The thermal cycler was set as follows: an initial incubation at 50°C for 2 min and 95°C for 5 min, followed by 40 cycle at 95°C for 30 s, 55°C for 30 s and 72°C for 30 s. The relative quantification of BnaHsfs transcription levels was determined by the methods described previously (Livak and Schmittgen, 2001).

To identify Hsf genes in B. napus, candidate genes were selected by using the conserved Hsf domain (PF00447) from the Pfam database to query the B. napus genome. Meanwhile, the amino acid sequences of 21 AtHsfs from Arabidopsis were used to protein BLAST the B. napus genome. A total of 64 BnaHsfs (BnaHsf01-BnaHsf64) were identified as members of the Hsf gene family in B. napus, through simultaneous consideration of the conservation of the DBD domain, the coiled-coil structure from the SMART database, and the CD-search of the NCBI database (Table 1).

Fifty-four BnaHsf genes were distributed unevenly among the 19 chromosomes of B. napus from A01 to A10 and C01 to C09; however, 10 members (including BnaHsf02) were located on unanchored scaffolds that could not be mapped to a specific chromosome. Most BnaHsf genes were localized to chromosome A03 and C03 (8 and 7 Hsf genes), while there was only one BnaHsf on chromosome A04 and A06 (Table 1). The deduced proteins of the BnaHsfs ranged from 226 amino acids (aa; BnaHsf56) to 517 aa (BnaHsf19) in length, with predicted isoelectric points (pI) varying from 4.61 (BnaHsf07) to 9.16 (BnaHsf57 and BnaHsf60) and molecular weight (MW) from 26.42 kDa (BnaHsf56) to 58.24 kDa (BnaHsf19) (Table 1).

Among the best-studied 21 Arabidopsis Hsfs, 15 members belong to group A and 5 members belong to Group B, and only one Hsf has been classified as part of group C (Nover et al., 2001). To explore the classification and the evolutionary characteristics of the BnaHsf genes, an unrooted phylogenetic tree was generated using protein sequences of BnaHsfs (Figure 1 and Supplementary Datasheet 1).

According to this phylogenetic analysis, the BnaHsfs genes were grouped into three classes, A (BnaHsf01-BnaHsf43), B (BnaHsf44-BnaHsf60), and C (BnaHsf61-BnaHsf64), as in Arabidopsis. Class A was the largest and consisted of 43 members from eight subclasses (A1–A8). In Arabidopsis, class A is further subdivided into 9 subclasses, A1-A9, with A9 (At5g54070) appearing as a single branch of the AtHsfs molecular phylogeny. However, no orthologous HsfA9 genes were found in B. napus, indicating that Hsf genes in this subgroup were lost. Class B consisted of 17 members from four subclasses (B1–B4), and class C was the smallest, containing only 4 members (Figure 1 and Table 1).

Multiple sequence alignment analysis of BnaHsfs proteins showed that a typical B. napus Hsf protein contains 5 conserved domains, including of DBD, OD, NLS, NES, and AHA domains in order from N-terminal to C-terminal. Among these, the highly structured DBD and OD domains are the most conserved sections in each BnaHsf. The OD domain (HR-A/B region) served as an important basis for classification of BnaHsfs. B. napus class B Hsf proteins, like other plant Hsf proteins, are compact and lack an insertion between the HR-A and HR-B regions (Figure 2), while an insertion of 21 aa in length was found in BnaHsfA and an insertion of 7 aa in length was found in BnaHsfC between the HR-A and HR-B regions. Thus, class A members are less conserved than class B and C members.

To study the structural diversity of BnaHsf genes, the exon/intron organization of individual BnaHsf genes was analyzed by comparing cDNA sequences with the corresponding genomic DNA sequence. The detailed gene structures are shown in Figure 3A. The number of introns ranges from zero to five in BnaHsf genes. Five introns were found in BnaHsf31, while none were found in the BnaHsf63/64 gene pair. Most BnaHsfs contained one or two introns (Figure 3A).

Conserved motif analysis was conducted by using MEME, and 15 motifs were detected in BnaHsf proteins (Figure 3B and Supplementary Figure 1). The DBD and OD domains, composed of motif 4, motif1 and motif 2, were the most conserved and were found in almost all 64 BnaHsf members, while motif 1 was absent in BnaHsf43. The NLS motif (motif 9) and NES motif (motif 15) were found in most class A and B members but not in class C BnaHsf proteins. Motif 3 is the insertion between the HR-A and HR-B regions that was found only in class A and class C members. The AHA motif (motif 6) was found in most class A members but not in classes B and C. In general, the structure of Hsf proteins was conserved throughout the BnaHsfs gene family.

To examine spatial and temporal expression profiles of BnaHsfs across different tissues and organs, an expression pattern map was drawn based on RNA-seq data (Supplementary Table 2) from twelve rapeseed tissues (leaf, root, stem, sepal, silique, pericarp, bud, stamen, ovule, new pistil, mature pistil and wilting pistil).

We found that BnaHsf genes were differently expressed among the subclasses in 14 tissues and at different developmental stages (Figure 4). BnaHsfA6 subclass (BnaHsf29-BnaHsf34) exhibited root specific expression but was hardly detected in other tissues. Subclass A7 (BnaHsf35-BnaHsf40) was also root specific but BnaHsf35 and BnaHsf37 were also detected in reproductive organs such as petal, stamen and pistil. The BnaHsf A2, A3, A4C and A5 subclasses were the most abundant BnaHsf genes and were constitutively expressed among all tissues, as were BnaHsf23 and -24 in subclass A4A, BnaHsf46 and -47 in subclass B1, and BnaHsf48 and -49 in subclass B2A. Subclass B3 (BnaHsf55-56) showed a higher transcription level in root, sepal and bud tissues. BnaHsf55 and -56 in subclass B4 were specifically expressed in ovule tissues. In contrast, class BnaHsfC was inactive in ovule tissues.

To determine the potential role BnaHsfs play in plant responses to different environmental stresses, the expression levels of BnaHsf genes under high temperature, drought and high CO₂ stresses were analyzed using RNA-seq data (Supplementary Table 3). The mRNA for these transcriptomic analyses were extracted from the leaves of rapeseed plants both in normal growth conditions and after 3 h of heat treatment, 3 days of drought treatment and 3 h of high CO₂ treatment. The results showed that BnaHsf genes were very sensitive to heat and drought stress (Figure 5). Members of subclasses A2, A4C, A5, and class B (except for the non-expressed subclasses B3 and B4) showed relatively higher basic transcription levels in leaf tissue under normal growth conditions. Among the subclasses with high basic expression levels, BnaHsf15 and -16, in subclass A2, were dramatically upregulated (>25-fold), becoming the predominant transcripts after 3 days of exposure to drought stress, but were suppressed after 3 h of heat treatment. B1 members, except for BnaHsf45, were strongly induced by drought stress, and BnaHsf46 and -47 were strongly induced by heat. All B2 members were significantly upregulated under drought conditions, but were only slightly induced by 3 h of heat treatment. As observed in the subclasses with low basic expression, a moderate induction was seen in the A1E subclass after exposure to drought and heat. Members of subclasses A3 and A4A were strongly induced by both drought and heat stress. Strikingly, the highest induction (>350-fold for BnaHsf36∼38) by drought was observed in A7A subclass, although expression of members of this subclass was hardly detectable under normal conditions. Heat treatment also resulted in a marked induction in A7A members. Members in class C were only upregulated by drought treatment. In the case of individual member genes, BnaHsf07 (in A1D) and BnaHsf40 (in A7B) were induced by drought stress, and BnaHsf43 in subclass A8 was strongly induced by heat stress.

Unlike drought and heat stress, high CO₂ treatment did not cause a significant effect on the transcription level of most of the BnaHsf family genes (Figure 5). However, three member genes, BnaHsf18 (in A3), BnaHsf21 (in A4A) and BnaHsf43 (in A8), were clearly induced (3 to 9-fold) by exposure to high CO₂ conditions. However, the expression of members of subclass A2 members (BnaHsf15/16) and of subclass B2 was largely suppressed by high CO₂ treatment.

Twelve BnaHsf genes from three main classes were selected for examination of their function under three stress conditions using quantitative Real-Time PCR (qRT-PCR). These genes included BnaHsf04 from subclass A1, BnaHsf15 and -16 from subclass A2, BnaHsf17 and -18 from subclass A3, BnaHsf20 and BnaHsf23 from subclass A4, BnaHsf42 from subclass A8, BnaHsf45 from subclass B1, BnaHsf46 from subclass B2 and BnaHsf61 from class C. qRT-PCR was carried out using rapeseed plants exposed to heat (0, 1, 3, and 6 h), drought (0, 1, 2, and 3 days), and high CO₂ (0, 1, 3, and 6 h).

The results of the qRT-PCR analyses were consistent with the expression patterns of selected BnaHsfs from RNA-seq data, and provided more details under progressive stresses (Figure 6 and Supplementary Table 4). The RNA-seq data showed that transcription levels of BnaHsf15-16 were decreased after exposure to 3 h heat (Figure 5), while in qRT-PCR results they were upregulated after 1 h of heat, and their expression then dropped to lower than basic levels (Figure 6). However, the strong induction of BnaHsf15-16 began with a marked reduction in transcription after 1 day of treatment under drought stress. In contrast to BnaHsf15-16, the expression levels of BnaHsf17-18 and BnaHsf43 were progressively induced by all prolonged drought, heat, and high CO₂ stress conditions. A similar pattern was found in BnaHsf43 under drought and heat conditions. For BnaHsf21-22, the induction in their transcription was weakened as heat stress processed, while under progressive drought, the induction by stress showed an opposite enhancement (Figure 6).

To identify the presence of putative regulatory cis-acting elements enriched in BnHsf genes, the promoter sequences upstream of their CDS were extracted and searched against the PLACE database (Higo et al., 1999). Analysis showed that HSEs were the most abundant cis-elements in promoter regions of BnHsf genes, including perfect type and active HSE variants (Table 2). Many of the subclasses of the BnHsf gene family possess the two types of HSE, although these are not present in subclasses A3 and A6A. Subclass A6A contains only one member, BnHsf29, which may be a pseudogene since its transcription cannot be detected in any tissues of rapeseed plants. The two members (BnHsf17 and -18) of subclass A3 were found to have three other major types of stress-related cis-elements present in the BnHsf family, STRE, DRE/CRT, and MYCATRD22. The STRE element was first found to be stress responsive in yeast and can serve as a specific binding site for HsfA1a in Arabidopsis (Martinez-Pastor et al., 1996; Haralampidis et al., 2002; Guo et al., 2008). STRE was present in most subclasses of BnaHsf except A1A, A6A, and A7. DRE/CRT and MYCATRD22, two types of cis-elements responsive to drought stress, also appeared in most BnHsf subclasses. Two other stress related cis-elements, ABREOSRAB21 and LTRE, were found in some of the BnaHsf family members. In addition, three CO₂-responsive elements (CCRE1/2/3) (Ohno et al., 2012; Tanaka et al., 2016) were observed in all BnHsf subclasses except A1A/B/D, A4A, A5, A6B, A7B, and B2A. The presence of these stress related cis-elements is likely responsible for the regulative expression patterns of BnHsf genes under drought, heat, and high CO₂ conditions. Moreover, some phytohormone responsive related cis-elements were enriched in promoter regions of BnHsf member genes, which may be involved in the stress acclimation response and development.

Brassica napus (rapeseed, genome AACC) is an amphidiploid species formed by recent interspecific hybridization between ancestors of B. oleracea (genome CC) and B. rapa (genome AA) (Chalhoub et al., 2014). In this study, we identified 64 Hsf genes in the genome of B. napus. Unlike yeast and animals, plants usually have many Hsf coding genes. There are 21 Hsf member genes in the model plant Arabidopsis thaliana (Nover et al., 2001), 25 members in rice (Chauhan et al., 2011), 56 Hsf genes in wheat (Xue et al., 2014). To date, BnaHsfs represent the largest Hsf gene family in plant species of which Hsf member genes were analyzed. The diversification of plant Hsfs is presumably the result of gene- and whole-genome duplications (WGD) at different points in evolution, followed by gene loss (Chalhoub et al., 2014). In the case of rapeseed, the allopolyploid process, that followed from the fusion of genomes A and C, might also play a crucial role in the expansion of the BnaHsf gene family. In addition, the large size of the BnaHsf family may have been needed for adaptation of rapeseed to diverse climatic zones.

Similar to other plant Hsf families, the modular structure of rapeseed Hsf proteins is well conserved. While in comparison with that in Arabidopsis, there is no Hsf member in subclass A9 subclass in rapeseed. This differs from other eudicot plants, most of whom have sub classA9 Hsfs. This subclass is also lost in the Hsf family of B. rapa (Huang et al., 2015). The DBD is characterized as a central domain for the Hsf protein: it specifically binds to HSEs in the target promoter region, and subsequently activates the transcription of associated heat-inducible genes. The DBD domain of plant Hsfs is encoded by two regions separated by an evolutionarily conserved intron, which was inserted immediately adjacent to the HTH DNA binding motif (Nover et al., 2001; Scharf et al., 2012). Most BnaHsf genes have this intron in their DBD domains; however, no intron was found in BnaHsf63 and -64 genes from class C, as shown by their gene structure (Figure 3A). As far as we know, the fact that is this highly conserved intron in not present in the DBD domain of a plant Hsf is unique. Furthermore, this fact may indicate that BnaHsf63 and -64 have a novel regulation pattern relative to other Hsf genes.

The functional diversification of BnaHsf family genes was also found in the tempo-spatial expression profile of these genes during development and abiotic stress treatments. Among the tissues at different developmental stages, subclasses A1A, A2/3, A4C, A5, A8, and B2 were found to be constitutively expressed at relatively high levels in all the tissues examined. While almost all member genes from subclasses A6/7 and B3/4 were hardly detected in any tissue, BnaHsf35 from subclass A7A showed a high level of expression in root tissue. Subclasses B1 and B2A also showed high levels of expression in root tissue. Members of class C also showed increased expression in root tissue, but were not expressed in ovule tissue. These results suggest that these BnaHsf genes may be involved in root development. Under abiotic stress, many BnaHsf genes were upregulated in response to drought treatment. The number of drought induced BnaHsf genes was comparable to those induced by heat. This suggests that Hsf genes may also play an important role in the response and the acclimation to drought stress in rapeseed. Furthermore, the most inducible BnaHsf genes were upregulated by both drought and heat treatment, as shown by the combination of RNA-seq and qRT-PCR data. While BnaHsf43 of subclass A8 was only induced by heat, BnaHsf07 of subclass A1D and members of class C were predominantly upregulated by drought.

According to the combined transcriptional analysis, heat inducible BnaHsf genes could be divided into three groups. The first group consisted of BnaHsf15 and -16 (subclass A2), BnaHsf47 (subclass B1) and BnaHsf50 (subclass B2A), in which the expression of member genes exhibited an immediate and strong induction after 1 h of heat to a high level of 40∼70-fold of that in non-stressed control, followed by a dramatic drop to the basic expression level after 3 and 6 h of heat treatment, even was slightly suppressed in subclass A2 after 3 and 6 h of heat. It may be that this group of genes governs early heat stress response. The second group contains BnaHsf21 and -22 of A4A, and BnaHsf46 of subclass B1. The transcription of the second group members showed also a fast and strong upregulation after 1 h of heat exposure. This upregulation of gene expression gradually declined after 3 and 6 h of heat, but still maintained at a high level relative to that in control. Genes from the second group might be involved in both early and late heat response. The last group greatly differed from the other two, comprised of BnaHsf17 and -18 of subclass A3 and BnaHsf43 of subclass A8. The genes from this third group were upregulated after 1 h of heat, and this induction was continuously enhanced as the stress treatment proceeded, finally peaking after 6 h of heat stress (46∼550-fold vs. control). The members of this group likely have some function to facilitate acclimation to prolonged heat stress.

Drought induced BnaHsf genes seemed to have a single expression pattern. The genes continuously increased transcription during exposure to drought, and reached peak expression after the 3 days drought treatment. BnaHsf genes also played an important role in response to high CO₂ treatment, as BnaHsf18 (A3), BnaHsf21 (A4A) and BnaHsf43 (A8) were strongly upregulated, while members of A2 and B2 subclasses were downregulated.

Regulatory element analysis revealed that there were many stress-related cis-acting elements enriched in the promoter regions of BnaHsf family genes. HSEs were found to be the dominant cis-elements (Table 2). Complex interactions may exist among BnaHsf genes, and these may come about via trans-acting regulation, as HSEs are marker binding sites for plant Hsf proteins. Previous work has supported this idea that HsfA1a/b target class B Hsf genes and are responsible for their induction during heat response in Arabidopsis (Lohmann et al., 2004; von Koskull-Döring et al., 2007), and that HsfA5 inhibits the activity of HsfA4 (Baniwal et al., 2007). Other abiotic stress-related cis-elements, including STRE, DRE/CRT, MYCATRD22, ABRE, CCRE, and LTRE were also major regulatory elements found in BnaHsf genes. The presence of these stress-related elements seemed to be correlated to the expression response of BnaHsfs to heat, drought, and high CO₂ treatments. For example, many drought related DRE/CRT and MYCATRD22 elements found upstream were associated with a marked induction of BnaHsf genes by drought stress. The two CCRE elements situated in promoter region also agreed with our observation of high induction levels of BnaHsf18 under high CO₂ treatment. Unlike BnHsf15 and -16 of subclass A2, the heat highly inducible genes BnHsf17 and -18 of subclass A3 do not have functional HSE elements, but rather an STRE element was found upstream of the target genes. The STRE element was identified to be stress responsive, and serves as a direct binding site for HsfA1a besides HSE in Arabidopsis (Guo et al., 2008). Furthermore, the deletion of STRE from the promoter of the Arabidopsis Hsp90-1 gene decreased its promoter activity under heat stress conditions (Haralampidis et al., 2002). These findings suggest that STRE also plays a crucial role in transcriptional regulation under heat conditions, as do HSEs. Moreover, the different heat induced expression patterns of subclass A3 and subclass A2 BnaHsf genes provides evidence for differential transcription regulation abilities of STRE and HSE element. Unexpectedly, rapeseed subclass A1A was not heat inducible, although HSE elements are found in the promoter. While HsfA1a serves as master regulator of thermotolerance in tomato (Mishra et al., 2002), it also functions actively in Arabidopsis under heat stress. These results may indicate differential gene regulation of rapeseed Hsf genes from those found in other plants, even those in the Brassicaceae.