In total, 104 sequences with conserved bHLH domains were identified as members of the P. sibirica bHLH family. The physical and chemical properties of the PsbHLHs are shown in Table 1. The lengths of the PsbHLHs range from 91–1159 amino acids (aa), and most of these proteins (74%) had lengths of 200–500 aa. The total number of atoms ranged from 1444 (PsbHLH16) to 17840 (PsbHLH18). PsbHLH18 had the largest molecular weight (Mw) of 128.17 kDa, while in other PsbHLHs it ranged from 10.27 kDa (PsbHLH16) to 85.93 kDa (PsbHLH101). The isoelectric point (pI) ranged from 4.6 (PsbHLH22) to 10.31 (PsbHLH32), with the average pI being 6.77. Of these, 66 PsbHLHs had pI <7, indicating acidity, while 38 PsbHLHs had pI >7, indicating alkalinity. The aliphatic index (AI) ranged from 52.34 (PsbHLH50) to 98.96 (PsbHLH80). The instability index (II) ranged from 21.81in PsbHLH16 to 82.05 in PsbHLH5. A total of 101 PsbHLHs greater than 40 were stable proteins, while three were less than 40 and belonged to unstable proteins. The predicted grand average of hydropathicity (GRAVY) values were all negative and ranged between −0.972 (PsbHLH88) and −0.271 (PsbHLH83), representing that all proteins have hydrophilic properties. Subcellular localization assays revealed that 102 PsbHLHs (98.08%) were located in the nucleus; however, PsbHLH47 was located in the chloroplast, and PsbHLH71 was found in the mitochondria.

Multiple sequence alignments were performed to analyze the sequence features and the presence of conserved amino acid residues in PsbHLHs. The bHLH structural domain of P. sibirica consisted of 66 amino acids. Visualizing the results of multiple sequence alignment allowed the identification of four conserved regions, one basic region (ERRR), two helical regions (IRLLLVP or ASLEAIYKLQQL), and a loop region (KD) (Figure 1). The size of the amino acid letters corresponding to each point represents the conservation of amino acids at each point in the bHLH transcription factor gene domain. The conservation of amino acid residues at 19 sites is greater than 50%. The conservation of Arg-11, Arg-12, Pro-27, Leu-22 and Leu-59 residues was greater than or equal to 90% implying the site was highly conserved. In addition, Helix1 is more conserved than Helix2.

Using the protein sequences of the 104 PsbHLHs and 135 A. thaliana bHLHs, a neighbor-joining phylogenetic tree was created to examine the evolutionary relationships and classification of PsbHLHs members. All of the bHLH proteins were divided into 26 subfamilies based on how they were categorized in A. thaliana’s bHLH proteins. The nomenclature of the PsbHLH subfamilies was consistent with that of the A. thaliana bHLH group. Among them, there were multiple members of AtbHLHs, in the Orphans, XIV, and XIII groups, but there were none of the PsbHLHs. As shown in Figure 2, the Ib subfamily had the largest number of family members, containing 12 PsbHLHs; and the smallest were subfamilies VI and XV, each containing only one PsbHLH. According to the statistics of the number of AtbHLH and PsbHLH members in each subfamily, the number of bHLHs in Va, XV, VIIIa, VIIIb, VIIIc,VI, IVb, IVd, III(d+e), IIIb, and III(a+c) were the same.

According to the genome annotation information, the 104 PsbHLHs were randomly distributed on eight chromosomes and numbered PsbHLH1–PsbHLH104 based on their chromosomal locations, as shown in Figure 3. The number of PsbHLHs on the chromosomes ranged from seven to 22. Chromosome 1 contained the highest number of PsbHLHs (22). The lowest number of PsbHLHs was present on chromosome 3, which had only seven PsbHLHs. Chromosomes 2, 7, and 8 contained 13 PsbHLHs each, whereas chromosomes 4, 5, and 6 contained 9, 15, and 12 PsbHLHs, respectively (Figure 3). In addition, except for chromosome 3, all the chromosomes exhibited member aggregation.

