Genome-wide identification of the bHLH transcription factor family in Rosa persica and response to low-temperature stress

Identification of members of the RbebHLH family

The genome sequence data, gene structure annotation information, and functional information for R. persica were obtained from our previous research (unpublished). The AtbHLH protein sequence was downloaded from the plant transcription factor database PlantTFDB (http://planttfdb.gao-lab.org). Using the TBtools Blast Wrapper tool (Chen et al., 2020), based on the bHLH gene family reference sequence of A. thaliana, possible sequences in R. persica were retrieved with an E-value threshold of <1 × 10⁻⁵. NCBI Batch CD Search (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) (Marchler-Bauer et al., 2015) was used for searches based on the conserved bHLH gene family structure domain. Detection of conserved structural domains in sequences by HMMER (https://www.ebi.ac.uk/Tools/hmmer/) and Pfam databases (http://pfam.xfam.org/). Genes that did not contain bHLH domains were excluded. TBtools Fasta Stats was used to evaluate the CDS length, the length of the protein, and chromosomal location information. Expasy PI (https://web.expasy.org/compute_) online tools were used to evaluate molecular weight (MW) and protein isoelectric (pI) points. The online tool WoLF PSORT (https://wolfpsort.hgc.jp) (Horton et al., 2007) was used for subcellular localization prediction. TBtools was used to visualize the chromosomal location and tandem repeat genes of RbebHLHs.

Multiple sequence alignment and gene structure analysis

MEGA (https://www.megasoftware.net/) and Jalview (https://www.jalview.org/) were used to generate a bHLH protein sequence alignment, evaluate known conserved sites, and explore new conserved sites. MEME was used to analyze the conserved motifs (Bailey & Elkan, 1994). The motif length was set to 6–200 amino acid residues, the maximum number of motifs was 10, and default values were used for other parameters. The NCBI CD-Search (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) online tools to predict the conserved motifs for functional annotation. Structural information was obtained using the R. persica genome annotation file, and the results were visualized using TBtools to draw the distribution of conserved motifs and a gene structure map.

Phylogenetic and gene duplication analysis

AtbHLH transcription factors from the A. thaliana TAIR database (https://www.arabidopsis.org) (Reiser et al., 2017) were downloaded, including data for proteins encoded by 162 AtbHLHs. The Muscle function in MEGA was used to generate a multiple sequence alignment of bHLH proteins of Arabidopsis and R. persica. ‘JTT+G’ was identified as the optimal amino acid substitution model for the phylogenetic analysis (Nei & Kumar, 2000). The FastTree function (Price, Dehal & Arkin, 2009) in Geneious Prime was used to construct the maximum likelihood evolutionary tree (1,000 bootstrap replicates), and iTOL (https://itol.embl.de) (Letunic & Bork, 2021) was used to display the phylogenetic tree. The chromosomal location of each gene was obtained from the GFF file of the R. persica genome. Additionally, genome data for A. thaliana, P. mume, P. persica, F. vesca, and R. chinensis were downloaded from the NCBI database. BLASTp (Mahram & Herbordt, 2015) and MCScanX (Wang et al., 2012) were used to analyze gene duplication events. Additionally, TBtools was used to calculate the non-synonymous substitution rate (K_a), synonymous substitution rate (K_s), and their ratio for homologous bHLH pairs in R. persica and R. chinensis to evaluate the type of selection pressure acting on protein-coding genes. Based on λ (i.e., the substitution rate per synonymous site per year), K_s values were converted to divergence times (T) in millions of years. Circos and collinearity diagrams were drawn using TBtools to visualize the gene density and the collinearity.

Cis-acting element analysis and transcriptional target gene network prediction

The GFF Sequence Extractor tool in TBtools was used to extract the 2,000 bp sequence upstream of the bHLH transcription factor start codon from the R. persica annotation file. PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) online tools were used to predict RbebHLH regulatory elements in gene promoter regions. The HeatMap in TBtools was used to draw a statistical map of the distribution of homeopathic elements. The JASPER database (https://jaspar.genereg.net/) was used to obtain transcription factor binding motif information, and MEME FIMO was used to predict target genes. Cytoscape was used to map the network of interactions between transcription factors and target genes.

Plant material and qRT-PCR analysis

Monocotyledonous rose material was harvested from Hutubi County, Xinjiang, China, next to National Highway 312. From mid-September to mid-April, the main roots (about 10 cm from the root base) and mature stems (about 15 cm from the stem base) were collected from Hutubi under natural overwintering conditions at intervals of 1 month, and more than 30 g of each part was sampled each time. The collected material was fixed in liquid nitrogen for 30 s, and then transported back to the laboratory on dry ice to maintain a low-temperature cold chain for storage in an ultra-low temperature refrigerator. Samples were used for RNA extraction, complementary DNA (cDNA) library construction, and transcriptome sequencing based on the R. persica reference genome. Raw data obtained by sequencing have not yet been published. The expression level was calculated and normalized according to the fragments per kilobase per million mapped reads (FPKM). The TBtools cartoon heat map was used to visualize the gene expression levels.

RNA was extracted using an RNA Extraction Kit (Omega Bio-tek, Norcross, GA, USA). cDNA was then synthesized using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Kusatsu, Japan). Quantitative real-time PCR(qRT-PCR) experiments were performed using the TB Green premix Ex Taq II Kit using the qTOWER2.2 system (Analytik Jena, Jena, Germany). Three biological and three technical replicates were performed for each sample. All primers are shown in Schedule 5. The reaction system consisted of 25 μL of TB Green Premix Ex Taq II (12.5 μL), forward and reverse primer mix (2 μL), cDNA sample (2 μL), and sterilized water (8.5 μL). The reaction program was as follows: initial denaturation at 95 °C for 30 s, followed by 40 repeated cycles at 95 °C for 5 s, 55 °C for 30 s, and 72 °C for 30 s. Relative gene expression levels were calculated using the 2^−∆∆Ct method (Zhang et al., 2021b) using the RbeGAPDH (glyceraldehyde-3-phosphate dehydrogenase) gene as an internal reference.

Physiological index measurement method

The relative conductivity was determined as follows: Wash the test material with tap water and absorb it, cut it into small pieces about 1 cm long, mix it well, put 1 g of plant in each test tube, and add 10 mL of deionized water. Place the tubes in a vacuum desiccator, pump for 10–15 min, deflate and shake lightly for 20 min, and measure the conductivity of the extracts with a conductivity meter. The conductivity of the deionized water blank was also measured. The conductivity was measured by placing the measured tubes in a boiling water bath at 100 °C for 15 min and then cooling down for 10 min before measuring the conductivity after boiling. Relative conductivity = (sample conductivity-blank conductivity)/(boiling conductivity-blank conductivity) × 100%. The malondialdehyde content was determined using the thiobarbituric acid colorimetric method (Hodges et al., 1999). Five biological replicates and three technical replicates were performed. Correlation analysis of physiological indicators and expression of selected genes was performed using Origin 2021.

Identification of members of the RbebHLH family

By a Blast search and verification against the NCBI CD database, 337 candidate bHLH sequences in the R. persica genome were obtained. Genes that did not contain the conserved bHLH domain were eliminated, resulting in 142 RbebHLHs, named RbebHLH1–RbebHLH142 (Table S1). Of these, 140 RbebHLHs were distributed on seven chromosomes. Two RbebHLHs failed to be assembled to the chromosomal level. Proteins were between 73 (RbebHLH33) and 960 (RbebHLH13) amino acids (aa), with a predicted MW of between 7643.82 (RbebHLH33) and 104532.56 (RbebHLH13). The predicted PI was between 4.42 (RbebHLH134) and 11.79 (RbebHLH111). The subcellular localization predicted by WoLF PSORT showed that 117 bHLHs were located in the nucleus, 12 in the chloroplast, six in the cytoplasm, three in the cytoplasm and nucleus, and one in the mitochondria, peroxisome, Golgi, and endoplasmic reticulum.

Multiple sequence alignment and gene structure analysis of RbebHLHs

The bHLH region is located at the C-terminus of the bHLH domain, which is mainly composed of a hydrophilic lipid helix composed of hydrophobic residues and a loop and contains two highly conserved sites (Leu23-Leu54). It mainly promotes the interaction between proteins and forms homodimers or heterodimers, which interact with the E-box element in the gene promoter region (Hernandez et al., 2007). The conserved domain of the RbebHLH protein obtained by simultaneous artificial construction using MUSCLE is shown in Fig. 1. The Basic domain contained four highly conserved amino acids (H5, E9, R10, R11, and R12), However, the bHLH region contained 18 relatively conserved amino acids (E14, R19, L23, R24, L26, V27, P28, K32, D34, A36, S37, L39, E41, A42, I43, Y45, K47, L49 and L57), among which R11, R12, L23, and L49 were highly conserved. These residues may play a crucial role in the binding of bHLHs to target genes and the formation of protein dimers.

The results of the RbebHLHs structural analysis are shown in Fig. 2. Motif prediction by MEME identified 10 conserved motifs. Motif1 and Motif2 were highly conserved in almost all RbebHLH proteins. In addition, each group had some degree of motif specificity. For example, Motif6 only existed in subfamily III and subfamily X, and the combination of “Motif9+Motif10+Motif5” was only detected in subfamily V. The specificity of the motifs may contribute to the functional differences of bHLHs in different groups. To clarify the evolution of the RbebHLHs family, the bHLH protein sequences of R. persica were evaluated by a phylogenetic analysis and exon-intron structure analysis. The RbebHLH proteins were divided into 10 subfamilies. The genes in each subfamily contained similar numbers of exons and introns with relatively conserved positions. The number of exons ranged from 1 (RbebHLH142, etc.) to 11 (RbebHLH123, etc.), and only about 10% of the genes contained eight or more exons. There were significant differences in the number of introns among different subfamilies. RbebHLH58 had the largest number of introns, and 18 genes did not contain introns. There were three genes containing only the 5′-UTR, five genes containing only the 3′-UTR, and 31 genes containing both the 5′-UTR and the 3′-UTR. These findings indicate that RbebHLH proteins in the same subfamily generally have a similar composition of conserved motifs, with some variation.

Phylogenetic and gene duplication analyses

We performed a phylogenetic analysis of bHLH protein sequences from R. persica (142 accessions) and A. thaliana (162 accessions). Based on the classification of Arabidopsis bHLH proteins (Heim et al., 2003), these 304 bHLHs proteins were classified into 26 subfamilies, named from Ia to XV (Fig. 3). We combined some of the similar subfamilies into 21 subfamilies and showed that subfamily XII has the highest number of RbebHLHs (24) in R. persica, and subfamilies II, IIIe, IX have the lowest number, all of which contain only one RbebHLH. It has been shown that proteins belonging to the same subfamily mostly have similar functions, and therefore the phylogenetic results contribute to the function of RbebHLHs proteins.

Chromosome mapping revealed that chromosome 6 had the largest number of RbebHLHs (26), chromosome 3 had the least (11), and chromosome 5 had the largest number of RbebHLHs. There were seven pairs of tandem duplications in RbebHLHs (RbebHLH15 and RbebHLH16, RbebHLH26 and RbebHLH27, RbebHLH52 and RbebHLH53, RbebHLH54 and RbebHLH55, RbebHLH65 and RbebHLH66, RbebHLH111 and RbebHLH112, and RbebHLH127 and RbebHLH128) (Fig. S1) and 118 pairs of segmental duplications (Fig. 4). Tandem duplications were present in three pairs on chromosome 6 and in one pair each on chromosome 1, 2, 4, and 5. The most segmental duplications were found on chromosome 2 (32), followed by chromosome 7 (20) (Table S2).

To further explore the evolution of the bHLH gene family, we performed a comparative genomic analysis and collinearity analysis among A. thaliana, Prunus. mume, Fragaria vesca, Rosa chinensis, and Prunus persica and R. persica. There were 127, 150, 151, 172, and 164 pairs of segmental duplications between R. persica and A. thaliana, P. mume, F. vesca, R. chinensis, and P. persica, respectively (Fig. 5). Substantial collinearity was detected between genes in R. persica and the other three species in Rosaceae, indicating that they were closely related (Table S3).

To determine the type of selection acting on the collinear genes in R. persica and R. chinensis, the ratio of the nonsynonymous substitution rate (K_a) to the synonymous substitution rate (K_s) was calculated (Table S4). In R. persica and R. chinensis, only the K_a/K_s ratio of RbebHLH111 and XM_024337436.2 was greater than 1.0 (2.88), suggesting that the gene was under positive selection. The remaining gene pairs had K_a/K_s ratios of less than 1.0, consistent with purifying selection. According to the K_s values, the divergence time of RbebHLHs and RcbHLHs was calculated (Table S4). Divergence time estimates ranged from 1.03 to 274.77 Mya with two distinct peaks, i.e., 115 gene pairs diverged less than 10 Mya and the remaining 44 genes diverged more than 80 Mya.

Cis-acting element analysis and transcriptional target gene network prediction

To further investigate the potential mechanism by which RbebHLHs contribute to the response to abiotic stress, cis-regulatory elements were predicted in the promoter regions. According to functional prediction results, cis-regulatory elements were divided into three categories: environmental stress elements, hormone response elements and development-related elements. The number and type of homeotropic elements for each gene were summarized in a heat map (Fig. 6). RbebHLH23 and RbebHLH87 contained the most cis-acting elements (37), while RbebHLH28 contained the least number of cis-acting elements (two). It can be found that most of the RbebHLHs contain both STRE and ARE elements in significantly greater numbers than other cis-acting elements. RbebHLH76 and RbebHLH87 contained the largest number of environmental stress-related elements (15), and RbebHLH28 and RbebHLH63 contained the fewest (one). In addition, we found that most of the RbebHLHs promoter regions have CGTCA-motif and TGACG-motif elements associated with the methyl jasmonate response, suggesting that this part of RbebHLHs may respond to low-temperature stress through the jasmonate signaling pathway.These results suggest that RbebHLHs are widely involved in various physiological and biochemical activities of plants and in response to environmental stimuli and stresses.

To understand the regulatory relationship between RbebHLHs and the C-repeat binding factors (CBF) family of key cold resistance genes, we used the JASPER database (http://jaspar.genereg.net/) to predict motif information and used MEME FIMO to predict target CBF genes (Fig. 7). A total of 21 RbebHLHs were predicted to interact with five CBFs. The 22 RbebHLHs included phytochrome-interacting factor1 (PIF1), PIF3, PIF4, 3 SPTs, 3 MYC2, 2 BIM2, DYT1, UNE10, etc. PIF3, PIF4, 3 SPTs, and two BIM2 interacted with four of the CBFs, and the maximum number of RbebHLHs that interacted with both CBF1 and CBF4 was 15. These results suggest that a close relationship exists between the RbebHLHs and CBF families and may regulate the plant response to low-temperature stress through interactions.

Expression characteristics of selected RbebHLHs during overwintering

To understand the response of RbebHLHs to low-temperature stress, we selected 30 genes with a high number of cis-acting elements in the promoter region and 21 genes with reciprocal relationships with the CBF family and analyzed their expression patterns during the overwintering period (September-April) (Fig. S3). A total of 45 RbebHLHs were analyzed for expression in roots and stems during overwintering, excluding duplicate genes. More RbebHLHs were highly expressed in roots in April and in stems in September, while the expression of RbebHLHs was generally low in roots and stems in December, January, and February. Both RbebHLH69 and RbebHLH86 were not expressed in roots and stems during overwintering, which may be related to the specificity of tissue expression. RbebHLH8, RbebHLH77, RbebHLH83, RbebHLH121 and RbebHLH124 showed high expression in the roots in April, whereas they were largely absent in the other months. The expression of RbebHLH12, RbebHLH20, RbebHLH23 and RbebHLH87 gradually decreased during overwintering, while the expression of most of the selected RbebHLHs showed a decreasing and then increasing trend during overwintering. These expression properties suggest that the functions of RbebHLHs have diverged to different degrees during evolution.

Given the results of the promoter analysis, we also focused on the expression of genes in the jasmonate signaling pathway with the RbebHLHs gene. Expression heatmap analysis revealed that JAZ6, CBF5, RbebHLH20, RbebHLH87, RbebHLH90 and RbebHLH130 had similar expression characteristics during overwintering, all with high expression in September. JAZ8, JAZ10, CBF1, RbebHLH78 and RbebHLH131 were highly expressed in April (Fig. 8). To verify the accuracy of the transcriptome data, we selected 4 months in the root for qRT-PCR analysis (Table S4). The qRT-PCR results were in general agreement with the transcriptome data, suggesting that RbebHLHs may regulate the response of R. persica to low-temperature stress through the jasmonic acid signaling pathway.

Expression of RbebHLHs correlates with plant cold tolerance

To investigate the relationship between RbebHLHs expression and the cold tolerance of R. persica, we measured the relative conductivity and MDA of R. persica roots and stems during overwintering (Fig. 9). Relative conductivity and malondialdehyde content are indicators of cell membrane integrity and lipid peroxidation, respectively. The conductivity of both roots and stems showed a trend of increasing and then decreasing, with the lowest conductivity in September and a peak in January for roots. Stem conductivity was lowest in October and peaked in February. The malondialdehyde content also showed a trend of increasing and then decreasing, with the lowest malondialdehyde content in the root part in April and the peak in February. The relative conductivity was higher in December, January and February than in other months, while malondialdehyde content was higher in January, February and March. This indicates that R. persica was more severely affected by low-temperature stress in these months. Correlation analysis of the expression of the selected 45 RbebHLHs with relative conductivity and malondialdehyde content showed that six genes were significantly correlated with malondialdehyde content and seven genes were significantly correlated with relative conductivity (Fig. 9C). Among them, RbebHLH20, RbebHLH40, RbebHLH58 and RbebHLH113 were significantly negatively correlated with both malondialdehyde and relative conductivity, while RbebHLH130 and RbebHLH115 were significantly negatively correlated. These results suggest that RbebHLHs play an important role in regulating R. persica in response to low-temperature stress.