Identification of bHLH family genes in barley
Proteins containing bHLH domain in barley were identified by searching the whole barley genome. Totally 168 bHLH genes were obtained in barley genome after removing the redundant protein sequences (Additional file 1 and 2: Table S1 and S2). By EXPASY analysis, we found the amino acid length of HvbHLH proteins ranged from 36 to 938, accordingly the molecular weight (MW) ranging from 9.63 to 102.01 kDa (Additional file 1: Table S1). The isoelectric point (pI) varied greatly, from 4.63 to 11.9 (Additional file 1: Table S1). While the MW and pI of 11 HvbHLH proteins were not predicted, due to existence of some consecutive undefined amino acids (Additional file 1: Table S1).
Phylogenetic analysis of bHLH proteins
In order to analyze the evolutionary relationship of the bHLH genes between barley and other species, the unrooted phylogenetic tree was constructed using the full-length amino acid sequences of all bHLH in Arabidopsis, rice and barley by the neighbor-joining (NJ) method. The result showed that 450 bHLH proteins from Arabidopsis, rice and barley were divided into 27 subgroups (Subgroup 1 to 27; Fig. 1) based on the classification of bHLH proteins in Arabidopsis and rice [6-7]. Except Subgroup 20, all members of bHLH family in barley could be found in all the other subgroups (Fig. 1). Subgroups 3 and 9 were the smallest, each only having three members in barley, rice and Arabidopsis. Subgroup 27 was the largest, having 20, 16 and 17 members in barley, rice and Arabidopsis, respectively. Another NJ phylogenetic tree was also constructed based on HvbHLH protein sequences, and HvbHLHs were classified into 25 subfamilies (Subfamily A to Y, Fig. 2a). Classification in HvbHLH was generally consistent with those in Arabidopsis and rice phylogenic tree (Fig. 1 and 2a).
Analysis of motifs and exon/intron structures
Ten conserved motifs containing 11 to 57 amino acids were predicted in HvbHLH proteins through MEME software (Additional file 3: Fig. S1). HvbHLH proteins consist of 1- 5 motifs, varying with different subfamilies (Fig. 2b). However, all of them contained at least of the Motif 1 or 2, and the proteins in the same subfamily exhibited the similar motif composition (Figs. 2a and b). For example, Motifs 2, 3, and 10 were identified in all 6 members of Subfamily K, and Motifs 1, 2, and 7 were identified in 12 of the 15 members of Subfamily X. Furthermore, some motifs are unique to one or more Subfamilies. For instance, Motif 5 was specifically shared by each member in the Subfamily D, and motif 3 was only found in Subfamilies B, K and L. This conservation of the motif composition patterns in each subfamily might indicate their similar biological functions in the same subfamily [31].
The alignment of genome and coding sequences, and also exon/intron structure within the coding sequence in bHLH genes was performed using the TBtools [32]. The structure differed greatly among HvbHLH genes, but most members of the same subfamily had the similar exon/intron structure (Figs. 2a and c). For instance, the members in the Subfamily A had the relatively simple structure, possessing only 1-2 exons, while the members in the Subfamily T had 6-8 exons, being much more complex in its structure (Fig. 2c).
Chromosomal distribution and gene duplication of HvbHLH genes
According to their physical positions, 161 HvHLH genes were nonrandomly mapped on seven barley chromosomes, with most located at the both ends of the chromosomes (Fig. 3). More than one-third of the HvbHLH genes were located on Chr3 (30 members) and Chr5 (26 members), and only 12 genes was located on Chr1 (Fig. 3). In addition, 7 HvbHLHs could not be conclusively localized to a chromosome (Fig. 3).
Totally 18 segmental duplication events were identified among HvbHLHs (Fig. 4, Additional file 4: Table S3). Among them, 16 events happened between diverse chromosomes, while only 2 events were detected within a same chromosome. Nine genes were identified as tandem duplication genes (Fig. 4, Additional file 4: Table S3), and the genes located on chromosome 3 showed the highest numbers (7) of tandem duplications, while two genes (HORVU4Hr1G087590 and HORVU4Hr1G087610) were tandem duplicated on Chr4 (Fig. 4). Each pair of duplication genes were derived from the same subfamily (Fig. 1 and 4). Obviously, some HvbHLH genes were generated due to gene duplication, and the segmental duplication may play vital roles in the HvbHLH gene family expansion in barley genome.
Stress-related cis-elements in the promoter of HvbHLH genes
To understand the potential regulatory mechanisms of HvbHLH genes, cis-elements were analyzed using 1.5-kb upstream sequences from ATG through PlantCARE. Various cis-elements involved in abiotic stress responses were identified, and mainly they could be divided into two groups. One is elements responsive to hormones, including abscisic acid (ABRE), jasmonic acid (CGTCA and TGACG motif), auxin (AuxRR-core and TGA-element), gibberellins (GARE motif, TATC-box and P-box,), salicylic acid (SARE and TCA-elements). The other one is responsive to environmental stress, including low-temperature responsive element (LTR), anaerobic induction element (ARE), MYB binding site involved in drought-inducibility (MBS), and defense and stress responsive elements (TC-rich repeats)(Fig. 5, Additional file 5: Table S4). All HvbHLH genes had light responsive cis-elements, belonging to the most commonly predicted cis-elements in the promoter of HvbHLH genes, and 155 (92.26%) HvbHLH genes contained the G-box elements.most genes had more than one kind of cis-elements, suggesting that both groups could respond to multiple stresses.
Expression profiles of HvbHLHs in different plant tissues
To further understand the functions of HvbHLHs, the expression profiles of HvbHLH genes in 8 plant tissues of a barley cultivar ‘Morex’ were analyzed based on the available transcriptomic data [33]. The expression of 149 HvbHLH genes were identified in all 8 tissues, and their expression levels varied considerably with genes and tissues (Fig. 6 and Additional file 6: Table S5). Some HvbHLHs in the same subfamily showed the similar expression pattern, while some others displayed the different patterns (Fig. 6). Totally 128 (128/149, 85.91%), 125 (125/149, 83.89%), 118 (118/149, 79.19%), and 101 (101/149, 67.79%) HvbHLH genes expressed in four-day-old embryos, roots and shoots, and internodes of six-leaf-old seedlings, respectively (Fig. 6). In addition, 118 (118/149, 79.19%) and 110 (110/149, 67.79%) HvbHLH genes expressed in young fluoresce (5 mm and 1-1.5 cm) and developing grains (5- and 15-day post anthesis, DPA), respectively (Fig. 6).
Some HvbHLHs highly expressed in the specific tissue(s). For instance, the genes HORVU2Hr1G044230 and HORVU6Hr1G020520 showed higher expression in caryopsis (5 and 15 DPA) than other tissues (Fig. 6), suggesting their significant roles in grain growth and development. HORVU2Hr1G044230 and HORVU1Hr1G072810 had higher expression in embryos and roots, while HORVU2Hr1G044230, HORVU7Hr1G050530 and HORVU1Hr1G054260 had higher expression in shoots (Fig. 6), indicating these genes may participate in growth and development of vegetative organs. Notably, HORVU2Hr1G044230 displayed extremely high expression in all the tissues (Fig. 6).
Expression profiles of HvbHLHs under low nitrogen stress
The expression patterns of HvbHLH genes in response to LN stress were analyzed in the roots of the two barley genotypes (XZ149, LN tolerant and XZ56, LN sensitive) using RNA-seq (Additional files 7 and 8: Tables S6 and S7). Consequently, 16 differentially expressed genes (DEGs) encoding HvbHLH proteins were identified using pair-wise comparison for each genotype under LN stress (Fig. 7, Additional file 9: Table S8). Remarkably, no DEG was found in the two barley genotypes under LN stress at 12 d, suggesting HvbHLHs may play a regulating role at early stress. Furthermore, all the DEGs were changed only at one time point, except MLOC_74557 (HORVU6Hr1G064820) and MLOC_21066 (HORVU3Hr1G000150), which was down-regulated in XZ149 both at 6 h and 48h after stress (Fig. 7).
Eleven DEGs were down-regulated in their expression under LN stress, and their response patterns showed the marked difference. Six DEGs responded quickly to LN stress (at 6 h), and MLOC_37666 (HORVU2Hr1G108480), MLOC_72940 (HORVU2Hr1G117820), MLOC_21066 (HORVU3Hr1G000150) were down-regulated only in XZ149, while MLOC_67937 (HORVU5Hr1G002090) was down-regulated only in XZ56 (Fig. 7). The other five DEGs were changed at 48 h after LN stress, and among them MLOC_12832 (HORVU5Hr1G066530), MLOC_66385 (HORVU4Hr1G061760) were down-regulated only in XZ149, MLOC_64390 (HORVU3Hr1G027630) only in XZ56, while MLOC_62844 (HORVU2Hr1G066100) and MLOC_43080 (HORVU3Hr1G000150) in both XZ149 and XZ56 (Fig. 7).
Four DEGs were up-regulated under LN stress, and three of them responded at 6 h, with MLOC_72280 (HORVU2Hr1G114070) only in XZ149, and MLOC_36351 (HORVU3Hr1G108680) and MLOC_76486 (HORVU5Hr1G018100) in both XZ149 and XZ56. MLOC_55964 (HORVU4Hr1G075250) was up-regulated at 48 h in both XZ149 and XZ6(Fig. 7). Notably, the change fold of all the up-regulated DEGs was larger in the tolerant genotype XZ149 relative to the sensitive genotype XZ56 (Fig. 7).
To confirm the transcriptomic results and understand the expression patterns of HvbHLHs, the dynamic response to LN stress were examined using real time PCR for the 6 identified DEGs at 8 time points after LN stress. The expression patterns of HvbHLHs at 6 h, 48 h and 12 d after LN stress were highly consistent with those obtained by the transcriptomic analysis. (Additional file 10: Fig. S2). Under LN stress, different HvbHLHs exhibited distinct dynamic expression patterns, suggesting that HvbHLHs may be involved in regulation of LN tolerance in different ways. For instance, MLOC_72280 (HORVU2Hr1G114070) was up-regulated during 1 h -2 d under LN stress, but down-regulated at 3 d -12 d in XZ149, while it was up-regulated at 1 h, 3 h, 6 h, and 2 d, and down-regulated at 24 h and 3 d - 12 d in XZ56 (Additional file 11: Fig. S3). For MLOC_37666 (HORVU2Hr1G108480), its expression generally showed the decrease with the time of exposed LN stress in XZ149, but increase at 1 h, 3 h, and 6 d under LN stress in XZ56 (Additional file 11: Fig. S3).