Gene duplication is considered the primary driving force in the evolution of genomes. Analysis of the segmental and tandem replication of the gene family plays a significant role in explaining the process of gene family expansion. In this study, we investigated the duplication events in the PsbHLH gene family. Gene duplication events were observed on all the chromosomes. Genes with segmental duplications were most frequently found on chromosome 1 (11 genes), followed by chromosome 5 (8 genes). In total, 48 gene duplication events representing approximately 46% (48 of 104) of the total PsbHLHs were found. A total of 34 genes underwent segmental duplication, which accounted for 65% of all syntenic relationships and played an essential role in the evolution of the PsbHLHs. These segmentally duplicated PsbHLHs formed 19 gene clusters that were linked to each other by colored lines, as presented in Figure 4. Among them, PsbHLH10, PsbHLH41, and PsbHLH65 are duplicate gene pairs with each other. PsbHLH18 forms two duplicate genes with PsbHLH86 and PsbHLH93. Moreover, 16 (35%) PsbHLHs, including PsbHLH16/17, PsbHLH32/33, PsbHLH33/34, PsbHLH55/56, PsbHLH56/57, PsbHLH57/58, PsbHLH58/59, PsbHLH73/74, PsbHLH80/81, and PsbHLH97/98, were demonstrated to be tandem duplicated genes. Gene family members formed by tandem replication are closely arranged on the same chromosome and form gene clusters with related sequences and functions. Each gene pair with duplication events belonged to the same subfamily, except for PsbHLH18 and PsbHLH86 segmental duplication gene pairs belonging to subfamilies Ib and XIII(a+b), respectively.

The selection pressure on a gene throughout evolution can be determined by Ka/Ks ratio. The Ka/Ks values for the PsbHLH gene family were calculated (Table 2). Ka/Ks >1 indicates positive selection, Ka/Ks<1 indicates purified selection, and Ka/Ks=1 indicates neutral selection. All duplicated PsbHLH pairs with Ka/Ks < 1 might have undergone purifying selection, suggesting that the majority PsbHLHs evolved slowly. Therefore, the Ks value was used to calculate the divergence time of duplicated gene pairs (Table 2). The divergence time of these gene pairs ranged from 12.22 to 78.54 Mya. The earliest gene duplication event occurred at 78.54 Mya.

To learn more about the evolutionary process of the PsbHLH gene family, a syntenic map of the bHLH genes of four Rosaceae species (Prunus salicina, P. mume, P. persica, and P. avium) and the model plant A. thalian was constructed. Red lines represent homologous gene pairs in P. sibirica and other plant genomes. Between P. sibirica and P. avium, 93 orthologous gene pairs were identified, followed by P. persica, P. mume, P. salicina, and A. thalian, with 89, 88, 86, and 66 homologous gene pairs, respectively. P. sibirica had more homologous genes with the four plants in the Rosaceae family, which is indicative of the taxonomic relationships between the species (Figure 5).

The prediction of gene structure plays an important role in understanding the evolution of gene family members. The exon-intron structures within PsbHLHs were investigated to determine their genetic features and structural characteristics. Different subfamilies had different intron/exon patterns, but the PsbHLHs sequences within the same subfamily had similar numbers of exons and introns, especially for PsbHLHs on homologous branches (Figure 6). The number of introns ranged from 0 to 17. No introns were present in PsbHLH19, PsbHLH40, PsbHLH69, PsbHLH75, or PsbHLH89. Exon analysis revealed that the number of exons ranged from 1 to 16. In addition, 34 PsbHLHs had no untranslated region (UTR). Gene structure analysis revealed differences in the structural composition of PsbHLHs.

Twenty motifs were identified in the PsbHLHs using Motif Elicitation (MEME). As shown in Figure 6, different PsbHLHs featured various amounts of conserved motifs, ranging from 1 to 12. Motif 1 and motif 2 constitute the conserved domain of bHLH, motif 1 contains one basic and one helical region, whereas motif2 contains one loop and another helical region (Table S1). Except for PsbHLH2, PsbHLH3, PsbHLH16, PsbHLH27, and PsbHLH103, almost all of the PsbHLHs contained motif1 and motif 2, and both of these conserved domains were close to one another. Interestingly, PsbHLH18 has two conserved domains (motif 1 and motif 2). PsbHLHs grouped in the same subfamily frequently shared consistent motif patterns. In addition, motifs 19, 12, 13, and 17 occurred only in subfamilies XII, VIIIb, III (d+e), and Ib, respectively.

