The published 158 bHLH genes in A. thaliana were used as a query sequence to run the BLAST (Basic Local Alignment Search Tool) program on the protein dataset of chestnut with default parameters. Furthermore, we used the latest Hidden Markov Model (HMM) file of the bHLH domain (PF00010) to search against the all proteins of chestnut using HMMER 3.0 (Evalue < 1e-5). After determine the existence of the bHLH domain using the NCBI Batch CD Search tool, we obtained 94 bHLH genes in chestnut, which were renamed CmbHLH1 to CmbHLH94 based their order on chromosome. The ExPasy website was used to obtain their physicochemical properties, and we used the eggNOG-mapper (Cantalapiedra et al., 2021) to obtain GO and KEGG annotation. We performed GO and KEGG enrichment analysis and visualization using TBtools (Chen et al., 2020). The subcellular localization of CmbHLH genes were predicted by using the online tool Cell-Ploc (Chou and Shen, 2010), and the full length coding sequence (CDS) without stop codon was cloned into the pAN580 and N-terminal fused with the green fluorescent protein (GFP) to verify the location of CmbHLH9, CmbHLH25 and CmbHLH55. The fusion constructs then transiently expressed into A. thaliana protoplasts.

MEGA 7.0 (Kumar et al., 2016) was used to construct a phylogenetic tree for the bHLH genes in chestnut with the maximum likelihood method, and we used the “Find Best DNA/Protein Models (ML)” function to find the best amino acid substitution model (partial deletion 95%). The final parameters were as follows: Jones–Taylor–Thornton (JTT) model; Gamma Distributed (G); Partial deletion 95%; 1,000 bootstrap replications. A phylogenetic tree of CmbHLH genes was constructed with neighbor-joining method using MEGA 7.0 with parameters: Poisson model, pairwise deletion, and 1,000 bootstrap replications. Additionally, Bayesian method was used to construct a phylogenetic tree of CmbHLHs. Briefly, MAFFT 7.0 (Katoh and Standley, 2013) and ModelFinder (Kalyaanamoorthy et al., 2017) were used to perform multiple sequence alignments and select optimum protein substitution model, respectively. MrBayes 3.2 (Ronquist et al., 2012) was used to construct Bayesian tree based on the best-fit model: JTT + F + I + G4. TBtools was used to visualize the structure of CmbHLH genes. The conserved domain and motif information of the CmbHLH genes were obtained from the Batch-CDD and MEME online programs, respectively. The combined analysis of phylogenetic tree, conserved domain and motif distribution was used to explore the sequence characteristics in different subgroups. The psRNATarget (Dai and Zhao, 2011) and Simple Sequence Repeat Identification Tool (SSRIT) on Gramene website (Tello-Ruiz et al., 2021) were used to predict miRNA targets and SSRs of the CmbHLH genes, respectively. The 3D structure models of all CmbHLH proteins were predicted through Swiss Model.

To further explore the evolution of the CmbHLH genes, the Multiple Collinearity Scan Toolkit (MCScanX) program (Wang et al., 2012) was used to perform a collinearity analysis within chestnut and with A. thaliana, rice, oak and grape genomes. Notably, we used a detailed method as in our previous study to determine the duplication model of the CmbHLH genes (Yu et al., 2022a). Based on the algorithm of the MCScanX program (Wang et al., 2012), the protein-coding genes from chestnut genome was compared against itself, and the duplicate genes were first defined as i) “dispersed duplicates”; ii) If the distance or gene rank between two duplicate genes is less than 20, they are re-labeled as “proximal duplicates”; iii) If the duplicate genes are adjacent (gene rank is 1), they are re-labeled as “tandem duplicates”; iv) The duplicated genes in the collinearity regions are re-labeled as “WGD or segmental ". The MCScanX program was run to directly obtain CmbHLH genes from proximal, dispersed and tandem duplication. Specifically, based on the collinearity results within the chestnut genome, we drawn the homologous collinearity dot-plot within the chestnut genome. Further, non-synonymous (Ka) and synonymous substitution sites (Ks) values of homologous gene pairs on the homologous collinear blocks were calculated by KaKs_Calculator (Wang et al., 2010), and the median Ks values was calculated by writing script (Yu et al., 2023). The median Ks values of homologous collinear blocks containing the CmbHLH genes were labeled on the dot-plots. Finally, based on the distribution of Ks values corresponding to the whole genome duplication (WGD) event that occurred in the chestnut genome in our earlier stage (Yu et al., 2023), combined with the complementarity of homologous collinear blocks, we finally determined the CmbHLH genes formed by WGD event.

