Mechanism of tobacco hairy root development and functional study of the HSF gene family based on comparative transcriptomics - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice BMC Plant Biol . 2026 Mar 9;26:685. doi: 10.1186/s12870-026-08531-9 Search in PMC Search in PubMed View in NLM Catalog Add to search Mechanism of tobacco hairy root development and functional study of the HSF gene family based on comparative transcriptomics Xiaozong Wu Xiaozong Wu 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Xiaozong Wu 1 , Zhitao Qi Zhitao Qi 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Zhitao Qi 1 , Lei Wang Lei Wang 2 China Tobacco Yunnan Industrial CO., LTD. , Kunming, 650000 China Find articles by Lei Wang 2 , Yaolan Zhuang Yaolan Zhuang 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Yaolan Zhuang 1 , Yixuan Xue Yixuan Xue 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Yixuan Xue 1 , Haoyu Yan Haoyu Yan 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Haoyu Yan 1 , Jifan He Jifan He 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Jifan He 1 , Peilin Li Peilin Li 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Peilin Li 1 , Zhiwen Zhu Zhiwen Zhu 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Zhiwen Zhu 1 , Guiliang Tang Guiliang Tang 3 Department of Biological Sciences, Michigan Technological University, Houghton, MI USA Find articles by Guiliang Tang 3 , Meng Li Meng Li 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Meng Li 1, ✉ , Chaonan Shi Chaonan Shi 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China Find articles by Chaonan Shi 1, ✉ Author information Article notes Copyright and License information 1 International Joint Research Laboratory for Utilization of Plant Functional Components, College of Tobacco Science and Engineering, Zhengzhou University of Light Industry, Zhengzhou, 450001 China 2 China Tobacco Yunnan Industrial CO., LTD. , Kunming, 650000 China 3 Department of Biological Sciences, Michigan Technological University, Houghton, MI USA ✉ Corresponding author. Received 2025 Oct 1; Accepted 2026 Mar 5; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13085611 PMID: 41803700 Abstract Hairy roots induced by Agrobacterium rhizogenes exhibit unique phenotypic traits such as rapid proliferation, hormone-independent growth. To further elucidate the transcriptomic differences between hairy roots and the roots, stems, and leaves of normal plants, as well as their potential mechanisms of formation, this study conducted a comparative transcriptomic analysis. The results indicated that differentially expressed genes (DEGs) were significantly enriched in pathways related to energy metabolism, stress response, and secondary metabolism. Then weighted gene co-expression network analysis (WGCNA) analysis identified two key co-expression modules. Moreover, among the differentially expressed transcription factors (TFs), the HSF (heat shock transcription factor) family showed the most significant enrichment. Further analysis of 17 HSFs revealed that all of them were significantly correlated with plant hormones. Among these, two were specifically expressed in hairy roots, and nine exhibited both constitutive high expression and responsiveness to salicylic acid (SA), melatonin and drought treatments. These results establish an HSF-mediated regulatory network, offering valuable theoretical insights and potential targets for breeding stress-resistant germplasm with enhanced root systems, as well as for the synthetic biology-based engineering of hairy roots. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s12870-026-08531-9. Keywords: Differentially expressed genes (DEGs), Gene expression pattern, Gene regulatory network, Heat shock transcription factors (HSFs) family, Tabacco hairy roots Introduction The differentiation and development of plant roots, stems, and leaves are precisely regulated by complex genetic networks. Studies have shown that multiple important gene families play critical roles in these regulatory process. In Arabidopsis , the DA1 protein cleaves WUS to regulate shoot meristem function and organ size [ 6 ]. Similarly, in maize, the TSH4 gene defines meristem fate through a dual negative feedback mechanism and has been implicated in domestication [ 11 ]. During root development, PLT transcription factor (TF) gradients help establish the stem cell niche, while factors such as BRAVO, PLT3, and WOX5 collectively maintain root apical stem cell homeostasis in Arabidopsis [ 34 ]. Root morphogenesis is further modulated by ARF and IAA family genes through auxin signaling pathways [ 32 ]. Moreover, heat shock transcription factors (HSFs) have been shown to participate in the regulation of crown root development in rice [ 47 ]. These findings highlights the evolutionary conservation and diversification of key regulatory genes. By comparing gene expression profiles across different tissues, developmental stages, and stress conditions, researchers can systematically identify key regulatory genes involved in plant growth and development, stress responses, and secondary metabolism. Transcriptomic approaches have elucidated gene functions across species-such as OsWOX11 -mediated crown root development in rice [ 47 ], nitrogen (N) uptake-related ZmNRT1.1B in maize [ 10 ], and the roles of MYB/bHLH TF families in Arabidopsis anthocyanin biosynthesis [ 28 ]. Additionally, in tobacco, transcriptomic studies following inoculation with Phytophthora nicotianae revealed that methyl ferulate, a key root exudate whose content is regulated by NtCOMT10 , plays a critical role in plant defense by rapidly enriching rhizosphere Bacillus populations to enhance resistance against pathogen infection [ 25 ]. These studies collectively demonstrate the broad application of comparative transcriptomics in elucidating the genetic basis of plant traits. Compared to normal root systems, hairy roots induced by the transformation of rol genes from the Ri plasmid T-DNA of Agrobacterium rhizogenes ( A. rhizogenes ) exhibit rapid proliferation due to abnormal cell differentiation, a highly branched morphology, and the capacity for hormone-autonomous growth [ 44 ]. The rolB gene, a key component of the A. rhizogenes T-DNA region, is crucial for transformation efficiency, with its high expression triggering abundant hairy root (HR) formation alongside contributions from rolA , rolC , and rolD [ 46 ]. This makes the HR system an ideal model for studying root development and secondary metabolism, with applications in producing compounds like astragalosides [ 15 ]. Transcriptomic analyses have proven instrumental in dissecting the underlying regulatory networks. For instance, studies in peanut HR identified AhHDA1’s role in modulating flavonoid and phenylpropanoid pathways [ 31 ]. Similarly, study in tobacco HR revealed that MPK4 overexpression enhances nicotine accumulation [ 23 ], and identified auxin response factors (ARFs) NtARF7 and NtARF19 as key mediators of rolB gene signaling [ 4 ]. These studies have advanced our understanding of root differentiation and identified target genes for engineering hairy roots to enhance medicinal compound production via synthetic biology. Therefore, this study tests the hypothesis that hairy root development is driven by specific transcriptomic reprogramming events involving auxin signaling and secondary metabolic pathways. The aim is to systematically investigate differentially expressed genes (DEGs) and regulatory networks between hairy roots and normal tissues, with objectives to elucidate associated pathways, identify core modules via weighted gene co-expression network analysis (WGCNA), and characterize key TFs. Collectively, this research elucidates the systemic impact of A. rhizogenes -mediated transformation on plant gene expression, laying a molecular foundation for optimizing hairy root culture systems and enhancing bioactive compound production. Materials and methods Plant material preparation Tobacco ( Nicotiana tabacum cv. K326) seedlings were cultivated in trays under controlled conditions (24 ± 2 °C, 16 h light/8 h dark) for one month. Healthy, pest-free seedlings were selected, and leaf segments (approx. 1 cm × 1 cm) were excised as explants. These were surface-sterilized under aseptic conditions by sequential treatment with 75% ethanol (1 min) and 5% sodium hypochlorite (5 min), followed by three rinses with distilled water prior to Agrobacterium -mediated infection. Induction of tobacco hairy roots A. rhizogenes strain A4 (stored in the laboratory of the Tobacco Science and Engineering College of Zhengzhou Light Industry University) was used for hairy root induction. The strain was cultured in LB liquid medium containing 50 mg/L kanamycin (Kan, Bomeibio, China) and streptomycin sulfate (Bomeibio, China) with shaking (28 °C, 150 r/min). Bacterial density was monitored using an ultraviolet spectrophotometer (Shanghai Yili Analytical Instrument Co., Ltd., China) at 600 nm until OD600 reached 0.6–0.8. The culture was then centrifuged (4 °C, 4000 r/min, 15 min), the supernatant discarded, and the pellet resuspended in 1/2 MS liquid medium to obtain the infection suspension. Sterilized explants were immersed in this suspension for 10 min, blotted dry with sterile filter paper, and placed on 1/2 MS solid medium for co-culture (24 °C, dark, 48 h). Following co-culture, explants were rinsed three times with 1/2 MS liquid medium containing 100 mg/L cefotaxime (Cef, Bomeibio, China), then four times with sterile water, and transferred to 1/2 MS solid medium supplemented with 100 mg/L Cef. Culture conditions were identical to those for tobacco seedlings. Hairy roots emerged after approximately two weeks. All procedures were performed under sterile conditions. Extraction of RNA Samples were taken from tobacco hairy roots (R), as well as the roots (G), stems (J) and leaves (Y) of the same batch of tobacco seedlings. Each group of samples had three biological replicates. After sampling, the samples were immediately placed in liquid nitrogen and then stored at -80℃. Total RNA was extracted using a RNA extraction kit (TRIzol, Thermo Fisher Scientific). The concentration and purity of the RNA samples were detected using the Nanodrop ND-2000 instrument (Thermo Fisher Scientific, USA), and the RNA integrity index (RIN) was determined using the Agilent Bioanalyzer 5300 (Agilent Technologies, CA, USA). Construction and sequencing of cDNA library Qualified RNA samples are adsorbed onto Oligo(dT) magnetic beads to enrich mRNA. They are randomly fragmented into fragments of approximately 300 base pairs (bp), and then a six-base random primer is added. Under the action of reverse transcriptase, the mRNA was used as a template to synthesize the first strand cDNA, and then the second strand was synthesized to form double-stranded cDNA. After double-stranded purification, the ends were repaired and capped with A tails, and finally, PCR amplification was performed to obtain the cDNA library. The library was subsequently sequenced on the Illumina NovaSeq 6000 platform. Transcriptome data analysis Raw sequencing data were quality-controlled using Fastp (v0.19.5) [ 9 ]. Clean reads were aligned to the reference genome ( https://solgenomics.net/organism/Nicotiana_tabacum/genome ) [ 12 ] with Hisat2 (v2.1.0) [ 18 ] using default parameters. Gene functions were annotated via eggNOG ( http://eggnog5.embl.de/#/app/home ) [ 16 ], Gene Ontology (GO) ( http://www.geneontology.org ) [ 1 ], and Kyoto Encyclopedia of Genes and Genomes (KEGG) ( http://www.genome.jp/kegg/ ) [ 19 ]. Transcripts were assembled using Cufflinks (v2.2.1) [ 35 ] or StringTie (v2.1.2) [ 26 ] and compared with known transcripts to identify novel transcripts. Expression levels were quantified as fragments per kilobase per million mapped reads (FPKM). Differential expression analysis between groups was performed using DESeq2 (v1.24.0) [ 21 ], DEGseq (v1.38.0) [ 42 ], or EdgeR (v3.24.3) [ 30 ]. Genes with p -adjust < 0.05 and |log2FC| ≥ 1 were considered DEGs. Results were visualized as bar charts, Venn diagrams, and volcano plots via the Majorbio online platform ( https://www.majorbio.com/ ) [ 29 ]. GO and KEGG enrichment analysis The GO database is mainly used for classifying and annotating genes and their products. Based on this database, genes can be classified and annotated in CC (cellular component), MF (molecular function), and BP (biological process). The KEGG database enables the linkage of individual gene sequences to biological pathways. The DEGs screened out were subjected to GO and KEGG enrichment analysis through the Majorbio online platform. The method used was Fisher exact test. When the corrected p -value was less than 0.05, it was considered that this function was significantly enriched, and the results were visualized to explore the related functions of the DEGs. WGCNA WGCNA is a system biology method based on gene expression levels [ 20 , 27 ]. In this study, WGCNA was performed using the OECloud online tools ( https://cloud.oebiotech.com ) based on the DEGs that overlapped in the three comparison groups. Filter out genes with low fluctuation in expression variation (standard deviation ≤ 0.5). In order to approximate the conditions of a scale-free network distribution, this study employed an adaptive method to select the weight parameter power for the adjacency matrix. Following data import into the system, the power value was ultimately determined to be 30. TF family analysis To identify key TF families regulating hairy root formation, all differentially expressed TFs were quantified and enriched; the most significantly enriched family was selected for further analysis. Gene function annotation was performed via eggNOG-mapper ( http://eggnog-mapper.embl.de/ ) [ 8 ] using default parameters. GO annotation files (latest version) were obtained with TBtools (v2.332) [ 7 ] and integrated with functional annotation results. Gene set enrichment analysis (GSEA) using the hypergeometric test was conducted on the most enriched TF families, with all expressed genes as background. Enrichment factor was calculated for each term, and FDR < 0.05 was set as the significance threshold. Enrichment results were visualized via the bioinformatics online platform ( https://www.bioinformatics.com.cn/ ) [ 36 ]. To explore potential regulatory interactions, Pearson correlation analysis was performed between expression levels of differentially expressed HSFs and key root development regulators (WOX and ARF). |r| > 0.8 and p < 0.05 were considered strong and significant correlations. Correlation results were visualized as a heatmap using the Metware Cloud platform ( https://cloud.metware.cn ) [ 40 ]. Analysis of members of the HSF TF family Based on the expression profiles of differentially expressed HSF genes, correlation coefficients among family members were calculated and visualized as a circular graph. For phylogenetic analysis, HSF protein sequences from Arabidopsis thaliana , Solanum lycopersicum , and Oryza sativa L. were downloaded from National Center for Biotechnology Information (NCBI - https://www.ncbi.nlm.nih.gov/ ). Selected tobacco HSF sequences were aligned with these homologs using MEGA11 [ 39 ], and a Neighbor-Joining (NJ) phylogenetic tree was constructed and visualized with iTOL ( https://itol.embl.de/ ) [ 22 ]. Gene structure information for the 17 tobacco HSFs was extracted from the genome annotation file. Motif analysis was performed using the Multiple Em for Motif Elicitation (MEME) online tool ( http://memesuite.org/ ) [ 3 ] with the number of motifs set to 20 (other parameters default), and results were visualized using TBtools. Multispecies collinearity analysis and prediction of protein three-dimensional (3D) structure and cis-acting elements Genomic data of N. sylvestris , N. tomentosiformis , N. attenuata , N. benthamiana , A. thaliana , N. tabacum , and O. sativa L. were downloaded from NCBI. Collinearity analysis between the 17 tobacco HSFs and each species was performed using TBtools. For protein structure prediction, each HSF protein sequence was submitted to the SWISS-MODEL online platform ( https://swissmodel.expasy.org/ ) [ 43 ]; the template with the highest sequence identity was selected as the predicted 3D model. Promoter regions (2000 bp upstream of the transcription start site) of the 17 HSFs were extracted using TBtools. Cis -regulatory elements were predicted via PlantCARE ( https://bioinformatics.psb.ugent.be/webtools/plantcare/html/ ) [ 24 ], and results were visualized using TBtools. TF binding site prediction and protein-protein network construction HSF TF binding sites were predicted using the PlantRegMap database ([ 38 ]; https://plantregmap.gao-lab.org/ ). Target gene promoter sequences were input with a significance threshold of p ≤ 1 × 10⁻⁴. Based on the predicted TF-target regulatory relationships, a visual network was constructed using Cytoscape (v3.10.4) [ 33 ]. A protein-protein interaction (PPI) network was generated via the STRING database ( https://cn.string-db.org/ ). Amino acid sequences of target proteins were submitted with species specified; interactions were filtered at a confidence score ≥ 0.7. The resulting network data were downloaded and imported into Cytoscape for visualization and topological analysis. Expression profiles and genetic variation analysis The gene expression profile data used in this study were derived from previously published research, aiming to reveal the expression patterns of differentially expressed HSFs in hairy roots under abiotic stress and exogenous hormone regulation. Specifically, the data include transcriptome data from tobacco plants treated with exogenous salicylic acid (SA) at 0 h, 1 h, 4 h, and 12 h time points [ 41 ], as well as transcriptome data from tobacco plants subjected to drought stress and exogenous melatonin treatment [ 5 ]. Results Quality evaluation of RNA sequencing data and systematic screening of DEGs Transcriptome sequencing generated an average of 45.6 million clean reads per sample after quality control, achieving an average mapping rate of 91.6% (Supplementary Table S1). The Q20 values of all samples exceeded 94.7%, and the Q30 values exceeded 90.7% (Supplementary Table S1). Additionally, 68.9% of the transcripts were longer than 1000 bp (Figure S1). Collectively, the findings confirm the high-quality sequencing data of all samples. Furthermore, functional annotations revealed 89,483 expressed genes, including 81,307 known genes and 8176 novel genes. Additionally, 182,235 expressed transcripts were detected, including 135,537 known transcripts and 46,698 new transcripts. Then principal component analysis (PCA) demonstrated clear segregation between samples from different groups, whereas samples within the same group exhibited tight clustering (Fig. 1 A), highlighting significant inter-group divergence and high intra-group consistency. These results were further supported by the sample correlation clustering heatmap (Fig. 1 B). Fig. 1. Open in a new tab Clustering and differential gene expression analysis of transcriptome data. A PCA of different transcriptome samples. B Hierarchical clustering analysis of transcriptome samples. C Comparative analysis of DEG counts between hairy roots and normal plant roots/stems/leaves. D - E Venn diagrams of upregulated ( D ) and downregulated ( E ) genes for different comparison groups A comprehensive analysis of gene expression differences was conducted among various plant tissues. Compared to normal roots (G), stems (J), and leaves (Y), hairy roots exhibited 21,096 (11,151 upregulated and 9,945 downregulated), 21,787 (10,380 upregulated and 11,407 downregulated), and 23,997 (12,719 upregulated and 11,278 downregulated) DEGs, respectively (Fig. 1 C). To further elucidate the overlap and uniqueness of these DEGs of different tissue comparisons, a Venn diagram analysis was performed. This analysis revealed that 8860 DEGs (5670 upregulated and 3190 downregulated) were commonly differentially expressed among all the tissues examined (Fig. 1 D-E). Additionally, volcano plots were generated to visualize the DEGs between tobacco hairy roots and other tissues, including leaves, stems, and roots (Figure S2). These plots clearly indicated significant upregulation and downregulation of the DEGs, highlighting the distinct gene expression profiles of different tissue types. GO enrichment analysis highlights key molecular features of DEGs Enrichment analysis of all DEGs based on the GO database revealed that these DEGs were predominantly and significantly enriched in biological processes related to carbohydrate metabolic pathways (such as sucrose biosynthesis and fructose-1,6-bisphosphate metabolism), as well as the metabolism and degradation of amino sugars and chitin (including chitinase activity and binding functions) (Fig. 2 A-C). At the molecular function level, the DEGs were significantly enriched in photosystem-related components and thylakoid membrane-associated genes, suggesting that hairy roots exhibit distinct mechanisms of light energy capture and conversion compared to conventional plant roots, stems, and leaves (Fig. 2 A-C). Further analysis demonstrated significant enrichment of these DEGs across multiple functional categories: stress response systems (encompassing glutathione metabolism and reactive oxygen species (ROS) regulation), protein homeostasis maintenance (including proteasome core complex assembly, protein-folding chaperones, and autobinding functions), transmembrane signal transduction mechanisms (particularly receptor kinase activity and signal receptor functionality), and flavin adenine dinucleotide (FAD)-binding proteins (Fig. 2 A-C). This coordinais and isoquinoline alkaloid biosynthesis demonstrates hairy roots will be potential to serve as biological synthesis factory for specialized metabolites. Fig. 2. Open in a new tab Bar plot of GO enrichment analysis of all DEGs. A Bar plot of GO enrichment analysis of DEGs between hairy roots and normal plant roots. B Bar plot of GO enrichment analysis of DEGs between hairy roots and normal plant stems. C Bar plot of GO enrichment analysis of DEGs between hairy roots and normal plant leaves Identification of distinct metabolic and signaling pathways through KEGG analysis of DEGs The KEGG enrichment analysis of all DEGs revealed that these DEGs exhibited significant enrichment characteristics in signaling processes related to energy metabolism, stress response, and secondary metabolism. Then the significant enrichment of energy metabolism pathways (glycolysis/gluconeogenesis, the pentose phosphate pathway, the TCA cycle, and pyruvate metabolism) suggests a potential remodeling of energy metabolic flux in hairy roots (Fig. 3 A-C). Furthermore, as a unique in vitro culture system, hairy roots exhibit marked activation of photosynthesis-related pathways (including photosynthetic antenna proteins and carbon fixation) (Fig. 3 A-C), suggesting the potential evolution of an atypical light-harvesting and energy conversion system in this tissue. These finding is particularly remarkable given that root tissues conventionally lack significant photosynthetic capacity. Additionally, the enrichment of secondary metabolic pathways such as phenylpropanoid biosynthesis and isoquinoline alkaloid biosynthesis demonstrates hairy roots has the potential to serve as biofactories for specialized metabolites (Fig. 3 A-C). Based on these enrichment analysis results, we have elucidated the molecular basis for the unique characteristics of hairy roots: they maintain certain root-specific traits (such as sugar metabolism), acquire light-responsive potential typically associated with stem/leaf tissues, and evolve distinctive defense and biosynthetic capabilities. Fig. 3. Open in a new tab KEGG enrichment analysis of all DEGs. KEGG enrichment analysis results of DEGs between hairy roots and normal plant roots ( A ), hairy roots and normal plant stems ( B ), hairy roots and normal plant leaves ( C ) WGCNA and functional module identification of differential genes in hairy root systems Prior to WGCNA analysis, 8860 genes were screened, and those with low expression variability (standard deviation ≤ 0.5) were removed, resulting in 8554 genes retained for subsequent analysis. Based on the predetermined power value, a weighted gene co-expression network was constructed, and these genes were grouped into three modules (midnightblue and pink, the grey module is considered non-informative and thus excluded from further interpretation) (Fig. 4 A). The midnightblue module comprises 3033 genes and demonstrates a significantly positive correlation with the phenotype ( r = 0.73). In contrast, the pink module, containing 5464 genes, exhibits a significantly negative correlation with the phenotype ( r = − 0.87) (Fig. 4 B). Fig. 4. Open in a new tab Analysis of WGCNA Results. A In the hierarchical clustering dendrogram, the upper tree diagram clusters genes based on weighted correlation coefficients. The lower part illustrates the process of module identification and merging: the first row (Dynamic Tree Cut) shows the initially identified gene modules, where each color represents a preliminary module; the second row (merged dynamic) displays the final gene modules obtained by merging highly correlated preliminary modules. B Gene co-expression module–trait correlation heatmap analyzes the relationship between module eigengenes and target samples using Pearson correlation coefficients. Orange indicates a significant positive correlation ( r > 0), while blue indicates a significant negative correlation ( r < 0). Values in parentheses represent the corresponding p -values ( p < 0.05 considered significant). C Co-expression network of the top 50 hub genes in the pink module. The lines in the figure indicate the degree of connectivity between genes, with larger circles representing higher connectivity. The more connections a gene has with surrounding nodes, the more central and important its position within the network. D Co-expression network of the top 50 hub genes in the midnightblue module. The lines in the figure indicate the degree of connectivity between genes, with larger circles representing higher connectivity. The more connections a gene has with surrounding nodes, the more central and important its position within the network Co-expression networks were constructed using the top 50 genes from two modules highly correlated with specific phenotypes, leading to the identification of hub genes. Five hub genes were identified in the pink module, including Nitab4.5_0000073g0500 (UDP-glucosyltransferase), Nitab4.5_0003969g0060 (heat shock protein), Nitab4.5_0006998g0030 (cation/H⁺ exchanger), Nitab4.5_0000745g0040 (MYB TF), and Nitab4.5_0003098g0040 (cysteine/histidine-rich zinc finger protein) (Fig. 4 C). In the midnightblue module, two hub genes were identified: Nitab4.5_0001133g0020 (BURP domain-containing protein) and Nitab4.5_0016580g0020 (receptor-like kinase) (Fig. 4 D). A total of 33 TFs were predicted to bind to the promoter regions of seven core genes, with the C2H2 (seven), Dof (six), and HSF (five) families being the most abundant. Nitab4.5_0001003g0080 (AP2/ERF) and Nitab4.5_0003182g0070 (MIKC_MADS) exhibited high connectivity in the regulatory network (Figure S3). Transcriptional regulator profiling of common DEGs in hairy roots and root-stem-leaf systems Through further analysis of transcriptome data from hairy roots and normal root, stem, and leaf tissues, we systematically identified differentially expressed TFs (DEG-TFs) and performed gene family enrichment analysis. Then we identified 17 differentially expressed HSFs (15 of them were specifically upregulated in hairy roots) (Fig. 5 A), with this family showing the highest enrichment among all differentially expressed TFs (Fig. 5 B). GO enrichment analysis revealed that these HSFs were significantly enriched in DNA-binding transcriptional regulation activity (molecular function) and intracellular macromolecule synthesis/metabolism-related terms (biological process) (Fig. 5 C). Previous study identified auxin response factor ( ARF ) and WUSCHEL related homeobox ( WOX ) genes are two key gene families that co-regulate root organogenesis [ 45 , 46 ]. Thus, we performed correlation analysis between the HSF and WOX/ARF family members involved in root development. Among these, 12 HSFs showed significant co-expression with two WOXs (including Nitab4.5_0011394g0010 and Nitab4.5_0000495g0100 ), and 14 HSFs showed significant co-expression with 11 ARFs related to root development ( p < 0.05) (Figs. 5 D-E). Of these WOX and ARF genes associated with HSF expression, four different expression ARFs ( Nitab4.5_0003923g0040 , Nitab4.5_0000899g0130 , Nitab4.5_0000304g0070 , Nitab4.5_0000798g0140 ) were clustered into the midnightblue module as the two downregulated HSFs , and three ARFs ( Nitab4.5_0009330g0010 , Nitab4.5_0000476g0080 , Nitab4.5_0004657g0030 ) were clustered into the pink module as the 15 upregulated HSFs. These findings suggest that these HSF TFs may participate in the transcriptional regulation of root development through direct modulation or modular synergistic effects. Fig. 5. Open in a new tab Analysis of differentially expressed TFs. A Classification and quantification of all differentially expressed TFs. B Gene set enrichment analysis of all differentially expressed TFs. C GO enrichment analysis of the differentially expressed HSF gene family. D - E Correlation analysis between differentially expressed HSFs and key root development factors: auxin response factor ( ARF ) ( D ) and WUSCHEL-related homeobox ( WOX ) ( E ) Functional prediction and evolutionary analysis of differentially expressed HSFs As key transcriptional regulators, HSFs are essential for maintaining cellular homeostasis under abiotic stresses [ 48 ], and we conducted structural and evolutionary analyses of 17 HSFs identified from transcriptome data to investigate their potential roles in hairy root formation. Through correlation analysis of these HSFs , it was found that two genes downregulated in hairy roots showed significant negative correlations with 15 genes upregulated in hairy roots (Fig. 6 A). This result is consistent with the expression trends of these genes in the roots, stems, and leaves of normal plants and hairy roots. Evolutionary analysis classified these 17 HSFs into five subgroups, with the two downregulated genes independently distributed in the first and fifth subgroups (Fig. 6 B). Conservative motif analysis revealed that all these HSFs contain the two key conserved domains, motif 1 and motif 3, then with the exception of Nitab4.5_0006360g0010 , all have a simplified two-exon structure (Fig. 6 C). Collinearity analysis further revealed the evolutionary trajectory of tobacco HSFs: N. tabacum shares 10 and 11 homologous gene pairs with its ancestral species N. sylvestris and N. tomentosiformis , respectively, while only two pairs are retained with the distantly related species N. attenuata (Fig. 6 D). This discrepancy not only confirms the dual ancestral origin of cultivated tobacco but also suggests functional divergence of the HSF family during speciation. Interspecies collinearity analysis identified four and five homologous gene pairs between tobacco HSFs and those in Arabidopsis and rice, respectively (Fig. 6 D). Fig. 6. Open in a new tab Analysis of 17 differential HSF gene families. A Correlation among the 17 HSFs . Red lines between gene IDs indicate a positive correlation, while purple lines indicate a negative correlation. B Phylogenetic tree analysis of HSF gene families in Arabidopsis , rice, tomato, and tobacco. C Motif and gene structure analysis. D Synteny analysis of HSFs in tobacco with their ancestral species, closely related species, as well as Arabidopsis and rice. E Analysis of cis -acting elements in the promoter regions. F 3D protein structure analysis Analysis of promoter cis -acting elements revealed that the promoter regions of these HSFs are enriched with various regulatory elements related to stress response, hormone regulation, and metabolism (Fig. 6 E), which may be associated with their pleiotropic functions in hairy root development. Further prediction revealed that 76 TFs could bind to the promoter regions of 17 HSFs , with a higher number of ERF (19), MYB (12), and HSF (seven) families. Nitab4.5_0001817g0020 (ERF), Nitab4.5_0001270g0190 (C2H2), and Nitab4.5_0003182g0070 (MIKC_MADS) showed high connectivity in the regulatory network (Figure S4). 3D structure prediction showed that HSF proteins within the same subfamily exhibit highly similar structures, consistent with the gene structure analysis (Fig. 6 F). Given that the occurrence of hairy roots is highly dependent on the dynamic balance of endogenous hormone levels, this study further focuses on the hormone signaling genes that show significant responses in hairy roots and analyzes their correlation with HSFs . The results indicate that the formation of hairy roots involves coordinated reprogramming of multiple hormone signaling pathways: key response factors (auxin ARF , cytokinin ARR , ethylene EIN3 , gibberellin PIF4 , and abscisic acid (ABA) PYL ) were significantly downregulated, while multiple metabolic and signaling components (such as auxin GH , cytokinin HK3 , ethylene EIN2 , and SA pathway genes) were generally upregulated (Fig. 7 A). Further analysis revealed that these upregulated hormone-related genes showed significant positive correlations with all 15 HSFs , suggesting that HSFs may participate in integrating hormone signals through extensive interactions, thereby cooperatively regulating hairy root development (Fig. 7 B). Fig. 7. Open in a new tab Analysis of differentially expressed hormone-related genes in hairy roots. A Heatmap of expression levels of DEGs involved in hormone pathways in hairy roots. The annotations on the right include the gene ID, corresponding gene annotation, and the associated hormone signaling. The color gradient from blue to red represents expression abundance from low to high. B Heatmap of correlations between differentially expressed hormone pathway genes and 17 HSFs in hairy roots. Pearson correlation coefficients were used for the analysis (* indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001). The color gradient from green to white indicates the strength of correlation, with green representing negative correlation and white representing positive correlation Screening key hairy roots regulating candidate genes through expression profiling and whole-genome resequencing Through a systematic analysis of the expression patterns and genetic variation characteristics of 17 HSFs in tobacco hairy roots, we identified two key categories of HSFs . The first category includes two genes that are specifically expressed in hairy roots but show almost no expression in normal plant tissues ( Nitab4.5_0002158g0140 and Nitab4.5_0002782g0100 ) (Fig. 8 ), suggesting their potential direct involvement in the initiation or maintenance of hairy roots. The second category comprises nine genes ( Nitab4.5_0000082g0150 , Nitab4.5_0000459g0060 , Nitab4.5_0001519g0010 , Nitab4.5_0002328g0020 , Nitab4.5_0004411g0050 , Nitab4.5_0004441g0010 , Nitab4.5_0005947g0030 , Nitab4.5_0006360g0010 , and Nitab4.5_0006852g0070 ) that are not only highly expressed in hairy roots but also responsive to various treatments such as SA, drought and melatonin (Fig. 8 ), demonstrating their broad responsiveness to both endogenous and exogenous signals. Furthermore, by integrating whole-genome resequencing data, we discovered abundant structural variations in the intronic regions of these key genes. Specifically, Nitab4.5_0002328g0020 contains 14 single-nucleotide polymorphisms (SNPs) and four InDels, including an 11 bp deletion, while Nitab4.5_0002904g0100 harbors one SNP and seven InDels, including a 27 bp insertion (Supplementary table S2). These variations may regulate gene expression by affecting splicing efficiency, transcriptional stability, or the activity of regulatory elements, ultimately contributing to the formation of the hairy root phenotype. Fig. 8. Open in a new tab Analysis of the expression patterns of 17 HSF genes. The symbols R, G, S, and Y represent samples from different tissues: R denotes hairy roots, while G, S, and Y denote the root, stem, and leaf of normal plants, respectively. The abbreviations SA0, SA1, SA4, and SA12 means the time points of 0, 1, 4, and 12 h after SA treatment, respectively. The group of CK, MEL, D, and MEL_D, respectively represent the control group (CK_1, CK_2, and CK_3), melatonin (MEL_1, MEL_2, MEL_3), drought (D_1, D_2, D_3) and the combination of drought and melatonin (MEL_D_1, MEL_D_2, and MEL_D_3) treatment Topological analysis of the HSF-mediated regulatory network for hairy root formation Through intersection analysis of key candidate genes, we found that during hairy root induction, Nitab4.5_0010472g0010 can regulate downstream core genes while also feedback-regulating the HSF family itself, whereas Nitab4.5_0005947g0030 acts as a key upstream regulator of the HSF family (Figure S5). Further construction of a protein-protein interaction network, including core genes, the HSF family, and their significantly co-expressed genes (such as ARF and hormone-related genes), revealed that Nitab4.5_0001098g0070 ( PYL8 ), Nitab4.5_0003335g0010 ( OST1 ), Nitab4.5_0003710g0010 ( PYL6 ), Nitab4.5_0000786g0150 (protein kinase), and Nitab4.5_0013277g0010 ( ABI2 ) exhibited the highest connectivity (Figure S6). These genes encode proteins primarily involved in the ABA signaling pathway and kinase cascades, indicating that ABA signaling plays a central role in regulating hairy root formation. This result is consistent with previous findings showing significant correlations between HSFs and hormone-related genes, suggesting that HSFs may participate in the regulation of hairy root formation by mediating the ABA signaling pathway. Discussions Plant hairy roots are characterized by strong genetic stability, rapid growth rate, and short production cycle, making them an important bioreactor and effective tool for the large-scale production of secondary metabolites such as flavonoids [ 44 ]. In this study, DEGs identified through comparative transcriptomics between hairy roots and normal plant tissues were significantly enriched in energy metabolism, stress response, and secondary metabolism (such as the pentose phosphate pathway, the TCA cycle, etc.), indicating that while retaining some root functions, hairy roots undergo extensive metabolic reprogramming and structural adaptive changes. Previous study investigated the influence of the rolB gene on the ARF TF gene family. They found that the rolB gene could promote the development of tobacco hairy roots by selectively inducing NtARF7 and NtARF19 . We also identified a series of ARFs that were specifically highly expressed in hairy roots. Compared with previous findings on the auxin-signal-responsive rol genes [ 4 ], our study further reveals extensive transcriptional reprogramming in hairy roots. This reprogramming involves multiple pathways, including metabolism, secondary biosynthesis, and stress response, with the HSF family identified as playing a central regulatory role. We further identified a series of ARFs that are specifically highly expressed in hairy roots. These ARFs showed a highly significant correlation with the differentially expressed HSF TFs screened in this study. These results further support the central role of auxin signaling in hairy root development. Compared to the studies by Qin et al., [ 28 ] and Strotmann et al., [ 34 ] on the regulation of root development and metabolism by WOX and MYB/bHLH TFs in model plants, this study also identified significant differential expression of multiple TF family members, such as MYB, bHLH, WRKY, NAC, and HSF, in hairy roots. Furthermore, this research conducted an in-depth analysis of the HSF TFs, which showed the most significant enrichment in the detection results, and screened several key HSF candidate genes involved in the formation of hairy roots. HSFs are a class of regulatory proteins that play a central role in plant stress responses and are widely involved in responses to various biotic and abiotic stresses [ 13 ]. Studies have shown that the upstream open reading frame (uORF) in the 5’ untranslated region (5’UTR) of the HsfA1a gene regulates the expression of WOX11 at the translational level, thereby promoting the development of crown roots in rice [ 47 ]. Additionally, overexpression of HSFB2b significantly enhances the elongation of soybean hairy roots and their tolerance to salt stress [ 2 ]. In this study, a systematic analysis of differentially expressed TFs between tobacco hairy roots and normal root, stem, and leaf tissues revealed that the HSF gene family was the most significantly enriched in hairy roots, suggesting its potential key role in the development or functional regulation of hairy roots. The formation of hairy roots is characterized by the coordinated reprogramming of multiple hormone signaling pathways: key response factors (such as ARF , ARR , EIN3 , PIF4 , and PYL ) are significantly downregulated, suggesting the suppression of classical hormone signal transduction pathways, while hormone metabolism and signaling components (such as GH3 , HK3 , EIN2 , and SA pathway genes) are generally upregulated, indicating the specific activation of hormone synthesis, modification, and signal feedback loops. Further analysis revealed that all upregulated hormone-related genes show significant positive correlations with all 15 HSFs , suggesting that the HSF family may act as an integrative node in the hormone signaling network, coordinating cross-talk among multiple pathways through extensive interactions and collectively driving transcriptional reprogramming related to hairy root development. To investigate the response mechanisms of tobacco roots to the biotic stress induced by A. rhizogenes infection, this study simulated three potential physiological responses-defense signal activation, osmotic/oxidative stress, and endogenous adaptive regulation by applying treatments with drought [ 14 ], and exogenous melatonin [ 17 ], SA [ 37 ], respectively. In this study, we successfully identified four HSF s specifically expressed in hairy roots. Additionally, we found three genes that not only exhibit high basal expression levels in this system but are also significantly induced by various treatments, including SA, drought, and melatonin. These findings indicate that these genes may possess dual functions in multi-stress response and developmental regulation. Notably, abundant SNP and large-fragment InDel variations were identified in the intronic regions of genes Nitab4.5_0002328g0020 and Nitab4.5_0002904g0100 . These structural variations may affect transcriptional stability, mRNA splicing efficiency, or protein function, thereby participating in the molecular mechanisms regulating hairy root morphogenesis. Based on the current results, future research could focus on the four HSFs specifically expressed in hairy roots identified through screening. Functional validation using gene editing and overexpression techniques can be employed to clarify their specific roles in hairy root morphogenesis and their upstream and downstream regulatory relationships. Secondly, the cross-regulatory mechanisms of the three HSFs involved in both development and multi-stress responses across different signaling pathways provide a theoretical foundation for studying plant development and stress resistance mechanisms. The abundant intronic variations (such as SNPs/InDels) present in these genes may influence their splicing regulation and expression efficiency, offering potential applications for molecular marker-assisted selective breeding. Finally, applying these candidate genes to optimize root architecture and enhance stress resistance in major crops (such as soybeans and rice) holds promise for developing new germplasms with more robust root systems and stronger stress tolerance. This would provide a molecular basis and genetic resources for sustainable agriculture, while also offering theoretical support and new targets for future efforts in synthetic biology aimed at genetically engineering hairy roots to improve their growth capacity and target compound production. Conclusion In this study, comparative transcriptomic analysis of tobacco hairy roots and normal root, stem, and leaf tissue samples revealed that DEGs were significantly enriched in signaling pathways related to energy metabolism, stress response, and secondary metabolism, such as glycolysis/gluconeogenesis and the TCA cycle in energy metabolism. Enrichment analysis of differentially expressed TFs indicated that the HSF gene family was the most significantly enriched. Further phylogenetic classification and expression profile variation analysis of 17 differentially expressed HSFs led to the successful identification of two hairy root-specific expressed genes, as well as nine genes that not only exhibited constitutive high expression in this system but also responded to various induction treatments such as SA, melatonin and drought treatments. Furthermore, we found that HSFs exhibited highly significant correlations with hormone pathway genes. We also constructed a protein-protein interaction network that includes hub genes, the HSF family, and their significantly co-expressed genes, such as ARF and hormone-related genes. This study provides a theoretical basis and potential targets for developing new germplasms with more robust root systems and enhanced stress resistance, as well as for directed improvement of hairy root traits using synthetic biology approaches. Supplementary Information 12870_2026_8531_MOESM1_ESM.png (361.1KB, png) Supplementary Material 1. Supplementary Figure 1. Distribution of transcript length. 12870_2026_8531_MOESM2_ESM.png (1,021.1KB, png) Supplementary Material 2. Supplementary Figure 2. Volcano plot of DEGs from comparative analysis between hairy roots and roots/stems/leaves of normal plants. 12870_2026_8531_MOESM3_ESM.png (1MB, png) Supplementary Material 3. Supplementary Figure 3. Regulatory network of upstream TFs for WGCNA hub genes. In the figure, the red, purple, and green circles represent the WGCNA hub genes, the two TFs with the highest predicted connectivity, and other predicted TFs, respectively. The size of the circles corresponds to the level of connectivity. 12870_2026_8531_MOESM4_ESM.png (2MB, png) Supplementary Material 4. Supplementary Figure 4. Regulatory network of upstream TFs for HSFs . In the figure, the purple, green, and pink circles represent the HSFs , the three TFs with the highest predicted connectivity, and other predicted TFs, respectively. The size of the circles corresponds to the level of connectivity. 12870_2026_8531_MOESM5_ESM.png (30.4KB, png) Supplementary Material 5. Supplementary Figure 5. Venn diagram illustrating the overlap among the HSF gene set, the hub gene set, the predicted upstream TF set for HSFs , and the predicted upstream TF set for hub genes. 12870_2026_8531_MOESM6_ESM.png (790.5KB, png) Supplementary Material 6. Supplementary Figure 6. Protein interaction network of HSF family members, co-expressed core genes, and ARF TFs and hormone-responsive genes screened through co-expression associations. 12870_2026_8531_MOESM7_ESM.xlsx (12.2KB, xlsx) Supplementary Material 7. Supplementary table 1. Quality control of transcriptomic data. Supplementary table 2. Detailed information of different expresses HSF gene. Authors’ contributions Xiaozong Wu and Chaonan Shi designed the study. Zhitao Qi, Yaolan Zhuang and Haoyu Yan were responsible for the planting of tobacco seedlings and the induction of tobacco hairy roots. Lei Wang, Jifan He and Peilin Li participated in the analysis of transcriptomic data and the visualization of figures. Lei Wang and Zhiwen Zhu performed the gene family analysis. Yixuan Xue conducted the correlation analysis among genes and the protein-protein interaction network analysis. Zhitao Qi and Meng Li contributed to the drafting of the initial manuscript, while Xiaozong Wu, Chaonan Shi and Guiliang Tang were involved in revising and finalizing the article. Funding This research was funded by the Henan Province Science and Technology Research Project (Project No.: 252102110270 and 242102110334), Henan Provincial Natural Science Foundation (Project No.: 252300420221), Doctoral fund project of Zhengzhou University of Light Industry, China (Project No.: 2024BSJJ011). Data availability The authors declare that all data involved in this study are publicly available without restrictions. The transcriptome sequencing data generated in this study have been deposited in the NCBI Sequence Read Archive database under the accession number PRJNA1328999. Declarations Ethics approval and consent to participate Our study did not involve any human or animal subjects, material, or data. We declare that the plant material in the experiment was collected and studied by relevant institutional, national, and international guidelines and legislation. Consent for publication Not applicable. Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Contributor Information Meng Li, Email: [email protected]. Chaonan Shi, Email: [email protected]. References 1. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 2. Bian XH, Li W, Niu CF, Wei W, Hu Y, Han JQ, et al. A class B heat shock factor selected for during soybean domestication contributes to salt tolerance by promoting flavonoid biosynthesis. New Phytol. 2020;225(1):268–83. [ DOI ] [ PubMed ] [ Google Scholar ] 3. Bailey TL, Johnson J, Grant CE, Noble WS. The MEME Suite. Nucleic Acids Res. 2015;43(W1):W39–49. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 4. Bose R, Sengupta M, Basu D, Jha S. The rolB-transgenic Nicotiana tabacum plants exhibit upregulated ARF7 and ARF19 gene expression. Plant Direct. 2022;6(6):e414. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 5. Chen Z, Jia W, Li S, Xu J, Xu Z. Enhancement of Nicotiana tabacum Resistance Against Dehydration-Induced Leaf Senescence via Metabolite/Phytohormone-Gene Regulatory Networks Modulated by Melatonin. Front Plant Sci. 2021;12:686062. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 6. Cui G, Li Y, Zheng L, Smith C, Bevan MW, Li Y. The peptidase DA1 cleaves and destabilizes WUSCHEL to control shoot apical meristem size. Nat Commun. 2024;15(1):4627. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 7. Chen C, Wu Y, Li J, Wang X, Zeng Z, Xu J, et al. TBtools-II: A one for all, all for one bioinformatics platform for biological big-data mining. Mol Plant. 2023;16(11):1733–42. [ DOI ] [ PubMed ] [ Google Scholar ] 8. Cantalapiedra CP, Hernández-Plaza A, Letunic I, Bork P, Huerta-Cepas J. eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol Biol Evol. 2021;38(12):5825–9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 9. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 10. Cao H, Liu Z, Guo J, Jia Z, Shi Y, Kang K, et al. ZmNRT1.1B (ZmNPF6.6) determines nitrogen use efficiency via regulation of nitrate transport and signalling in maize. Plant Biotechnol J. 2024;22(2):316–29. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Dong Z, Hu G, Chen Q, Shemyakina EA, Chau G, Whipple CJ, et al. A regulatory network controlling developmental boundaries and meristem fates contributed to maize domestication. Nat Genet. 2024;56(11):2528–37. [ DOI ] [ PubMed ] [ Google Scholar ] 12. Edwards KD, Fernandez-Pozo N, Drake-Stowe K, Humphry M, Evans AD, Bombarely A, et al. A reference genome for Nicotiana tabacum enables map-based cloning of homeologous loci implicated in nitrogen utilization efficiency. BMC Genomics. 2017;18(1):448. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 13. Fang Y, Liao H, Wei Y, Yin J, Cha J, Liu X, et al. OsCDPK24 and OsCDPK28 phosphorylate heat shock factor OsHSFA4d to orchestrate abiotic and biotic stress responses in rice. Nat Commun. 2025;16(1):6485. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 14. Gorgues L, Li X, Maurel C, Martinière A, Nacry P. Root osmotic sensing from local perception to systemic responses. Stress Biol. 2022;2(1):36. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 15. Hwang C, Yan S, Choe Y, Yun C, Xu S, Im M, et al. Efficient hairy root induction system of Astragalus membranaceus and significant enhancement of astragalosides via overexpressing AmUGT15 . Plant Cell Rep. 2024;43(12):285. [ DOI ] [ PubMed ] [ Google Scholar ] 16. Hernández-Plaza A, Szklarczyk D, Botas J, Cantalapiedra CP, Giner-Lamia J, Mende DR, et al. eggNOG 6.0: enabling comparative genomics across 12 535 organisms. Nucleic Acids Res. 2023;51(D1):D389–94. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 17. Jiang L, Yuan Z, Yan W, Tang P, Yuan P, Zheng P, et al. Transcriptomic and metabolomic analyses unveil TaASMT3-mediated wheat resistance against stripe rust by promoting melatonin biosynthesis. Plant J. 2025;122(2):e70182. [ DOI ] [ PubMed ] [ Google Scholar ] 18. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 19. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 20. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 21. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 22. Letunic I, Bork P. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 2021;49(W1):W293–6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 23. Liu X, Singh SK, Patra B, Liu Y, Wang B, Wang J, et al. Protein phosphatase NtPP2C2b and MAP kinase NtMPK4 act in concert to modulate nicotine biosynthesis. J Exp Bot. 2021;72(5):1661–76. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 24. Lescot M, Déhais P, Thijs G, Marchal K, Moreau Y, Van de Peer Y, et al. PlantCARE, a database of plant cis -acting regulatory elements and a portal to tools for in silico analysis of promoter sequences. Nucleic Acids Res. 2002;30(1):325–7. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Ma S, Chen Q, Zheng Y, Ren T, He R, Cheng L, et al. A tale for two roles: Root-secreted methyl ferulate inhibits P. nicotianae and enriches the rhizosphere Bacillus against black shank disease in tobacco. Microbiome. 2025;13(1):33. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 27. Pan L, Chen Y, Ren Z, Khojely DM, Wang S, Li Y, et al. Using WGCNA and transcriptome profiling to identify hub genes for salt stress tolerance in germinating soybean seeds. Front Plant Sci. 2025;16:1569565. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 28. Qin J, Zhao C, Wang S, Gao N, Wang X, Na X, et al. PIF4-PAP1 interaction affects MYB-bHLH-WD40 complex formation and anthocyanin accumulation in Arabidopsis. J Plant Physiol. 2022;268:153558. [ DOI ] [ PubMed ] [ Google Scholar ] 29. Ren Y, Yu G, Shi C, Liu L, Guo Q, Han C, et al. Majorbio Cloud: A one-stop, comprehensive bioinformatic platform for multiomics analyses. Imeta. 2022;1(2):e12. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 30. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 31. Su L, Liu S, Liu X, Zhang B, Li M, Zeng L, et al. Transcriptome profiling reveals histone deacetylase 1 gene overexpression improves flavonoid, isoflavonoid, and phenylpropanoid metabolism in Arachis hypogaea hairy roots. PeerJ. 2021;9:e10976. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 32. Su X, Zhang X, Luo J, Wang Y, Feng B, Yang Y, et al. The IAA7-ARF7-ARF19 auxin signaling module plays diverse roles in Arabidopsis growth and development. Planta. 2025;262(1):12. [ DOI ] [ PubMed ] [ Google Scholar ] 33. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 34. Strotmann VI, García-Gómez ML, Stahl Y. Root stem cell homeostasis in Arabidopsis involves cell-type specific transcription factor complexes. EMBO Rep. 2025;26(9):2323–46. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 35. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 36. Tang D, Chen M, Huang X, Zhang G, Zeng L, Zhang G, et al. SRplot: A free online platform for data visualization and graphing. PLoS ONE. 2023;18(11):e0294236. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 37. Tang B, Lu J, Leontovyčová H, Hoffmann G, Rowe JH, O’Donnell SF, et al. AM. SALICYLIC ACID SENSOR1 reveals the propagation of an SA hormone surge during plant pathogen advance. Science. 2025;390(6769):188–94. [ DOI ] [ PubMed ] [ Google Scholar ] 38. Tian F, Yang DC, Meng YQ, Jin J, Gao G. PlantRegMap: charting functional regulatory maps in plants. Nucleic Acids Res. 2020;48(D1):D1104–13. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 39. Tamura K, Stecher G, Kumar S. MEGA11: Molecular Evolutionary Genetics Analysis Version 11. Mol Biol Evol. 2021;38(7):3022–7. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 40. Wang P, Wu X, Shi Z, Tao S, Liu Z, Qi K, et al. A large-scale proteogenomic atlas of pear. Mol Plant. 2023;16(3):599–615. [ DOI ] [ PubMed ] [ Google Scholar ] 41. Wu XZ, Guo WW, Zhu ZW, Li M, Xu JB, Zhu RJ, et al. Studies on flavonoids biosynthesis genes expression induced by salicylic acid in Nicotiana tabacum L. J light Ind. 2025;40(2):80–9. [ Google Scholar ] 42. Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8. [ DOI ] [ PubMed ] [ Google Scholar ] 43. Waterhouse A, Bertoni M, Bienert S, Studer G, Tauriello G, Gumienny R, et al. SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 2018;46(W1):W296–303. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 44. Zhu YT, Zhu X, Wang LH, Wang YL, Liao CL, Zhao M, et al. Plant hairy roots: Induction, applications, limitations and prospects. Ind Crops Prod. 2024;219:119104. 45. Zhang T, Ge Y, Cai G, Pan X, Xu L. WOX-ARF modules initiate different types of roots. Cell Rep. 2023;42(8):112966. [ DOI ] [ PubMed ] [ Google Scholar ] 46. Zhang C, Guo X, Wang H, Dai X, Yan B, Wang S, et al. Induction and metabolomic analysis of hairy roots of Atractylodes lancea . Appl Microbiol Biotechnol. 2023;107(21):6655–70. [ DOI ] [ PubMed ] [ Google Scholar ] 47. Zhang T, Xiang Y, Ye M, Yuan M, Xu G, Zhou DX, et al. The uORF-HsfA1a-WOX11 module controls crown root development in rice. New Phytol. 2025;247(2):760–73. [ DOI ] [ PubMed ] [ Google Scholar ] 48. Zhang H, Zhu J, Gong Z, Zhu JK. Abiotic stress responses in plants. Nat Rev Genet. 2022;23(2):104–19. [ DOI ] [ PubMed ] [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials 12870_2026_8531_MOESM1_ESM.png (361.1KB, png) Supplementary Material 1. Supplementary Figure 1. Distribution of transcript length. 12870_2026_8531_MOESM2_ESM.png (1,021.1KB, png) Supplementary Material 2. Supplementary Figure 2. Volcano plot of DEGs from comparative analysis between hairy roots and roots/stems/leaves of normal plants. 12870_2026_8531_MOESM3_ESM.png (1MB, png) Supplementary Material 3. Supplementary Figure 3. Regulatory network of upstream TFs for WGCNA hub genes. In the figure, the red, purple, and green circles represent the WGCNA hub genes, the two TFs with the highest predicted connectivity, and other predicted TFs, respectively. The size of the circles corresponds to the level of connectivity. 12870_2026_8531_MOESM4_ESM.png (2MB, png) Supplementary Material 4. Supplementary Figure 4. Regulatory network of upstream TFs for HSFs . In the figure, the purple, green, and pink circles represent the HSFs , the three TFs with the highest predicted connectivity, and other predicted TFs, respectively. The size of the circles corresponds to the level of connectivity. 12870_2026_8531_MOESM5_ESM.png (30.4KB, png) Supplementary Material 5. Supplementary Figure 5. Venn diagram illustrating the overlap among the HSF gene set, the hub gene set, the predicted upstream TF set for HSFs , and the predicted upstream TF set for hub genes. 12870_2026_8531_MOESM6_ESM.png (790.5KB, png) Supplementary Material 6. Supplementary Figure 6. Protein interaction network of HSF family members, co-expressed core genes, and ARF TFs and hormone-responsive genes screened through co-expression associations. 12870_2026_8531_MOESM7_ESM.xlsx (12.2KB, xlsx) Supplementary Material 7. Supplementary table 1. Quality control of transcriptomic data. Supplementary table 2. Detailed information of different expresses HSF gene. Data Availability Statement The authors declare that all data involved in this study are publicly available without restrictions. The transcriptome sequencing data generated in this study have been deposited in the NCBI Sequence Read Archive database under the accession number PRJNA1328999. Articles from BMC Plant Biology are provided here courtesy of BMC ACTIONS View on publisher site PDF (9.8 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top