The Chinese chestnut (Castanea mollissima) were collected from Yongren County, Chuxiong Prefecture, Yunnan Province, China. There are six varieties, namely ‘Yongfeng1’ (YF1), ‘Yongren zao’ (YRZ), ‘Yongfeng 2’ (YF2), ‘Yimen 2’ (YM2), ‘Yimen 1’ (YM1), and ‘Yanlong’ (YL), respectively. Developing spiny bur with nuts were collected at 10-day intervals, starting at 20 days until spiny bur matured into ripe fruit with cracked bulbs. Plant tissues were collected and immediately frozen in liquid nitrogen. There were three periods in total. The Chinese chestnuts collected in this study were cultivars, and Southwest Forestry University approved the experiment. The formal identification of the plant material was undertaken by Prof. Zhiheng Zhao. Our experimental research complied with local legislation, national and international guidelines.

To identify the bHLH genes in Chinese chestnut, the Hidden Markov Model (HMM) profile of the Pfam bHLH domain (PF00010) was downloaded from the pfam database (http://pfam.xfam.org/), and HMMER software (v3.1b2) (Potter et al., 2018) was performed to find the CmbHLH genes against in the Chinese chestnut genome (https://www.ncbi.nlm.nih.gov/search/all/?term=PRJNA912750). All CmbHLH candidates were further confirmed in the CD-search (NCBI database) (Marchler-Bauer et al., 2015) to verify the bHLH domain existing in each sequences, and any redundancy was manually removed. The conserved motifs of CmbHLH proteins were identified in Multipe EM for Motif Elicitation (MEME, http://meme-suite.org/).

The bHLH sequences of Arabidopsis thaliana and Zea mays were downloaded from PlantTFDB (http://planttfdb.cbi.pku.edu.cn/). Using the ClustalX in MEGA X for multiple sequence alignment of constructed with 124 CmbHLHs, 225 AthbHLHs, 323 TabHLHs, and 308 ZmbHLHs (Meng et al., 2020), and Maximum Likelihood method was used to construct a phylogenetic tree with 1000 bootstrap sampling iterations in MEGA11 (Kumar et al., 2018). The amino acid sequences of bHLH proteins were aligned using MAFFT, and a phylogenetic tree was constructed using MEGA X based on the JTT substitution model. The resulting tree was visualized and annotated using iTOL (https://itol.embl.de/).

The 2,000 bp upstream promoter sequences of 124 chestnut CmbHLH transcription factor genes were extracted using TBtools. Cis-acting elements were identified via the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The elements were subsequently filtered and grouped based on similar regulatory functions. A heatmap was generated using TBtools, and the number of cis-acting elements in each functional category was statistically analyzed.

The exon–intron structures of the identified genes were analyzed using the Gene Structure Display Server (GSDS) v2.0 (http://gsds.gao-lab.org/). The coding sequences (CDS) and corresponding genomic sequences were aligned to determine exon and intron positions. The resulting gene structure models were visualized to examine the structural diversity among gene family members.

In this study, the dynamic changes of CmbHLH genes in six chestnuts varieties were analyzed during fruit development and gall-information. The RNA-seq data of six chestnuts varieties flower and fruit development came from our previous study (PRJNA574282). And the RNA-seq data of chestnut in response to Dryocosmus kuriphilus Yasumatsu were downloaded from NCBI (PRJNA512447) (Zhu et al., 2019). Clean data were obtained from raw sequencing reads using fastp v0.20.1 (Chen et al., 2018), which removes adapter sequences and trims reads containing poly-N and low-quality bases (reads with more than 50% of bases having a Phred quality score ≤20). Clean reads from all samples were mapped to the C. mollissima genome (https://www.hardwoodgenomics.org/Genome-assembly/1962958) by using HISAT2 software for further analysis (Kim et al., 2019). Then all transcripts were quantified by Salmon 1.4 (Patro et al., 2017). Gene expression profiling was performed using FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values. The raw FPKM values generated by Salmon v1.4 were used directly for expression pattern analysis without further normalization, as FPKM inherently accounts for transcript length and sequencing depth. Gene co-expression networks were constructed based on FPKM expression values using the WGCNA (Weighted Gene Co-expression Network Analysis) R package (v1.71). An appropriate soft-thresholding power was selected based on the scale-free topology criterion. Modules of co-expressed genes were identified using hierarchical clustering and dynamic tree cutting. The module–trait relationships were calculated to identify modules associated with specific conditions or developmental stages. The resulting networks were visualized using Cytoscape (v3.9.1).

124 ZjbHLH protein sequences were used as queries, and protein-protein interactions were predicted by the STRING website (https://string-db.org/). The orthologs of Arabidopsis thaliana and Zea mays were selected as references. Finally, an interaction network among CmbHLHs was constructed in this study.

Total RNA was extracted using a plant RNA extraction kit (Vazyme Biotech) according to the manufacturer’s instructions. The construction of a cDNA library was performed by reverse transcribing 1 mg RNA samples using the 5×HiScript^® Reverse Transcriptase (vazymes) and 4×gDNA (vazymes) kits, following the manufacturer’s protocol. The expression of representative genes was analyzed by qRT-PCR in at least four biological replicates. The primer sequence is in Supplementary Table S1 .

Actin was chosen as the reference genes. Relative expression levels of qRT-PCR data were calculated using the 2^-ΔΔCt method, and data are presented as mean ± standard deviation (SD). Statistical significance was tested using one-way ANOVA, corrected for multiple comparisons using Tukey’s HSD. A P value < 0.05 was considered statistically significant. All statistical analyses were performed using GraphPad Prism software.

A total of 124 non-redundant bHLH proteins were identified in the Castanea mollissima genome based on the Hidden Markov Model search and confirmed that them have typical bHLH domain. These 124 putative bHLH protein were renamed as shown in Supplementary Table S1 . The encoded bHLH protein ranged from 86 (CmbHLH94) to 1241 (CmbHLH19) amino acids (aa), with Open Reading Frame (ORF) ranging from 258 bp to 3723bp in length. Based on the different physical and chemical properties of CmbHLH, it is suggested that they are likely have multiple functions.

Based on the maximum likelihood (ML) phylogenetic tree in CmbHLH proteins, 124 CmbHLH were divided into eight subgroups and were named as group I to VIII. Among them, CmbHLH_I and CmbHLH_VIII had the fewest members, with seven and eight genes respectively. Within each subgroup, there was a high level of concordance in gene length and physicochemical properties among members of the same gene family, implying their potential co-localization on either the same or adjacent chromosomal loci ( Figure 1 ). The phylogenetic analysis suggests that CmbHLH_I and CmbHLH_VIII potentially represent the most ancient subgroups in C. mollissima, while the gene family CmbHLH_II could be considered more recent ( Figure 1 ). Additionally, the CmbHLH_VI subgroup exhibited the highest degree of differentiation and expansion, suggesting a lower level of conservation compared to the other subgroups, demonstrating that CmbHLH_VI is a dynamically evolving bHLH subgroup in C. mollissima ( Figure 1 ). In order to futher explore the function and classification of CmbHLH protein, the ML phylogenetic tree was constructed with 124 CmbHLHs, 225 AthbHLHs from Arabidopsis thaliana, 323 TabHLHs from Triticum_aestivum, and 308 ZmbHLHs from Zea mays ( Figure 2 ). The CmbHLH proteins were divided into 8 subgroups, named from I to VIII, according to the classification of Arabidopsis thaliana and Zea mays. The phylogenetic tree showed that the majority of orthologous bHLH genes are closely related in all three species, which indicates that most bHLHs arose prior to species divergence. CmbHLHs, which are relatively distant among species and relatively closely related within species, may have originated from the duplication event of paralogous genes.

The 124 CmbHLHs were mapped to 12 different chromosomes and 4 different contigs ( Figure 3 ). Among them, the CmbHLHs genes located on chromosomes GWHANWH00000079 and GWHANWH00000090 were the most, with 19 and 14 bHLHs genes, respectively ( Figure 3 ). Interestingly, some genes are tightly clustered together. Interestingly, some of the genes were tightly packed into clusters to form tandem repeats (CmbHLH68 and CmbHLH99; CmbHLH47 and CmbHLH52), most of these genes belong to the same subfamily. A previous study analyzing repetitive events in rice and Arabidopsis indicated that some bHLH subfamily members are most likely derived from repetitive events. Also, analysis of the gene structure of all CmbHLHs shows that the bHLHs gene structure of C. mollissima contained 1–5 motifs, indicating that the gene family was highly conserved in terms of gene structure ( Figure 4 ). it showed that motif-1 and motif-2 ware typical domains within the bHLH gene family in C. mollissimat. The subfamily of CmbHLHs VIII, with the exception of CmbHLH_106, exclusively contains motif-1 and motif-2, providing further evidence for the high conservation of CmbHLHs VIII within the C. mollissimat bHLH gene family ( Figure 4 ).

Multiple sequence alignment and MEME analysis showed that all CmbHLHs genes share the same conserved domains motif-1 and motif-2 ( Figure 4 ). In addition, 18 CmHLHs genes such as CmbHLH63, CmbHLH60, CmbHLH74, CmbHLH128, CmbHLH44, CmbHLH80, CmbHLH79, CmBEE3, CmbHLH59a, CmbHLH59b, CmbHLH122, CmbHLH62a, CmbHLH130a, CmbHLH137, CmbHLH62b, CmbHLH130b, CmbHLH31 and CmbHLH49 contain motif- 3 ( Figure 4 ). The structure and sequence of the CmHLHs gene family are highly conserved, with motif-1 and motif-2 forming the complete bHLH domain in the CmbHLHs genes family. Five residues (motif-1: Ala-1, Glu-2, Arg-3, Arg-5, Arg-6), five residues (motif-1: Arg-12, Leu-16, Leu-19, Val-20, Pro-21), two residues (motif-2: Gly-1, Leu-3), and six residues (motif-2: Ala-9, Leu-10, Leu-13, Leu -15, Val-17, Ala-20) constitute the basic region of the bHLH domain, the first helical region, the loop region and the second helical region ( Figure 5 ), respectively. Motif-3,4,5 are also highly conserved domains. Of the 21 amino acid sequence, 15 residues are more than 50% amino acid conserved. Seven of these residues (motif-3: Arg-9, Ala-10, Arg-12, Gly-13, Ala-15, Thr-16, Ser-20) are 100% amino acid conserved suggesting that the 18 CmbHLHs containing motif-3 may be a very important and functionally similar transcription factor in chestnut.

Cis-acting element analysis was performed on the bHLH gene family in chestnut to explore the potential functions of CmbHLH transcription factors. All 124 CmbHLH transcription factors contained cis-elements in their promoter regions related to responses to biotic and abiotic stresses, growth and development, and hormone signaling ( Supplementary Figure 1 ). The promoter regions of the CmbHLH transcription factors harbored multiple anaerobic induction elements (ARE), light-responsive elements, and drought-responsive elements (MBS). Notably, CmbHLH_77, CmbHLH_15, and CmbHLH_83 contained a relatively high number of ARE elements, while 60 CmbHLH transcription factors possessed a larger number of MBS elements.

Among common hormone-responsive elements, CmbHLH_35 contained 10 methyl jasmonate-responsive elements (TGACG-motif), and both CmbHLH_4 and CmbHLH_55 contained 8 TGACG-motifs. Additionally, 45 CmbHLH transcription factors harbored multiple salicylic acid-responsive elements (TCA-element), while CmbHLH_35, CmbHLH_46, and CmbHLH_54 possessed multiple abscisic acid-responsive elements. In terms of growth and development, 45 CmbHLH family members contained growth-regulating G-box elements, and 9 CmTrihelix family members included meristem expression-related CAT-box elements. These findings suggest that CmTrihelix transcription factors play crucial roles in plant responses to abiotic stress, hormone signal transduction, and growth and developmental regulation.

Based on the orthologs in Arabidopsis and maize, STRING predicted interactions between CmbHLHs genes. Several important interactions are predicted in Figure 6 , where 100282922 (CmbHLH122) may be involved in cuticle development, regulation of stomatal movement, photoperiodism and flowering. And bHLH29 (CmbHLH29) may participate in the regulation of iron ion transport and the response to iron ion. While genes such as AT2G40200 (CmbHLH57), SPCH (CmbHLH98), bHLH104 (CmbHLH104), and AT2G14760 (CmbHLH84) are grouped into one module on the PPI network diagram. And it is speculated that the function of these genes may be related to stoma formation and regulation of osmotic pressure. The concentration of genes such as AT1G31050 (CmbHLH111), AT5G50915 (CmbHLH137), and BPEp, and UNE12 in one module may indicate that these transcription factors may be related to petal development and gamete formation. These results further proved the functional diversity of CmbHLHs genes. The predicted network provides some useful clues for functional studies, but further experimental evidences is needed.

Previous studies have shown that genes of bHLH transcription factor family genes are widely involved in plant development, fruit ripening, and abiotic stresss (Atchley et al., 1999; Massari and Murre, 2000; Carretero-Paulet et al., 2010; Zhang et al., 2018). To identify CmbHLHs genes that have a significant role in chestnut flower and fruit development, transcriptome data of 6 chestnut varieties at different stages of fruit development were analyzed. The results showed that 84 out of the 124 CmbHLHs genes were differentially expressed during chestnut fruit development (log2|FC|>=1, p<0.05) ( Figure 7 ). Among them, CmbHLH4, 7, 16, 38, 59, 64 and other genes were mainly expressed in the early and middle stages of chestnut fruit development and most of these genes belonged to subgroup V and VI. While genes such as CmbHLHs 11, 13, 15, 25, 93, 101, 122 were mainly expressed in the late stage of chestnut fruit development. The trends of these genes in different varieties were basically the same in all 6 varieties. In general, the expression of CmbHLHs exhibited stage-specific patterns during the maturation process of chestnut kernels. Among them, the expression of CmbHLHs in YF1, YM1, and YRZ showed similarity. There were a large number of CmbHLH-specific high expressions in stage III, and only a few CmbHLH-specific high expressions in stage I and II. YM2 and YL showed different expression patterns. CmbHLHs showed more specific high expression in stage I, but less in stage II. The results show that these genes are essential transcription factors for fruit development in chestnut, and that they play an important role in fruit development.

A co-expression network was constructed based on the correlation of expression levels of bHLH and other genes during fruit development. The network contained 1045 nodes (96.17% of other genes and 3.83% of bHLH transcription factors) and 1272 edges (87.74% of pairs were positive correlated and 12.26% were negative correlated) ( Figure 8 ). In general, we consider the expression of transcription factors and target genes to be positively correlated. The two most highly connected transcription factor are CmbHLH82 and CmbHLH57, which are positively correlated with the vast majority of other genes. From the expression patterns of CmbHLH82 and CmbHLH57 in YF1, YRZ and YM1, it can be seen that their expression decreases with the increase of development, which indicates that transcription factors are crucial in the gene initiation stage.

In order to further analyze the function of CmbHLHs genes in resistance to D. kuriphilus infestation, RNA-seq data of CmbHLHs genes from our previous studies were used to investigate the dynamic changes during four leaf gall-formation stages caused by D. kuriphilus. We identified 15 CmbHLHs genes that were progressively less expressed with D. kuriphilus infestation. These genes included CmbHLH164, CmbHLH42, CmbHLH163, CmbHLH124a, CmbHLH96a, CmbHLH60. They continued to be down-regulated in expression under the influence of D. kuriphilus infestation. At the same time, 21 CmbHLHs genes showed an up-regulated trend in expression with the infested of D. kuriphilus. These genes include CmbHLH63, CmbHLH51. Overall, most CmbHLHs genes showed down-regulated expression when infested by D. kuriphilus. The persistent down-regulation of CmbHLH42 may indicate that D. kuriphilus affects downstream anthocyanin biosynthesis by affecting the expression of CmbHLH42. In addition, genes such as CmbHLH63 and CmbHLH51 showed up-regulated expression when chestnut leaves are infected by D. kuriphilus. This may indicate that CmbHLH63 and CmbHLH51 may play an important role in the process of chestnut resistance to biological stress. The down-regulated expression of these genes may be closely related to the formation of galls in chestnut leaves when they are invaded by D. kuriphilus ( Figure 9 ).

To further evaluate the potential functionality of these genes, five genes from eight subfamilies were selected and performed qRT-PCR analysis to examine the expression patterns of these representative genes throughout the four stages of chestnut seed development ( Figure 10 ). The CmbHLH genes showed different expression patterns at different stages of seed development, suggesting its diverse regulatory roles. All genes were expressed in different stages, and five genes were significantly up regulated at I_stage and peaked at II_stage. Subsequently, during the III_stage and IV_stage, there was a declining trend in gene expression levels. However, they remained significantly higher than the control level. These five representative genes play a pivotal role in seed development ( Figure 10 ).