Based on the structure annotation of the chestnut genome, the 2000 bp upstream sequences of 94 CmbHLH genes were extracted from the chestnut genome sequence. The upstream sequences of CmbHLH genes were submitted to PlantCARE software for identification of cis-acting elements and function categorization. We used Tbtools to visualize the predicted cis-acting elements, except for those members without specific function.

In order to explore the potential function of the CmbHLH genes in different tissues of chestnut, public database resources were utilized. RNA-seq data of chestnut buds at different stages (20, 25 and 30 days post-anthesis), chestnut nuts at different stages (70, 82 and 94 days post-anthesis), chestnut fertile/abortive ovules at different stages (15-July, 20-July and 25-July) were obtained from NCBI, and corresponding accession numbers were in Supplementary Table S1. The above transcriptome data from chestnut buds, nuts and fertile/abortive ovules were all three biological replicates. All RNA-seq reads were mapped on the reference genome (N11-1 from Castanea Genome Database) and the Kallisto software was used to quantify transcriptome data into Transcripts Per Kilobase of exon model per Million mapped reads (TPM). Tbtools was used to map the heatmaps for exhibiting the expression level of CmbHLH genes. The “Normalized” (scale method) in TBtools was used to normalize the expression, and the heatmaps were created by TBtools based on the transformed data of log₂ (TPM + 1) values.

The samples used for real-time quantitative PCR experiments were the nuts of “Yan Long” chestnut planted by Hebei Normal University of Science and Technology at 70, 82, and 94 days post-anthesis with three biological replicates. All samples were immediately frozen in liquid nitrogen and stored at −80°C for RNA extraction. The HiPure Total RNA kit (Magen, R4111, China) was used to extracted and isolated total RNA. PrimeScript RT Master Mix (Takara Biotechnology Co., Dalian, China) was used to reverse transcribe RNA into single stranded cDNA. Real-time quantitative PCR experiments were performed on ABI 7500 Real-Time PCR system (Applied Biosystems Inc., Foster City, CA, United States) with TB Green Premix Ex Taq (Takara). The relative gene expression values were calculated using the comparative CT (2^-△△CT) method (Livak and Schmittgen, 2001), and 18S gene of chestnut was used an internal control.

We identified 94 bHLH genes in chestnut genome using the BLASTP and HMMER 3.0 programs, and they were designated as CmbHLH1 to CmbHLH94 based on their relative position on chromosome and scaffold (Figure 1A; Supplementary Table S2). Eighty-eight CmbHLH genes were randomly distributed on all 12 chromosomes and other six CmbHLH genes were located on five unanchored scaffolds. Chromosome 8 had the highest number of CmbHLH members (13), while chromosome 6 had the fewest CmbHLH members (3) (Figure 1A). The analysis of all CmbHLH protein sequences showed that they had extensive variation in physicochemical properties (Supplementary Table S2). The length of CmbHLH proteins varied from 159 to 711 aa, molecular weight ranged from18.0 to 79.2 KDa. Most CmbHLH proteins were weakly acidic, because their average theoretical isoelectric point (pI) was 6.48. Only the instability index of CmbHLH42 protein was less than 40, and only the grand average of hydropathicity (GRAVY) value of CmbHLH49 protein was greater than zero, which meant that most CmbHLH proteins were unstable and hydrophilic. The aliphatic index of CmbHLH proteins varied from 49.33 to 109.53. A total of 6,831 microRNA (miRNA) targets were identified in CmbHLHs, indicating that corresponding miRNAs can activate or inhibit gene expression (Supplementary Table S3). In addition, 180 simple sequence repeats (SSRs) were predicted in the CmbHLHs, and the 3D structure models of CmbHLH proteins were predicted (Supplementary Tables S4, S5; Supplementary Figure S1). All these results provided valuable information reference for further research on the CmbHLH genes.

