Genome-wide identification of bHLH gene family and its response to cadmium stress in Populus × canescens

Identification of the basic/helix-loop-helix family genes in P. canescens

The genomic sequence data of P. canescens were sourced from the Aspen database (https://www.aspendb.org/downloads and the annotation version was sPta717alba_v2). bHLH transcription factors of Arabidopsis thaliana (AtbHLHs) were obtained from Carretero-Paulet et al. (2010), and those of Populus trichocarpa (PtbHLHs) were retrieved from NCBI’s database. Hidden Markov model (HMM) files for the bHLH domain (PF0010) were acquired from interProScan (https://www.ebi.ac.uk/interpro/entry/pfam/PF00010/), which were then utilized to identify bHLH proteins in P. canescens using the SPDE software (Xu et al., 2021). BLASTp searches were also performing using AtbHLHs against the P. canescens amino acid sequence data. To confirm the presence of the bHLH domain (E-value <1e⁻⁵), CD-search (https://www.ncbi.nlm.nih.gov/) and SMART (http://smart.embl-heidelberg.de/) were employed to analyze the identified PcbHLH sequences. The molecular weights (kDa) and isoelectric points (pI) of the PcbHLHs were calculated using the SPDE software.

Chromosomal locations of PcbHLHs

Chromosomal locations of the PcbHLH genes were determined using the P. canescens database (https://www.aspendb.org/downloads), and a distribution map was created with TBtools (Chen et al., 2023).

Multiple sequence alignments of PcbHLHs

Multiple sequence alignments were conducted using ClustalX in MEGA7, followed by visualization with Jalview (version 1.8.3). The variable sequences at the N- and C-terminal regions were excised, preserving the conserved domains in the central region. Sequence logos for the bHLHs were created by submitting the multiple alignment sequences to a specialized online platform (https://weblogo.berkeley.edu/logo.cgi). Employing the criteria delineated by Massari & Murre (2000), we classified the PcbHLH proteins into two principal categories based on the sequence information within the N-terminal region of the bHLH domains: DNA binding and Non-DNA binding (containing fewer than four basic amino acids). The DNA binding bHLHs were further categorized into two groups: E-box binding and Non-E-box binding (based on the presence of only Glu12 or Arg15). Consequently, E-box binding were further subdivided into two subgroups: G-box binding (characterized by His/Lys8, Glu12 and Arg16) and non-G-box binding (defined by the presence of only Glu12 and Arg16).

Phylogenetic analysis of PcbHLHs

To elucidate the evolutionary relationships between P. canescens and Arabidopsis bHLHs, we performed multiple sequences alignment of 170 complete PcbHLH sequences and 167 complete AtbHLH sequences using MEGA7. Subsequently, a Maximum Likelihood phylogenetic tree was constructed in MEGA7, with a bootstrap value set to 1,000 replicates for reliability. The phylogenetic tree was visualized using Itol (https://itol.embl.de/), which provided a clear depiction of the bHLHs’ evolutionary history.

Gene structure and protein motifs analysis

The exon/intron organization and splicing phase of the predicted PcbHLHs were analyzed using the Gff3 annotation files from P. canescens genome. This data were then graphically represented using TBtools. To identify conserved motifs within the PcbHLH proteins, we employed MEME (https://meme-suite.org/meme/tools/meme, version 5.5.4) to identify a maximum of ten motifs with an optimal width ranging from 10 to 100 amino acids. The resulting phylogenetic trees, gene structures, and conserved motifs were integrated and visualized within TBtools for comprehensive analysis.

Cis-acting element analysis in PcbHLH promoters

Promoter regions of PcbHLH genes, defined as 2000bp upstream sequences, were analyzed for cis-acting elements using PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/). TBtools facilitated the visualization of these elements, and a heat map was generated using Rstudio (version 4.2.2; RStudio Team, 2024) to represents the distribution of various cis-acting elements, providing insight into the regulatory mechanism of these genes.

Gene duplication and synteny analysis of PcbHLH genes

Gene duplication and synteny analysis provided insights into the evolutionary dynamics of bHLH genes in poplar and other plants. This analysis was performed using the Multiple Collinearity Scan toolkit (MCScanX), with results visualized in TBtools. Additionally, the nonsynonymous (ka) and synonymous (ks) substitution rates of bHLH gene pairs were calculated, offering further evidence of the evolutionary pressures acting on these genes.

Plant growth conditions and and Cd treatment

P. canescens seedlings, derived from micropropagation (Li et al., 2023), were cultured hydroponically at 25 °C under a photoperiod of 16 h light and 8 h dark within an artificial climate chamber. After acclimatization in half-strength Hoagland nutrient medium for one month, seedlings were treated with 10 µM CdCl₂ to induce Cd stress for 168 h. A concentration of 10 µM CdCl₂ was sufficient to elicit a transcription factor response in P. cansecens without causing harm. Samples of root, stems and leaves from individual P. canescens plants were collected at 0, 6, 12, 24, 48, 96, and 168 h post-treatment, with all samples collected simultaneously to minimize the impact of diurnal rhythms. All samples were collected in triplicate to negate the impact of biological variability.

RNA extraction and qRT-PCR analysis

Total RNA was extracted from root, stem, and leaf samples of P. canescens seedlings using the RNAprep Pure Plant Plus Kit (TIANGEN, Beijing, China). First-strand cDNA synthesis was performed using the PrimeScript™ RT Master Mix (TaKaRa, Dalian, China). Quantitative primers were designed using TBtools and qPCR was conducted with the 2 × Q3 SYBR qPCR Master Mix (TOLOBIO, Shanghai, China) on a 7300 Real-Time PCR System (Applied Biosystems, CA, United States). Technical replicates were conducted three times for each sample. Relative expression levels were quantified using the 2^−ΔΔCT method with EF1B (accession number FJ372570) gene serving as a reference (Wildhagen et al., 2010).

Identification and classification of bHLH protein in P. canescens

The initial bioinformatic screening of P. canescens using HMMER with hidden Markov models and BLASTp yielded a total of 179 predicted bHLH proteins. Subsequent filtering, employing CD-Search and SMART to exclude proteins with incomplete domains, culminated in the identification of 170 putative bHLH proteins (Table S1). These proteins were designated as PcbHLH1 to PcbHLH170, based on their chromosome locations. An analysis of their biochemical properties (Table S2) indicated a range of amino acid lengths from 90 for PcbHLH43 and PcbHLH154 to 740 for PcbHLH46. The molecular weights spanned from 10.29 kDa to 71.95 kDa, and the predicted isoelectric points ranged from 4.63 to 10.18. All proteins exhibited a grand average of hydropathy values below zero, indicative of their hydrophilic nature.

Chromosomal localization of PcbHLH genes

Distribution mapping revealed that 169 of the 170 PcbHLH genes are unevenly distributed across the P. canescens genome, with Chr2 harboring the highest number (17 genes). In contrast, only four PcbHLH genes were mapped to Chr3 and Chr7 (Fig. 1). Notably, PcbHLH170 remained unassigned to any chromosomes, which may be attributed to incomplete genome assembly.

Multiple sequence alignments of PcbHLHs

A comprehensive multiple sequence alignment of the 170 PcbHLH protein sequences was performed, revealing that the bHLH domain encompasses four conserved regions: a basic region, two helical regions, and a loop region. The basic region is characterized by specific residues, including His-5, Glu-9, Arg-10, Arg-12, and Arg-13, while the first helix region is defined by Ile-16, Leu-23, Leu-27, Val-28, and Pro-29. The loop region is predominantly composed of Lys-35 and Asp-41, and the second helix region comprises Ala-43, Leu-46, Glu-48, Ala-49, Ile-50, Tyr-52, and Leu-56. Notably, Fig. 2 demonstrated the high conservation of Leu-23 across the 170 PcbHLH amino acid sequences, underscoring its pivotal role in facilitating dimerization among PcbHLH proteins.

The basic region of the PcbHLH proteins is crucial for DNA-binding. As per Massari’s classification, these proteins are divided into two main groups based on sequence information within the basic region of the bHLH domains. The majority, comprising 165 PcbHLH proteins, were identified as DNA binders, while a smaller group of five was classified as non-DNA binders (Massari & Murre, 2000) (Table 1). The DNA-binding proteins are further classified into E-box binders (including G-box binder) and non-E-box binders. Residues conservation suggested that 144 PcbHLHs are likely to bind E-boxs, with 109 of these identified as G-box binders, contingent on the presence of His/Lys-5, Glu-9, and Arg-13 in the basic region.

Phylogenetic analysis of PcbHLH genes

To elucidate the evolutionary relationships between PcbHLHs and the AtbHLHs, a maximum-likelihood phylogenetic tree was constructed. The analysis of clade support values and tree topology led to the identification of 22 subfamilies. As shown in Fig. 3, each subfamily includes representatives from both A. thaliana and P. canescens, indicating a high degree of conservation in the bHLH domains through their evolution. The largest subfamily, Clade III, contains 37 members, whereas the smallest, Clade XIV, consists of only five members.

Analysis of gene structure and conserved motif of bHLH family

Motif analysis using MEME identified ten types of putative conserved protein motifs within the PcbHLH family, furthering our understanding of their conservation and diversity. As depicted in Fig. 4, members within the same subfamily exhibit similar motif structure, suggesting shared structural and functional characteristics. Gene structure analysis provided insights into the evolutionary relationships among PcbHLH family members. Of the 170 family members, 11 were found to lack introns and were clustered within two subfamilies. The remaining 159 members, which contained introns, showed significant structural similarities, reinforcing the notion of close evolutionary ties within their respective subfamilies.

Analysis of cis-acting elements in PcbHLH promoters

Cis-acting elements within the promoter regions are crucial for classifying subfamilies and functional characterizing members of the PcbHLH family. A 2,000 base pair (bp) upstream of the transcriptional start sites were identified as the promoter regions for this analysis. Several functionally significant cis-acting elements were identified and categorized into three groups: TC-rich repeats (cis-acting element involved in defense and stress responsiveness), LTR (cis-acting element involved in low-temperature responsiveness), MBS (MYB binding site involved in drought-inducibility) belong to stress-responsive. GCN4_motif (cis-regulatory element involved in endosperm expression), CAT-box (cis-acting regulatory element related to meristem expression), O2-site (cis-acting regulatory element involved in zein metabolism regulation) belong to plant development-related. AuxRR-core (cis-acting regulatory element involved in auxin responsiveness), P-box (gibberellin-responsive element), CGTCA/TGACG-motif (cis-acting regulatory element involved in the MeJA-responsiveness), ABRE (cis-acting element involved in the abscisic acid responsiveness), SARE (cis-acting element involved in salicylic acid responsiveness) belong to phytohormone responsive. Figure 5 illustrates that the majority of the PcbHLH members contain phytohormone-related cis-acting elements, implying their potential roles in various abiotic stress responses.

Gene duplication and synteny analysis of PcbHLHs

To discern the primary evolutionary forces shaping the PcbHLH gene family, gene duplication events within P. canescens were analyzed using TBtools. Figure 6A displays 92 segmental duplication gene pairs and eight tandem duplication gene pair (PcbHLH65/PcbHLH66, PcbHLH66/PcbHLH67, PcbHLH119/PcbHLH120, PcbHLH144/PcbHLH145, PcbHLH159/PcbHLH160, PcbHLH160/PcbHLH161, PcbHLH161/PcbHLH163, PcbHLH167/PcbHLH168). These findings suggest that gene duplication has been a significant mechanism in the evolution of PcbHLH genes. The Ka/Ks ratio for these duplicated genes, as detailed in Table S3, was consistently below 0.8, indicating a predominant purifying selection acting on the bHLH genes.

Comparative synteny analysis across species

To further elucidate the phylogenetic mechanisms underlying the PcbHLH family, we constructed syntenic maps comparing P. canescens with P. trichocarpa and Arabidopsis. A total of 144 orthologous bHLH gene pairs were identified between P. canescens and P. trichocarpa, and 93 orthologs between P. canescens and A. thaliana (Fig. 6B). The higher number of orthologous pairs between P. canescens and P. trichocarpa suggests a closer phylogenetic relationship compared to A. thaliana.

PcbHLH genes expression patterns in response to Cd stress

Prior research has established the role of certain clade IVc and Ib members as key regulators of Cd stress response in Arabidopsis (Hao et al., 2021). Leveraging these findings, we conducted a phylogenetic analysis of AtbHLHs and PcbHLHs to explore the functions of PcbHLH proteins. We meticulously selected 14 PcbHLH genes based on their functional homologs within the same subfamily for further analysis. Following Cd treatment, we observed transcriptional changes in these genes (Fig. 7A). In roots, genes such as PcbHLH148 and PcbHLH98 exhibited upregulated expression, while others like PcbHLH19 and PcbHLH50 showed downregulated expression. In stems, PcbHLH162 displayed decreased expression initially, with other genes upregulated up to 6 h. In leaves, all but PcbHLH61 exhibited upregulated expression, with a notable positive regulatory response observed for most PcbHLH genes, particularly in leaf tissues (Fig. 7B); PcbHLH96 was an exception, showing high expression levels in roots.