Cis-regulatory elements play an important role in regulating plant growth, development, and stress responses. To predict the regulatory pattern of PsbHLHs, the −2000 bp sequence upstream of ATG were analyzed in the promoter region of 104 PsbHLHs using PlantCARE (Figure 7). The discovery of numerous hormone- and stress-responsive elements suggests that PsbHLHs may possibly be crucial in hormone signal transduction and stress response. Cis-elements we screened were roughly divided into two categories: stress-responsive elements (TC-rich repeats and LTR) and hormone-responsive elements (ABRE, TGA-element, TCA-element, TGACG-motif, CGTCA-motif, GARE-motif, TATC-box, and P-box). Most cis-elements in the promoter regions of PsbHLHs were ABRE (289), followed by the TGACG motif (133) and the CGTCA motif (133). Additionally, 39 members had auxin-responsive element (TGA-element), 47 members had cis-acting element involved in defense and stress responsiveness (TC-rich repeats), 69 members had cis-acting element involved in gibberellin-responsiveness and gibberellin-responsive element (GARE-motif, TATC-boxes and P-box), 52 members had cis-acting element involved in low-temperature responsiveness (LTR), 55 members had cis-acting element that was responsive to salicylic acid (TCA-element), 89 members had cis-acting element that was responsive to the abscisic acid (ABRE), and 67 members had cis-acting regulatory element that was responsive to the MeJA (TGACG-motif and CGTCA-motif). In addition, PsbHLH9 had the most cis-elements (22), followed by PsbHLH86 (20).

To better understand the function of PsbHLHs as TFs, the EggNog website was used for gene ontology (GO) annotation analysis. Among the 104 PsbHLHs, 72 were annotated using GO (Figure 8). The results showed that 72 PsbHLHs were annotated to eight cellular components, 70.83% of the genes were located in cells and cell parts, 69.44% in organelles and less than 20% in other cellular components. Molecular functions include transcriptional regulator activity, binding, molecular transducer activity, and catalytic activity. PsbHLHs were annotated in 17 biological processes. All PsbHLHs were involved in the regulation of biological processes, and 71 genes (98.6%) were involved in cellular and metabolic processes. Thirty genes were involved in response to stimuli, 22 genes in the developmental processes, and 15 genes in reproductive processes. In addition, PsbHLH16 had the most annotations (16), whereas PsbHLH91 was annotated for 11 biological processes.

Twenty candidate genes were identified by integration and overall consideration of the results from the phylogenetic analysis, GO annotation, and cis-element analysis of PsbHLHs to further study the expression patterns and regulation processes in P. sibirica. These 20 genes expression levels were quantified indifferent tissues (roots, stems, leaves, pistils, and stamens) using the expression levels in petals as reference. Nine genes in pistils and 11 genes in the roots were highly expressed. Four genes were highly expressed in the stems and stamens. While, only two genes (PsbHLH54 and PsbHLH88) were highly expressed in leaves. PsbHLH39 was highly expressed in all five tissues. Five PsbHLHs (PsbHLHs19/29/38/44/66) were expressed in one tissue (Figure 9). PsbHLH39 was highly expressed (>60 times) in all tissues except the stamens. The expression levels of seven PsbHLHs were upregulated by more than 60 times in the roots, PsbHLH5/29/60 exceeded by 80 times; PsbHLH11 was upregulated by 80 times in the pistil; PsbHLH35 was upregulated by more than 60 times in the pistil and stem; and PsbHLH73 was upregulated by 90 times in the stamen.

In order to further explore the involvement of PsbHLHs in low-temperature tolerance, the expression levels of 20 genes at low temperature (−4°C; 0 h, 15 min, 30 min, 1 h, and 2 h) were detected by qRT-PCR. The relative expression level (0 h) in the control pistils without treatment was normalized to 1. The expression patterns of 11 PsbHLHs (PsbHLH29/35/38/39/40/42/44/46/57/60/89) showed upregulation, whereas nine PsbHLHs (PsbHLH5, 11, 19, 24, 66, 67, 73, 88, and 91) displayed a downregulation. In the upregulated group, all 11 genes were highly expressed, particularly PsbHLH44, which displayed the highest expression level under the low-temperature treatment. Seven PsbHLHs (PsbHLH29/35/38/39/40/42/89) genes were first upregulated at 15 min, followed by slight downregulation at 30 min, and then continued to be upregulated at 1 or 2 h of treatment. In addition, the expression levels of PsbHLH44/46/57/60 reached a maximum at 15 min or 1 h and then decreased with prolonged treatment time (Figure 10).