The classification of plant bHLH proteins had always been vague, but it was generally considered to be divided into 15–32 subgroups (Pires and Dolan, 2010b; Carretero-Paulet et al., 2010). We used the maximum likelihood method to construct an unrooted phylogenetic tree to explore the evolutionary relationships of CmbHLH genes (Figure 1B). The 94 CmbHLH proteins were classified into 19 subgroups based on the unrooted tree, which were named A to S (Figure 1B). Interestingly, subgroup J and K contained a single CmbHLH protein (CmbHLH16 and CmbHLH22), respectively. Except subgroup J and K, the number of CmbHLH genes in each subgroup varies greatly from 3 (subgroup D, H, Q and R) to 10 (subgroup E). In addition, we also constructed two phylogenetic tree of CmbHLH proteins with neighbor-joining and Bayesian method, respectively (Supplementary Figures S2, S3). The results showed that the phylogenetic trees of CmbHLHs constructed by the three methods were almost identical, with 94 CmbHLHs divided into 19 subgroups. Gene Ontology (GO) and The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed to characterize the potential functions of the identified CmbHLH genes. As shown in Figures 1C, D, CmbHLH genes were mainly involved in RNA biosynthetic process, positive regulation of gene expression, and transcription regulator activity, and these functions and processes were closely related to transcription factors. KEGG enrichment results indicated that they were mainly enriched in transcription factors, signal transduction, plant hormone signal transduction, environmental information processing, which were the main mechanisms by which bHLH family transcription factors regulate the expression of downstream genes.

The motif diversity and gene structure of CmbHLH genes were exhibited to reveal their evolution (Figure 2). The online MEME program was performed to identify motif patterns, and CmbHLH genes contained varying numbers of conserved motifs (Figure 2A). CmbHLH members in the same subgroup usually had similar motif arrangement patterns, while members in different subgroups had obvious differences. Although motifs 1 and 2 were present in almost all CmbHLH genes, some of the other motifs were detected in some specific subgroups. For example, motif 8 only can be detected in subgroup A, and motif 10 only can be detected in subgroup N. Interestingly, it was found that some CmbHLH genes in M subgroup had the same motif composition as the members in L subgroup, but some members of the M subgroup had completely different motif arrangements. This result implied that some CmbHLH genes had undergone dramatic changes in their non-motif sequences during evolution and formed the new subgroups without motif composition changing. Although the gene structure analysis indicated that there was extensive variation in the number of introns (0–11 introns) and gene structure of the CmbHLH genes, the members in the same subgroup had similar intron-exon structure (Figure 2B). For example, members in subgroup A had 0 or 1 intron, while members in subgroup Q had 3 introns. Interestingly, some CmbHLH members had exactly motif arrangement and extreme similar gene structure, such as CmbHLH39 and CmbHLH40, and subsequent analysis showed that they were formed from tandem duplication. These results indicated that the members in the same subgroup had similar motifs arrangement pattern and gene structure, which suggested that they may have similar functions.

The upstream 2,000 bp sequences of the CmbHLH genes were extracted to predict cis-acting elements in promoter regions, which can further explore the potential function and regulatory mechanism of the CmbHLHs (Figure 3; Supplementary Table S6). A total of 1,562 cis-acting elements were detected in promoter regions of CmbHLHs by the PlantCare software. Abundant cis-acting regulatory elements related to endosperm expression, meristem expression, and responses to gibberellin (GA) and auxin were identified in the upstream sequences of CmbHLH genes. This indicated that these genes may have potential functions in the morphogenesis of chestnut. Additionally, numerous cis-acting elements related to environment and stress were detected, such as cis-acting element involved in defense and stress responsiveness, cis-acting element involved in low-temperature responsiveness, cis-acting regulatory element essential for the anaerobic induction, and wound-responsive element. Notably, the cis-acting elements commonly considered important in relation to hormone response were detected, such as auxin-responsive element, gibberellin-responsive element, and salicylic acid responsive element, which were believed to be widely involved in all aspects of plant growth and development (Woodward and Bartel, 2005; Yamaguchi, 2008; Peng et al., 2021). The cis-acting elements involved in many functions were identified, which suggested that CmbHLH genes played diverse roles in chestnut.

