We applied the genome project (http://prunusmumegenome.bjfu.edu.cn) to download the whole genome data of P. mume (Zhang et al., 2012). The sequences of A. thaliana bHLHs (AtbHLHs) were obtained from Pires and Dolan (Pires and Dolan, 2010). The Pfam database (http://pfam.xfam.org, PF00010) was used to obtain the Hidden Markov Model (HMM) profile of the HLH domain and was searched bHLH proteins of P. mume with HMMER3 software (http://hmmer.janelia.org) (Sun et al., 2015). To ensure credibility, the E-value cut-off was set at 10^–5. SMART online software (http://smart.embl-heidelberg.de/) (Schultz et al., 1998) and NCBI CD-search (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (Marchler-Bauer et al., 2015) were used to confirm bHLH domains in presumptive PmbHLH proteins of P. mume.

The WoLF PSORT program (https://wolfpsort.hgc.jp/) was applied to predict the subcellular localization of PmbHLHs (Horton et al., 2007). The CDS length, molecular weights (MWs), theoretical isoelectric points (pI), amino-acid sequences (aa), the total number of positively charged residues (Arg + Lys), total number of negatively charged residues (Asp + Glu), grand averages of hydropathicity (GRAVYs), and instability index and aliphatic index of all predicted PmbHLHs were calculated using ExPASy (https://web.expasy.org/protparam/) (Artimo et al., 2012).

RAxML version eight software with maximum likelihood (ML) method was used to construct phylogenetic trees (Stamatakis, 2014). The optimal JTT (Jones-Taylor-Thornton) model of amino acid substitution was applied to construct ML phylogenetic trees. The calculated relationships of the phylogenetic tree were supported by performing 1,000 iterations of bootstrap test and visualized by iTOL (https://itol.embl.de/) (Letunic and Bork, 2006). Alignments of PmbHLHs domains were conducted by Muscle (https://www.ebi.ac.uk/Tools/msa/muscle/) (Edgar, 2004).

TBtools (Chen et al., 2020) and NCBI Batch CD-Search (https://www.ncbi.nlm.nih.gov/) (Marchler-Bauer et al., 2015) were applied to analyze and visualize gene structure and conserved domains. We applied the online MEME program to analyze motif structures of PmbHLH proteins (30 motifs were set as the maximum number), others using default parameters (Bailey et al., 2009). Jalview software (Clamp et al., 2004) and Weblogo3 (http://weblogo.berkeley.edu/logo.cgi) were used to visualize and analyze conserved domains.

We used the GDR database (https://www.rosaceae.org/) to retrieve genomes of Prunus avium and Prunus persica (Verde et al., 2013; Shirasawa et al., 2017). The P. mume genome database provided chromosomal distribution information of bHLH genes. Chromosomal location map, as well as duplication events of bHLH genes and mutation rates of Ka (nonsynonymous) versus Ks (synonymous) were predicted through TBtools (Chen et al., 2020). Syntenic relationship of bHLH genes was analyzed by MCScanX (Multiple Collinearity Scan toolkit) in P. mume, P. avium, and P. persica (Wang et al., 2012). We used the R circlize package to visualize the relationship between three varieties (Zhuo et al., 2018). The PmbHLH genes divergence time (T) was calculated through the equation: T = dS/(2λ × 106) Mya, where λ = 1.5 × 10⁻⁸ s for dicots (Lynch and Conery, 2000).

The transcriptome data of PmbHLHs in five tissues (bud, fruit, leaf, stem, and root) were downloaded from the NCBI Sequence Read Archive (accession number: SRP014885). Heat map of five tissues and low-temperature stress were drawn by TBtools (Chen et al., 2020) and iTOL (https://itol.embl.de/) based on FPKM values (Letunic and Bork, 2006).

Annual grafted seedlings of true mume ‘Beijing Yudie’ and apricot mei ‘Danfenghou’ were used for low-temperature treatment. These plant materials were exposed to the same soil humidity (60–70%), air humidity (65–70%), illumination time (12 h/12 h), and light intensity (150 μmol•m⁻²•s⁻¹). Stems of P. mume were treated for chilling treatment (4°C for 8 h and 5 days) and freezing treatment (−5°C for 1 h).

To further research the regulation process of PmbHLHs under low-temperature, we carried out detailed multi-stage chilling treatment and freezing treatment.

1) Chilling treatment: in September, we took annual branches of ‘Beijing Yudie” and put branches in a 4°C low-temperature incubator and undertook sampling at 0/2/4/8/16/24/48/72 h. Stem segments without bud points were used as sampling materials.

Total RNA of low-temperature treatments was isolated by RNA Extraction Kit (Takara, Beijing, China) based on references. First-strand cDNA synthesis was reversed with DNase-treated RNA (1 µg) through PrimeScript^TMRT Reagent Kit with gDNA Eraser (Takara, Dalian, China). The template used cDNA (2 µL) in a 10 µL qRT-PCR by TB Green Ⅱ Premix Ex Taq (Takara, Dalian, China). We applied the 2^ΔΔCt’ method to calculate relative expression levels and internal control used protein phosphatase 2A (PP2A)gene of P. mume (Wang et al., 2014). Supplementary Table S10 showed ten selected genes and specific primers. Each qRT-PCR was repeated at least three times.

Promoter cis-acting regulatory elements were analyzed by PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). The AraNet V2 tool (Lee et al., 2015) was used to construct the protein interaction network based on homologous proteins of PmbHLHs in Arabidopsis. The protein interaction network was visualized by Cytoscape (Shannon, 2003) and STRING software (http://string-db.org/).

A total of 95 non-redundant PmbHLHs were discovered based on HMMER software from P. mume genome. According to their location position, these genes were named PmbHLH01 to PmbHLH95. The length of PmbHLH proteins ranged between 91 (PmbHLH36) and 700 (PmbHLH45) amino acids and most of these genes (74%) had lengths of 200–400 aa. The presumptive isoelectric points (pI) ranged from 4.57 (PmbHLH14) to 10.1 (PmbHLH43). About molecular weight values, the smallest was 10.27 kDa (PmbHLH67) and the largest was 78.26 kDa (PmbHLH45). The instability index varied from 36.78 to 92.34, while only one (PmbHLH66) was considered a stable protein. Predicted GRAVY values were all negative, representing all genes that possessed hydrophilic characteristics. The aliphatic index, (Asp + Glu) value and (Arg + Lys) value of PmbHLHs showed diversity features. Analysis of the gene structure of 95 PmbHLH proteins showed that coding genes of most PmbHLH proteins (92.63%) have introns. Subcellular localization of PmbHLH proteins were mostly located in the nucleus (Supplementary Table S1).

To research the evolutionary relationship of PmbHLHs, a phylogenetic tree with ML method was established by full-length amino acid sequences of 95 PmbHLHs (Figure 1). To clearly describe the classification and potential functions of PmbHLHs, we performed another ML phylogenetic tree with 95 PmbHLH proteins and 162 AtbHLHs (Supplementary Figure S1 and Supplementary Table S2). Based on verified AtbHLHs (Pires and Dolan, 2010), 95 PmbHLH proteins were divided into 23 subfamilies and one orphan (Figure 1). XII was the largest subfamily, consisting of 11 members. Five subfamilies [IVc, Va, VIIIa, VIIIc (1), XIV] contained only one bHLH protein. This indicated that PmbHLHs were distributed unevenly in different subfamilies. To date, we know very little about the biological functions of PmbHLHs except for PmICE1 (Cao et al., 2014). However, plenty of bHLH proteins have been confirmed in Arabidopsis. Through Arabidopsis ortholog analysis, PmbHLHs might play an important role in abiotic stress, hormonal regulation, development, and so on (Supplementary Table S2).

The bHLH domain analysis showed that PmbHLH domains consisted of four conserved regions, including the basic region, two helix regions, and a loop region (Figure 2 and Supplementary Figure S2). Additionally, 23 amino acid residues were conserved (>50% consensus ratio) in their bHLH domains. Among residues, the Arg-15, Arg-16, Leu-26, Pro-31, Leu-61 were highly conserved with a greater than 90% consensus ratio. Among the 23 conserved amino acid residues, the basic region found six conserved residues (His-8, Ala-11, Glu-12, Arg-13, Arg-15, Arg-16), first helix region found seven conserved residues (Ile-19, Asn-20, Arg-22, Leu-26, Leu-29, Val-30, Pro-31), loop region found two conserved residues (Lys-44, Asp-45), second helix region found eight conserved residues (Ala-48, Ser-49, Leu-51, Ala-54, Ile-55, Tyr-57, Lys-59, Leu-61).

The DNA binding activity of target genes was decided by the bHLH domain in the basic region. According to the classification standard in A. thaliana (Toledo-Ortiz et al., 2003), PmbHLHs were defined as Non-DNA-binding proteins and DNA-binding proteins (Figure 2 and Table 1). In addition, DNA-binding proteins were divided into E-box-binding proteins (including G-box-binding proteins) and non-E-box binding proteins according to the existence of Glu-12 and Arg-15 (positions were corresponding to positions 13 and 16 in Arabidopsis). His/Lys-8, Glu-12. Arg-16 (positions were corresponding to positions 9, 13, and 17 in Arabidopsis) are responsible for the binding of the G-box. Based on the conservation of these residues, three PmbHLHs were classified to non-E-box-binding proteins for missing Glu-12/Arg-15 residues and 37 PmbHLHs were classed to E-box-binding proteins. Among 37 E-box-binding proteins, 24 as G-box-binding proteins, while 13 proteins missed the G-box-binding site. Furthermore, 55 of the 95 PmbHLHs were classed in non-DNA-binding proteins due to less than six amino acid residues in the basic region (Table 1 and Supplementary Table S3).

Further to the analysis features of PmbHLHs, we investigated intron/exon patterns according to the phylogenetic tree with PmbHLHs sequences (Figure 3A). Analysis of genomic DNA sequences showed that number of introns changing from zero to ten (Figure 3B). Most of them usually had one to eight introns, except PmbHLH69 and PmbHLH92. Different subfamilies had different intron/exon patterns, while the same classes were similar. Seven genes of 95 PmbHLHs were intron-less, five of these were in subfamily VIIIb. Seven subfamilies [Ib (1), IIIf, IVb, Vb, VIIIb, IX, XV] had a concentrated number of exons/introns.

We used MEME online tool to predict thirty conserved motifs of 95 PmbHLH proteins (Figure 3C and Supplementary Table S4). The number of PmbHLHs motifs was distinctive, ranging from 1 to 16. Most of PmbHLHs shared three to six motifs. Each PmbHLH protein contained motif one and motif two except for PmbHLH31, PmbHLH58, and PmbHLH69. Different subfamilies had unique motif combinations. For instance, subfamily IVb contained motifs 1, 2, and 26, and subfamily Ib (1) contained motifs 1, 2, 5, and 8. Some motifs were unique and existed in only one subfamily. Subfamily Ⅲ (d + e) contained the most motifs and motif 9, 20, and 25 were specific to it. Motif 11, 23, and 30 were respectively observed in subfamily XI, Ⅷ (2), and XI. These conserved motifs may play special functions. The same subfamily had similar motifs, implying these PmbHLHs might have similar functions. In contrast, the kinds of motifs showed little difference in the same subfamily. For example, PmbHLH07 had motif 21 except common motifs 1, 2, 5, and eight in subfamily IVa.

According to genome annotation information, 78 PmbHLH genes were mapped on eight chromosomes, while 17 PmbHLH genes were localized on unassembled genomic scaffolds (Figure 4A and Supplementary Figure S3). The distribution of PmbHLH genes on each chromosome was irregular. 18 PmbHLH genes (18.95%) were present on chromosome 2, which was the maximum number, whereas only two PmbHLH genes (2.11%) were located on chromosome 3. The proportion of PmbHLH genes on half of the chromosomes accounted for more than 10% but less than 20%.

Duplication events contained genome duplication, tandem duplication, segmental duplication, and transposon duplication, leading to plant evolution (Qiao et al., 2019). Duplication played an important role in PmbHLH gene expansion. Among 95 PmbHLHs, nine pairs of PmbHLH genes were described as tandem duplication, mapping on chromosome 1, 5, 6, 7, and scaffolds (Figure 4A and Supplementary Figure S3). The largest number of tandem duplications were distributed on Chromosome 6. By contrast, segmental duplication was absent in the PmbHLH gene family, which signified segmental duplication was not involved in gene expansion. The selection pressure of gene duplications was estimated by mutation rates of Ka (nonsynonymous) versus Ks (synonymous). We calculated the Ka/Ks of the PmbHLH gene family (Supplementary Table S5). Results showed that most of the Ka/Ks values of tandem duplication were <1 and changed from 0.270 to 0.949, implying a purifying selection during PmbHLH genes expansion. However, PmbHLH82-PmbHLH83 and PmbHLH85-PmbHLH86 only existed Ka value, which means that two gene pairs evolved by natural selection. The divergence time of PmbHLHs tandem duplication ranged from 3.67 to 64.06 Mya.

To further explore the evolutionary mechanism of the PmbHLH gene family, we constructed a syntenic map of P. mume associated with P. avium and P. persica. 41 syntenic orthologous gene pairs were distinguished between P. avium and P. mume, 59 pairs between P. persica and P. mume (Figure 4B and Figure 4C). This indicates that the P. avium, P. persica, and P. mume had a close relationship. Interestingly, we found that one PmbHLH gene only corresponds to one gene in these syntenic orthologous gene pairs. For further evolutionary studies, the divergence time of bHLH gene pairs was calculated in three varieties (Supplementary Table S6). Divergence time started 86.11 Mya to 0.67 Mya between P. persica and P. mume, and 1-2 Mya occurred in most duplicated events. In P. avium and P. mume, it began 60.81 Mya to 0.14 Mya and 0.5-2 Mya occurred in most of the duplicated events, which may indicate that the speciation time of these orthologous pairs was shorter in the two varieties. Additionally, higher syntenic genes appeared on chromosomes 2, 6, and 7.

Based on the FPKM values from RNA sequencing data, we investigated the expression pattern of PmbHLHs among different tissues and different low-temperature treatments. The heatmap of Figure 1 and Supplementary Figure S5 show that the PmbHLH gene family presented clear tissue-specific expression. Among 95 PmbHLHs, 53 genes were expressed in five tissues (bud, fruit, leaf, stem, and root), implying that these PmbHLHs may participate in the development and growth process of tissues. While, PmbHLH21 and PmbHLH32 (they belonged to subfamily XV and VIIa, respectively) lacked expression in all detected tissues. We discovered that nine PmbHLHs were expressed in only one tissue. They included two genes (PmbHLH54 and PmbHLH88) only in fruit, seven genes (PmbHLH02, PmbHLH05, PmbHLH19, PmbHLH42, PmbHLH43, PmbHLH53 and PmbHLH94) only in root. This probably indicated that the nine genes have special functions in fruit or root. Moreover, most of remaining PmbHLHs were expressed in three or four detected tissues (Supplementary Table S7). According to the phylogenetic tree analysis, PmbHLHs from the same subfamily had a similar expression pattern (Figure 1). For instance, IX subfamily (included PmbHLH01, PmbHLH13, PmbHLH25, PmbHLH38 and PmbHLH72) exhibited expression in all detected tissues. Furthermore, Ib (2) subfamily genes (included PmbHLH42, PmbHLH43, PmbHLH65, PmbHLH66, PmbHLH84, PmbHLH85, PmbHLH94, and PmbHLH95) were highly expressed in root, while PmbHLH33 and PmbHLH58 had a high expression level in leaf and stem respectively, implying that the functions of genes in the same family gradually differ in the evolution process.

Previous studies have proved that the cold resistance of apricot mei is stronger than that of true mume. Therefore, annual plants of apricot mei ‘Danfenghou’ and true mume ‘Beijing Yudie’ were treated with low temperatures. Based on transcriptome data, this study analyzed differences in the expression level of the PmbHLHs gene family during these three periods. FPKM value greater than one is an effective expression. 62 genes that were effectively expressed after low-temperature stress were detected. The FPKM value of these genes was plotted as a heat map.

As shown in Figure 5, according to different expression patterns in different periods, PmbHLHs genes in the two varieties were clustered into four groups (Supplementary Table S8). PmbHLH genes of the ‘Danfenghou’ group Ⅰ were highly expressed during freezing treatments and had a similar expression pattern to the ‘Beijing Yudie’ group Ⅲ. There are six genes (PmbHLH64, PmbHLH25, PmbHLH10, PmbHLH78, PmbHLH31, and PmbHLH47) in the two varieties with similar expression patterns, which can respond to freezing stress. The expression level of PmbHLHs in the ‘Danfenghou’ group Ⅱ and ‘Beijing Yudie’ group Ⅰ gradually decreased with the extension of treatment time. The expression profiles of PmbHLHs in the ‘Danfenghou’ group Ⅲ and ‘Beijing Yudie’ group Ⅱ were similar. 14 genes increased expression levels when two varieties were treated at 4°C for 8 h. They might be involved in the perception and transport of cold signals. The expression levels of ‘Danfenghou’ group Ⅱ and ‘Beijing Yudie’ group Ⅰ both increased at 4°C for 5 days, but there were no overlapping genes in the two cultivars. This indicated that the function of PmbHLH genes in the two cultivars after prolonged low-temperature acclimation has a difference.

According to FPKM value and fold log2 change value, 10 PmbHLH genes were screened out. Compared with control material that has not been treated with low temperature, the transcriptional expression levels of six PmbHLHs genes (PmbHLH25, PmbHLH28, PmbHLH38, PmbHLH40, PmbHLH57, and PmbHLH78) increased by 2–5 times at chilling treatment (4°C for 8 h and 5 days), which may be involved in the response of P. mume to chilling treatment. The expression levels of PmbHLH4, PmbHLH6, PmbHLH26, and PmbHLH46 gradually decreased with the extension of 4°C treatment time. Compared with the control, the highest reduction factor reached 5 times. They might negatively regulate downstream low-temperature response genes or proteins to participate in the chilling stress response of P. mume. We used transcriptome to further analyze the expression patterns of PmbHLHs genes in different varieties. The trend of expression levels with these 10 differential genes showed the same in strong cold-resistant ‘Danfenghou’ and weaker cold-resistant ‘Beijing Yudie’. However, there were differences in multiple variations among different varieties. For example, the expression of the PmbHLH28 gene in ‘Beijing Yudie’ was 4.8 times higher than that in the control, but only 1.1 times higher in the ‘Danfenghou’ group.

To further investigate PmbHLHs function in low temperature comprehensively, the 10 genes were detected by qRT-PCR experiments in which the stems of P. mume were treated for chilling treatment (4°C for 8 h and 5 days) and freezing treatment (−5°C for 1 h). The expression of PmbHLH genes was distinct in two varieties of P. mume (Figure 6A). In true mume ‘Beijing Yudie’, six genes were up-regulated in varying degrees under low-temperature stress. Remaining genes were down-regulated. Among up-regulated genes, PmbHLH25, PmbHLH40, PmbHLH46, and PmbHLH57 were highly expressed under 4°C and −5°C treatment. The greatest expression of PmbHLH25 and PmbHLH40 was found in 4°C treatments only. In addition, low temperature induced most of the genes to change significantly in apricot mei ‘Danfenghou’. PmbHLH46 and PmbHLH25 were up-regulated at 4°C treatment and at −5°C treatment, respectively. PmbHLH40 and PmbHLH57 were highly expressed under 4°C and −5°C treatment. In two varieties, four genes (PmbHLH25, PmbHLH38, PmbHLH40, and PmbHLH78) could be induced to high expression in true mume ‘Beijing Yudie’ and apricot mei ‘Danfenghou’. In addition, the expression of four PmbHLHs (PmbHLH38, PmbHLH40, PmbHLH57, and PmbHLH78) in apricot mei ‘Danfenghou’ higher than in true mume ‘Beijing Yudie’, further showing that apricot mei ‘Danfenghou’ was more resistant than true mume ‘Beijing Yudie’.

To explore the regulation process of PmbHLHs, we tested the expression pattern of 10 genes by qRT-PCR at chilling treatment (0/2/4/8/16/24/48/72 h) and freezing treatment (0/2/4/6/8/10/24/48 h) for different periods using true mume ‘Beijing Yudie’. As shown in Figure 6B, the expression level of PmbHLH04, PmbHLH25, PmbHLH26, PmbHLH46, and PmbHLH57 were reduced with chilling treatment. However, PmbHLH38 and PmbHLH40 were induced to high expression and up-regulation expression peaked on treating with 16 h. Under freezing treatment, the expression level of PmbHLH26, PmbHLH40, and PmbHLH78 changed slightly. PmbHLH38 and PmbHLH40 presented a trend of decreasing. After prolonging the freezing treatment time, PmbHLH04, PmbHLH06, PmbHLH25, PmbHLH28, and PmbHLH57 were induced to high expression. Especially, PmbHLH25 had the largest expression level increase (approximately 10-fold) at 48 h after being exposed to freezing conditions (Figure 6C). These results imply that PmbHLH genes may play a role in resisting low-temperature stress.

The promoter cis-elements (1500 bp) were analyzed through PlantCARE software and the predicted regulation mechanisms of PmbHLHs. The result showed that PmbHLH genes were abundant in abiotic and biotic elements (LTR, MBS, WRE3, WUN-motif, ARE), light-responsive elements (AAAC-motif, ATCT-motif, LAMP-element, Box 4, I-box, and G-box), plant growth and development-related elements (MSA-like, GCN4_motif, and RY-element). In addition, hormone-responsive elements (ABRE, TGACG-motif, G-box, MYC, P-box, and TATC-box) exhibited a wide range of positions in the promoter. This means that PmbHLHs may extensively participate in various physiological biochemistry pathways of P. mume (Supplementary Figure S4).

To further predict the functions of PmbHLHs, we constructed the interaction network using AraNet V2 based on homologous proteins of Arabidopsis (Lee et al., 2015). A total of 56 PmbHLH proteins had orthologs in Arabidopsis and predicted about 640 interaction protein pairs (Figure 7 and Supplementary Table S9). The interaction network showed that PmbHLHs might interact with MYB, bHLH, NAC, bZIP, HB, WRKY, ERF, and so on, implying PmbHLHs might exert functions by interacting with other genes. The above results predicted that PmbHLHs might be involved in the response of P. mume to low temperature. To further research ten PmbHLHs, we constructed and analyzed the interaction network of candidate ten genes. Among the ten PmbHLHs, PmbHLH04, PmbHLH26, and PmbHLH38 were absent in homologous proteins of Arabidopsis. The three candidate genes might be novel and potential proteins in response to low temperature. The other seven candidate PmbHLHs formed an interactive network centered on PmbHLH40, which showed high homology to JAM2. Moreover, they might strongly interact with MYB124, TIFY, NAC, WRKY, and IAA to drive function when suffered from low-temperature stress in P. mume. Supplementary Figure S6 shows that PmbHLH25 (AKS2) could directly interact with PmbHLH28 (AT5G57150). PmbHLH06 (bHLH93) could directly interact with PmbHLH38 (FBH4). This indicated that these PmbHLH might work by forming dimers under low-temperature stress. Overall, interaction networks could provide a crucial reference for investigating the regulation mechanism of PmbHLHs.