Based on these syntenic relationships, orthologous AtbHLHs from 104 PsbHLHs were selected in A. thaliana to predict interactions. Most of the PsbHLHs were predicted to exhibit interactions (Figure 11). Sixty bHLHs interacted with more than one bHLH protein, with PYE (PsbHLH63) having the largest number of interacting partners (24). Detailed information on the interactions is summarized in Table 3. In addition, the gene function of the homolog PsbHLHs was detected based on its verified gene function in A. thaliana. Since PsbHLH42 is homologous to the cold-inducible gene ICE1, it is very likely to take part in P. sibirica’s reaction to low-temperature stress.

Based on bioinformatics studies, expression pattern analysis, and homology with ICE1, the involvement of PsbHLH42 in the cold response was hypothesized. To investigate the function of PsbHLH42 in low-temperature stress, the pRI101-PsbHLH42-GFP vector was constructed and transiently transformed tissue culture seedlings of P. ussuriensis by the Agrobacterium-mediated method. An empty vector was transiently transformed into P. ussuriensis as the control. The qRT-PCR results demonstrated that the expression level of PsbHLH42 in the transgenic plants was significantly higher than that in the control plants after transient transformation, indicating that PsbHLH42 was successfully transferred into P. ussuriensis and further analysis could be carried out. The seedlings of the control and transgenic plants were subjected to 4°C stress for 6 h to evaluate the cold resistance of PsbHLH42. Before low-temperature treatment, the physiological signs did not alter significantly between the transgenic and control strains (Figure 12). After low-temperature treatment, compared with the control plants, the transgenic plants that overexpressed PsbHLH42 had higher soluble sugar and soluble protein contents, higher Superoxide dismutase (SOD) and peroxidase (POD) activities, and lower electrical conductivity and Malondialdehyde (MDA) content.