The collinearity analysis among chestnut, grape, A. thaliana, oak, and rice genomes was performed to further investigate the evolution of the CmbHLH genes (Figure 4). Interestingly, we detected 68 CmbHLHs in the collinearity regions of the chestnut and grape genomes, and there were only 51, 50 and 25 CmbHLHs in the collinear regions between chestnut and A. thaliana, chestnut and oak, chestnut and rice, respectively (Supplementary Tables S7–S10). Additionally, there were 88 orthologous gene pairs contained CmbHLHs between chestnut and grape genome, and only 74, 57 and 44 orthologous gene pairs were found between chestnut and A. thaliana, chestnut and oak, chestnut and rice, respectively (Supplementary Tables S7–S10). As we all know, grape had a relatively stable genome structure, which was often used as a good reference for understanding the genome or related gene evolution of other eudicot plants (Consortium et al., 2012; Yu et al., 2022a). Therefore, the chestnut and grape genomes had not experienced additional whole genome duplication (WGD) events after the eudicot common hexaploidization event (ECH) (Yu et al., 2022b; Hu et al., 2022), and the grape genome had a relatively stable genomic structure, which resulted in the above statistical results about chestnut and grape. The distant genetic relationship between chestnut and rice, and rice genome had experienced additional WGD after ECH, may result in the least CmbHLHs in their collinear regions, and the least orthologous gene pairs contained CmbHLHs between them.

Significantly, we can detect some CmbHLHs existed in three or more collinear gene pairs, such as CmbHLH1, which suggesting that more copies of their orthologous genes in A. thaliana, rice, oak and grape genomes had been retained (Figure 4). According to the gene balance hypothesis, the products of these CmbHLHs may participate in macromolecular complexes or signaling networks (Birchler and Veitia, 2007). The subsequent transcriptome analysis revealed that they may have important functions in the bud development and ovule fertility of chestnut. Additionally, the collinear blocks contained CmbHLHs generally including more gene pairs in chestnut and grape (average of 49 gene pairs) (Supplementary Table S11). Contrasty, there were only 19, 27 and 11 (average value) gene pairs per block contained CmbHLHs between chestnut and A. thaliana, between chestnut and oak, between chestnut and rice, respectively (Supplementary Tables S12–S14). We speculated that the differences in WGD experience and genetic relationships of these species may be the cause of these results.

Gene duplication event was an important motivation for the expansion of gene families and the occurrence of new functions (Liu et al., 2018; Quan et al., 2019; Yu et al., 2021). The Multiple Collinearity Scan toolkit (MCScanX) program (Wang et al., 2012) and writing script were used to conduct collinearity analysis of chestnut genome to determine the duplication model of CmbHLHs (Yu et al., 2021). MCScanX could identify most duplication types except WGD and segmental duplication, and the reason was that the plant genome undergone a large number of chromosome recombination after WGD, which made them difficult to distinguish (Figure 5) (Wang et al., 2018; Yu et al., 2022a; Wang et al., 2022; Yu et al., 2023). As in our previous study (Yu et al., 2023), by drawing the homologous collinearity dot-plot within chestnut genome labeled with the median Ks of the homologous blocks, combining their length and complementarity, WGD and segmental duplications were finally distinguished (Figure 6). The results indicated that most of the CmbHLHs formed from dispersed duplication, followed by WGD, and no CmbHLH member was considered as singletons (Supplementary Table S15). These results showed that dispersed duplication was the main reason for the increase in the member of the bHLH gene family in chestnut. Additionally, Ka/Ks values of homologous gene pairs were calculated, which indicated that all of CmbHLHs were experienced purifying selection pressure during evolution (Ka/Ks < 1) (Supplementary Table S16).