Sequences of P. sibirica (Prunus sibirica F106 Genome v1.0) were obtained from the Genome Database for Rosaceae (GDR, https://www.rosaceae.org/) (Jung et al., 2019). The hidden markov model (HMM) file of the bHLH domain (PF00010) was downloaded from the Pfam database (http://pfam.xfam.org/) (Mistry et al., 2021) and used for protein screening using HMMER software (version 3.3.2). The online CD-search tool (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (Lu et al., 2020) and Pfam were used to validate the bHLH domains of these protein sequences, and the physical and chemical parameters were calculated using ExPASy (https://www.expasy.org/) (Duvaud et al., 2021). The CELLO program (http://cello.life.nctu.edu.tw/) (Yu et al., 2006) was used to predict the subcellular localization.

MEGA software (version 11.0) and Jalview software (version 2.11.2.0) were used to visualize the sequences of the conserved domains in PsbHLHs. The sequence logo for the bHLH domain was created by submitting multiple alignment sequences to WebLogo (http://weblogo.berkeley.edu/logo.cgi) (Crooks et al., 2004).

The full-length bHLHs protein sequences from P. sibirica and A. thaliana were aligned using ClustalX (version 2.1). Thereafter, a phylogenetic tree was constructed based on the neighbor-joining (NJ) method using MEGA with the following parameter settings: Poisson’s model, pairwise deletion option, and 1000 bootstrap replicates provided for statistical reliability. A phylogenetic tree was constructed using iTOL (https://itol.embl.de/) (Letunic and Bork, 2021). The A. thaliana bHLH protein sequences were downloaded from TAIR (https://www.arabidopsis.org/) (Berardini et al., 2015). Subfamily grouping of PsbHLHs was performed according to the AtbHLH classification scheme.

Chromosomal locations were identified and plotted using TBtools and the P. sibirica genomic database. The nonsynonymous (Ka) and synonymous (Ks) substitution ratios were calculated using the Simple Ka/Ks Calculator in TBtools software to acquire the natural purification selection between target gene pairs. The divergence time was calculated using the formula T = Ks/2r, with Ks being the synonymous substitutions per site, and the r of dicotyledonous plants being 1.5×10⁻⁸ synonymous substitutions per site per year (Gaut and Doebley, 1997). Finally, sequences with complete bHLH domains were preserved and numbered PsbHLH1 through PsbHLH104 to correspond with their locations extending from top to bottom on the chromosome.

The GFF and genome data of A. thaliana, P. persica, P. mume, P. persica, and P. avium were derived from the GDR. Syntenic analysis maps of P. sibirica and other representative plants were constructed using TBtools.

The structure of PsbHLHs was analyzed using TBtools. The conserved motifs of PsbHLHs were identified using the online tool MEME (http://meme-suite.org) (Bailey et al., 2015). The sequence of each PsbHLHs was retrieved by TBtools according to the genomic full-length DNA sequences of PsbHLHs.

Cis-elements in the promoter (−2000 bp upstream of ATG) were predicted using PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). Graphs of the structures were drawn using TBtools software.

Eggnog (http://eggnog5.embl.de/) (Huerta-Cepas et al., 2019) was used for GO annotation, and the results were visualized using WeGo (https://wego.genomics.cn/) (Ye et al., 2018).

All PsbHLH sequences were submitted to STRING (version 11.0, http://string-db.org) (Szklarczyk et al., 2019), and A. thaliana was chosen as the reference organism. After the BLAST analysis, the orthologous genes of A. thaliana with the highest scores were used to construct a network.

P. sibirica clone No. 453 (Appraisal of Improved Cultivar of Forest Trees, S-SV-PS-002-2021, Liaoning Provincial Department of Forest and Grassland) from the National Forest Germplasm Resource Preservation Repository for Prunus species at Shenyang Agricultural University (Kazuo, Liaoning, China) was used as the experimental material. The clone No. 453 sprays were placed in an artificial climate box during their peak flowering period and subjected to stress. Samples were collected at 0 h (Control), 15 min, 30 min, 1 h, and 2 h for −4°C low temperature stress. Six different tissue samples (roots, stems, leaves, stamens, pistils, and petals) in a natural environment were collected. All samples were set up with three biological replicates, immediately frozen in liquid nitrogen, and stored at −80°C until required for further analysis.

Total RNA was extracted from the samples using RNA Extraction Kit (Takara, Beijing, China) following the manufacturer’s instructions. The extracted RNA was then reverse transcribed to cDNA using the PrimeScript™ RT Master Mix Perfect Real Time (Takara, Beijing, China), and stored at −80°C. Primer software (version 5.0) was used to design specific primers (Table S2), which were synthesized by SYNBIO Technologies (Suzhou, China). The 18S rRNA gene of P. sibirica was used as the internal reference (Kuchipudi et al., 2012). The qRT-PCR was performed using PerfectStart® Green qPCR SuperMix (TransGen Biotech, Beijing, China) in StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, United States) following the manufacturer’s instructions. The qRT-PCR conditions were as follows: 94°C for 30 s, followed by 45 cycles of 94°C for 5 s and 60°C for 15 s, and 72°C for 10 s. The 2^−ΔΔCT method was employed to calculate the relative expression levels of PsbHLHs (Livak and Schmittgen, 2001).

The full-length coding sequence of PsbHLH42 was cloned into a pRI101-GFP vector. The construct was transformed into Agrobacterium tumefaciens GV3101, which was used to transform the tissue culture seedlings of P. ussuriensis (Liu et al., 2021).

The physiological characteristics of pRI101-GFP empty vector and pRI101-PsbHLH42-GFP seedlings (leaves) were determined after 6 h of 4°C stress treatment. SOD and POD activities were quantified using nitrogen blue tetrazolium (NBT), and guaiacol, respectively. MDA levels were determined using the thiobarbituric acid-reactive substances method. The soluble sugar content (SS) was determined using anthrone colorimetry. Soluble protein content (SP) was determined using the Coomassie Brilliant Blue G-250 staining method. Cell membrane permeability (CMP) was determined using a conductivity meter. The CMP was first determined using a conductivity meter after mixing fresh leaf tissue (0.3 g) with distilled water (15 mL) had been mixed for 4 h. The test tubes were then placed in a boiling water bath for 30 min. After standing for 4 h, measurements were performed again using the same conductivity meter (Wang et al., 2021).