The transcriptome data of chestnut buds at different stages from NCBI were investigated to explore their potential function in chestnut bud development (Figure 7A). Some CmbHLH genes had high expression level at the all stages of chestnut bud development, suggesting their potential function in the chestnut bud development process. For example, the average TPM values of CmbHLH88 were 138.16, 165.03 and 155.20 at three stages, respectively. Notably, the expression level of some CmbHLH genes decreased (CmbHLH25) or increased (CmGRAS51) sharply as the buds continue to growth. The TPM value of CmbHLH25 in chestnut buds at 20 days post-anthesis was 107.91, while the TPM values at 25 and 30 days post-anthesis were 27.07 and 3.44, respectively. The TPM value of CmbHLH51 in chestnut buds at 20 days post-anthesis was 34.37, while the TPM values at 25 and 30 days post-anthesis were 71.61 and 109.33, respectively. These results indicated that CmbHLHs play an important role in chestnut bud development.

In addition, we analyzed transcriptome data from fertile/abortive ovules of chestnut at different developmental stages to characterize the possible role of the CmbHLH genes in ovule fertility (Figure 7B) (Du et al., 2021). We found that the expression levels of a few CmbHLHs in fertile/abortive ovules showed significant differences, such as CmbHLH55 and CmbHLH69. The expression level (TPM of all samples was more than 200) of CmbHLH55 in fertile ovules was significantly higher than that in abortive ovules at all three stages, while the expression of CmbHLH69 of fertile ovule (TPM of all samples were less than 50) was significantly lower than that of abortive ovules at all three stages. Additionally, the expression levels of some CmbHLH members had significantly changed during the growth of fertile ovules. For example, the TPM value of CmbHLH69 on 15-July was 86.24, while the TPM values on 20-July and 25-July reached 151.44 and 144.10, respectively. Notably, multitudinous of gibberellin (GA) responsive elements had been detected in their promoter regions, and GA was considered to be a plant hormone that was extremely important for plant reproductive organs (Figure 3) (Plackett et al., 2014). Cis-regulatory element involved in endosperm expression had also been identified in the promoter regions of these CmbHLHs, which indicated that they may have impact on the fertility or development of chestnut ovules.

Similarly, transcriptome data from chestnut nuts at different developmental stages were analyzed (Figure 7C) (Li et al., 2021). Obviously, the expression of some CmbHLH genes at 94 days post-anthesis of chestnut nut development had reduced dramatically, such as CmbHLH25 and CmbHLH78. For example, the average TPM values of CmbHLH25 were 274.67 at 70 days post-anthesis, and this values were 164.14 and 1.29 at 82 and 94 days post-anthesis, respectively. Considering that some bHLH genes had been verified to be involved in the regulation of secondary metabolite formation (Li et al., 2006) and that some cis-acting regulatory elements involved in the regulation of zein metabolism and flavonoid biosynthesis had been discovered (Figure 3; Supplementary Table S6), we believed that these CmbHLH members may be related to chestnut nuts ripening and related quality. To verify the expression patterns of CmbHLHs, six CmbHLH genes were analyzed using qRT-PCR, and the results were shown in Figure 8A. Overall, the qRT-PCR experiments analysis results were consistent with the RNA-seq analysis. For example, CmbHLH25 was significantly highly expressed in nuts at 70 days post-anthesis than in nuts at 82 and 94 days post-anthesis. The expression level of CmbHLH2 was highest at 70 days post-anthesis, and decreased to the lowest at 82 days post-anthesis, and increased at 94 days post-anthesis. The expression trends of other four CmbHLH genes analyzed using qRT-PCR during chestnut nut development were also consistent with the results of RNA-seq analysis. In addition, almost all CmbHLH proteins (85/94) were predicted in the nucleus, based on subcellular localization analysis (Supplementary Table S2). We instantaneously expressed GFP-CmbHLH9, GFP-CmbHLH25, and GFP-CmbHLH55 fusion proteins in A. thaliana protoplasts for subcellular localization. As shown in Figure 8B, CmbHLH25 and CmbHLH55 were co-localized to both the nucleus and cytoplasm, and CmbHLH9 was localized to the nucleus.