Transcriptomic and metabolomic analysis reveal the cold tolerance mechanism of common beans under cold stress
With the intensification of the global greenhouse effect, fluctuations in global climate have increased, leading to a higher frequency of extreme temperature events. Consequently, plants are experiencing increased exposure to both cold and heat stress . The impact of low temperatures on plants is becoming increasingly severe, resulting in hindered growth and development, reduced yield, and decreased quality. In extreme cases, it can lead to crop death, causing significant losses in agricultural production . Low temperature is a critical factor affecting the overall developmental stages of plant and crop phenology, often resulting in decreased crop yield and quality . According to statistics, cold stress has resulted in a 60% reduction in global yields of legume crops, such as chickpeas and soybeans, and a 70% decrease in mung beans . When subjected to low-temperature stress, plants undergo a series of physiological and biochemical reactions to mitigate the adverse effects of the environmental temperature. Overall, the physiological strategies of plants to cope with low-temperature stress primarily involve the cell membrane system, antioxidant system, osmotic regulation system, metabolic substances, and regulation of endogenous hormones . Plant metabolites result from gene expression, protein interactions, and various regulatory mechanisms, and they are more closely associated with the plant phenotype than mRNA transcripts and proteins. Therefore, plant metabolites are often analyzed to examine plant phenotypes and provide feedback on environmental stress, facilitating the discovery of specific patterns related to stress tolerance . Sugar acid metabolism and amino acid metabolism play significant roles in mediating plant resistance to low-temperature stress . Mu et al. found that low temperature induces the transcriptional expression of fructose-1,6-diphosphate aldolase (SiFBA5) in Saussurea invalucrata . They observed that overexpressing SiFBA5 tomatoes exhibited enhanced cold tolerance and photosynthetic efficiency. In addition, secondary metabolites play integral roles in physiological processes, including plant growth, development, and defense . Previous studies have shown that plant secondary metabolites respond to low-temperature stress and alleviate oxidative damage by scavenging reactive oxygen species (ROS) that accumulate due to low-temperature stress . Analysis of metabolites from the alpine plant under low-temperature conditions revealed a notable trend: as the temperature decreased, there was a decrease in amino acid accumulation within the leaves, accompanied by an increase in phenolic substances. This observation suggests that under low-temperature conditions, Saussurea involucrata may enhance its cold tolerance by augmenting the production of secondary metabolites, particularly phenolic substances . Meanwhile, research has demonstrated that exposure to a 15 ℃ treatment significantly augments potato stem diameter, root-to-shoot ratio, yield, and concentration of secondary metabolites, especially anthocyanin content, suggesting that appropriate low-temperature treatments can be advantageous for enhancing potato tuber pigmentation . Plant hormones serve as pivotal regulators not only in orchestrating active substances in plant physiological responses and governing plant growth, development, and differentiation but also in mediating stress responses under adverse environmental conditions . Previous studies have shown that plant hormones, acting as vital growth regulators, modulate physiological and biochemical traits under low-temperature stress . They alleviate oxidative damage by enhancing the accumulation of proline, antioxidants, secondary metabolites, and endogenous hormone content, thereby improving plant tolerance to low-temperature stress . Gao et al.‘s research suggests that brassinosteroid may function as an upstream signal of NO, inducing an increase in NO content within the plant body and subsequently inducing protein S-nitrosylation modification to alleviate damage to Chinese cabbage seedlings under low-temperature stress . Meanwhile, another study revealed that the exogenous application of methyl jasmonate (MeJA) could mitigate oxidative damage to Solanum lycopersicum under low-temperature stress by enhancing antioxidant enzyme activity and photosynthetic activity . In addition, transcription factors contribute significantly to plant cold tolerance. Research has shown that long hypocotyl 5 ( HY5 ) can directly regulate the transcription level of CBF or indirectly influence CBF expression through MYB15 , thereby precisely regulating tomato cold tolerance . In addition, functional enrichment analysis of the differentially expressed genes in important modules of Hordeum vulgare L. ‘s cold stress response revealed that these genes are involved in various key pathways related to cold tolerance in plants, such as the ABA signaling pathway, ROS signaling pathway, defense and protective proteins, and degrading proteins . In recent years, multi-omics analysis techniques have gained widespread application in the study of abiotic stress. Transcriptomics and metabolomics techniques provide a more comprehensive understanding by allowing detailed monitoring of metabolic regulation and molecular processes in plants exposed to biotic or abiotic stress environments . Cheng et al. employed multi-omics analysis to identify 18 significant metabolites, two key pathways, and six critical genes responding to low-temperature stress in Helicotrichon virescens . Zhao et al. conducted a comprehensive analysis integrating transcriptomics and metabolomics to explore the alterations in genes and metabolites of cold-tolerant wheat under low-temperature stress. Their findings revealed the pivotal roles of key pathways associated with ABA/JA signaling and proline biosynthesis in regulating wheat cold tolerance . Therefore, multi-omics methodologies offer a comprehensive approach to elucidate cellular life processes from diverse dimensions, enhancing our understanding of the potential mechanisms underlying plant stress resistance. Phaseolus vulgaris L. (common bean), native to Central and Southern America, is the most extensively planted, cultivated, and consumed legume globally. Renowned for its high protein content and abundant nutrients, it is a vital source of plant-based protein for human consumption . Common beans thrive in warm conditions and are susceptible to frost. Low temperatures below 10 ℃ can significantly impede their growth and development, rendering them a cold-sensitive vegetable crop . To delineate the variances between cold-resistant and cold-sensitive common bean materials, we conducted a comprehensive analysis integrating phenotype physiological assessments with multi-omics analysis to unveil significant changes in physiological parameters, genes, and metabolites. Through comparative analysis of physiological indicators, transcriptome, and metabolomic profiles of common beans subjected to low-temperature stress, we identified flavonoid metabolism and plant hormone signal transduction as key components of common bean response to low temperature. These findings offer valuable insights into the mechanisms underlying common bean cold tolerance and contribute to the optimal utilization of cold-tolerant resources in common bean breeding programs. 

Plant material and cold treatments The seeds of the two common bean materials were provided by the Vegetable Research Laboratory, College of Horticulture, Sichuan Agricultural University. The seeds of the cold-sensitive ‘Bai Bu Lao’ (BBL) and cold-tolerant ‘Wei Yuan’ (WY) common bean varieties were soaked at room temperature for 4 h and then germinated at 25 ℃ in the absence of light. After germination, the seeds were exposed to white light and subsequently sown in nutrient bowls. The substrate composition consisted of peat, vermiculite, and perlite in a 3:1:1 ratio. The cultivation conditions included a 12 h light cycle, a day/night temperature of 25 ± 2 ℃/18 ± 2 ℃, a light intensity of 300 µ mol/m 2 ·s, and a relative humidity of 60 ± 5%. Experimental treatments were conducted when the seedlings reached the two-leaf and one-heart stage. For the cold stress treatment, the temperature in the artificial climate box was set to 5 ℃ for both day and night). The light cycle was 12 h of light and 12 h of dark, with a light intensity of 100 µ mol/m 2 ·s and relative humidity of 60 ± 5%. Common bean leaves were collected after 0, 6, and 24 h of low-temperature treatment, with three biological replicates for each time point. The samples were frozen in liquid nitrogen and stored at -80 ℃ for future use. The entire experiment consists of six treatments, namely BBL0, BBL6, BBL24, WY0, WY6, and WY24. Physiological measurements Relative conductivity was measured using a conductivity meter (Shanghai INESA Scientific Instrument Co., Ltd, DDS-307). Malondialdehyde (MDA) content was determined using the thiobarbituric acid method . Leaf samples were crushed in a mortar after adding 5 ml of precooled phosphate buffer, the homogenate was centrifuged for 20 min after low temperature grinding. The supernatant as the enzyme crude extract were stored at 4 °C for superoxide dismutase (SOD) and peroxidase (POD) activities assay. SOD activity was measured using the nitroblue tetrazolium method, POD activity using the guaiacol method . cDNA library construction, sequencing, and data analysis Transcriptome sequencing was performed by Biomarker Technologies Co., Ltd. (Beijing, China). Each treatment included three biological replicates, resulting in 18 samples. First, RNA was extracted, followed by mRNA purification and fragmentation. cDNA was then synthesized, ligated, and the ligation products were purified. After fragment selection, the library was constructed. Following a quality inspection of the libraries, sequencing was performed in the PE150 mode using a high-throughput sequencing platform (Illumina NovaSeq 6000 platform, San Diego, USA). After sequencing, data analysis was conducted using the bioinformatics pipeline provided by BMKCloud ( www.biocloud.net ). The raw data was filtered to obtain clean data, which were then aligned to the common bean reference genome ( https://phytozome-next.jgi.doe.gov/info/Pvulgaris_v2_1 ) using HISAT2 software. Subsequently, the reads from the comparison pairs were assembled, and the transcriptome was reconstructed using String Tie for subsequent analysis. To detect differentially expressed genes, a | Fold Change ≥ 2 | and false discovery rate (FDR) < 0.01 were employed as screening criteria. The differentially expressed genes (DEGs) were compared against the Gene Ontology (GO) database ( http://www.geneontology.org/ ) to annotate their functional terms, and the number of genes associated with each GO term was tabulated . Functional categories enriched with significant DEGs compared to the genomic background were identified using FDR ≤ 0.05 as the threshold for significant enrichment. Next, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database ( http://www.genome.jp/tools/kaas/ ) was utilized to annotate and classify the DEGs for pathway function . The enrichment results were analyzed using the hypergeometric test method with the ClusterProfiler package, and visualized using bubble plots and bar charts. Significant enrichment of DEGs in pathways was determined using FDR ≤ 0.05 as the threshold. Real-Time quantitative polymerase chain reaction (RT-qPCR) validation Briefly, 12 genes were randomly selected from the differentially expressed genes (DEGs) to validate the sequencing results using RT-qPCR. Total RNA extraction was carried out using the RNA preparation kit from Tiangen Biochemical Technology Co., Ltd. (Beijing, China), followed by cDNA synthesis using the PrimeScript™ FAST RT reagent kit with gDNA Eraser kit (TaKaRa, Japan). Gene primers were designed using SnapGene software, and their specificity was confirmed using the NCBI database. The primers were synthesized by Shenggong Biotechnology Co., Ltd. (Shanghai, China), and the specific primer sequences are provided in Table . RT-qPCR analysis was performed using the 2× SYBR qPCR Mix from Jiangsu Baishimei Biotechnology Co., Ltd (Lianyungang, China) on a Bio-Rad CFX96 PCR (Bio-Rad, USA) instrument. The reaction volume was 20 µl, consisting of 10 µl 2× SYBR qPCR Mix, 0.5 µl each of forward and reverse primers (each 10 µ M), 1 µl cDNA, and 8 µl ddH 2 O. The RT-qPCR reaction followed a two-step amplification method. The reaction procedure included pre-denaturation at 95 ℃ for 30s, followed by denaturation at 95 ℃ for 10s, and annealing at 60 ℃ for 30s, for a total of 39 cycles. In addition, a dissolution curve analysis was performed. Notably, the expression data were normalized using Actin-11 as the reference gene , and the relative expression data were calculated using the 2 −∆∆Ct method . Metabolomic profiling The qualitative and quantitative analysis of metabolites was performed based on the self-built database GB-PLANT of Beijing Biomarker Technologies Co., LTD. The metabolites in the samples were analyzed using mass spectrometry for both qualitative and quantitative determination. Characteristic ions of each substance were selected through triple quadrupole screening, and the signal intensity of these characteristic ions was recorded in the detector. After obtaining the mass spectrometry data of metabolites from different samples, the peak areas of all substance mass spectrometry peaks were integrated, and the peak integration was corrected for the same metabolite across different samples. Metabolites Extraction: The sample extracts were analyzed using an UPLC-ESI-MS/MS system (UPLC, Waters Acquity I-Class PLUS; MS, Applied Biosystems QTRAP 6500+). The analytical conditions were as follows, UPLC: column, Waters HSS-T3 (1.8 μm, 2.1 mm * 100 mm); The mobile phase was consisted of solvent A, pure water with 0.1% formic acid and 5mM Ammonium acetate, and solvent B, acetonitrile with 0.1% formic acid. Sample measurements were performed with a gradient program that employed the starting conditions of 98% A, 2% B and kept for 1.5 min. Within 5.0 min, a linear gradient to 50% A, 50% B was programmed, Within 9.0 min, a linear gradient to 2% A, 98% B was programmed, and a composition of 2% A, 98% B was kept for 1 min. Subsequently, a composition of 98% A, 2% B was adjusted within 1 min and kept for 3 min. The flow velocity was set as 0.35 mL per minute; The column oven was set to 50 °C; The injection volume was 2 uL. The effluent was alternatively connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS. The sample extracts were analyzed using a UPLC-ESI-MS/MS system (UPLC, Waters Acquity I-Class PLUS; MS, Applied Biosystems QTRAP 6500+). Data analysis: After normalization of the original peak area information to the total peak area, subsequent analyses were performed. Principal component analysis (PCA) and Spearman correlation analysis were employed to assess the repeatability of the samples within groups and the quality control samples. The identified compounds were classified, and their pathway information was determined using KEGG, HMDB, and LIPID MAPS databases. Based on the grouping information, fold changes were calculated and compared. The significant differences for each compound were assessed using a t-test to calculate the p-value. OPLS-DA modeling was performed using the R language package “ropls”, and the reliability of the model was verified through 200 times permutation tests. The VIP value of the model was calculated using multiple cross-validation. Differential metabolites (DEMs) were screened by combining the fold changes, p-values, and VIP values of the OPLS-DA model. The screening criteria were FC > 1, p-value < 0.05, and VIP > 1. The significance of KEGG pathway enrichment for the differential metabolites was calculated using the hypergeometric distribution test. Integrative analysis of transcriptome and metabolome The DEGs and DEMs with a Pearson correlation coefficient (r) > |0.8| was selected to establish a correlation network, which was visually analyzed using Cytoscape (v.3.10.2) software . Statistical analyses Physiological data were organized using Excel 2019 software. Single factor analysis (ANOVA) was performed using SPSS 27.0 software, with multiple comparisons conducted using the Duncan method. The significance threshold was set at P < 0.05. Figures were generated using Origin 2021 software. Data are represented as the mean ± standard error (SE) of three biological replicates. 

The seeds of the two common bean materials were provided by the Vegetable Research Laboratory, College of Horticulture, Sichuan Agricultural University. The seeds of the cold-sensitive ‘Bai Bu Lao’ (BBL) and cold-tolerant ‘Wei Yuan’ (WY) common bean varieties were soaked at room temperature for 4 h and then germinated at 25 ℃ in the absence of light. After germination, the seeds were exposed to white light and subsequently sown in nutrient bowls. The substrate composition consisted of peat, vermiculite, and perlite in a 3:1:1 ratio. The cultivation conditions included a 12 h light cycle, a day/night temperature of 25 ± 2 ℃/18 ± 2 ℃, a light intensity of 300 µ mol/m 2 ·s, and a relative humidity of 60 ± 5%. Experimental treatments were conducted when the seedlings reached the two-leaf and one-heart stage. For the cold stress treatment, the temperature in the artificial climate box was set to 5 ℃ for both day and night). The light cycle was 12 h of light and 12 h of dark, with a light intensity of 100 µ mol/m 2 ·s and relative humidity of 60 ± 5%. Common bean leaves were collected after 0, 6, and 24 h of low-temperature treatment, with three biological replicates for each time point. The samples were frozen in liquid nitrogen and stored at -80 ℃ for future use. The entire experiment consists of six treatments, namely BBL0, BBL6, BBL24, WY0, WY6, and WY24. 

Relative conductivity was measured using a conductivity meter (Shanghai INESA Scientific Instrument Co., Ltd, DDS-307). Malondialdehyde (MDA) content was determined using the thiobarbituric acid method . Leaf samples were crushed in a mortar after adding 5 ml of precooled phosphate buffer, the homogenate was centrifuged for 20 min after low temperature grinding. The supernatant as the enzyme crude extract were stored at 4 °C for superoxide dismutase (SOD) and peroxidase (POD) activities assay. SOD activity was measured using the nitroblue tetrazolium method, POD activity using the guaiacol method . 

Transcriptome sequencing was performed by Biomarker Technologies Co., Ltd. (Beijing, China). Each treatment included three biological replicates, resulting in 18 samples. First, RNA was extracted, followed by mRNA purification and fragmentation. cDNA was then synthesized, ligated, and the ligation products were purified. After fragment selection, the library was constructed. Following a quality inspection of the libraries, sequencing was performed in the PE150 mode using a high-throughput sequencing platform (Illumina NovaSeq 6000 platform, San Diego, USA). After sequencing, data analysis was conducted using the bioinformatics pipeline provided by BMKCloud ( www.biocloud.net ). The raw data was filtered to obtain clean data, which were then aligned to the common bean reference genome ( https://phytozome-next.jgi.doe.gov/info/Pvulgaris_v2_1 ) using HISAT2 software. Subsequently, the reads from the comparison pairs were assembled, and the transcriptome was reconstructed using String Tie for subsequent analysis. To detect differentially expressed genes, a | Fold Change ≥ 2 | and false discovery rate (FDR) < 0.01 were employed as screening criteria. The differentially expressed genes (DEGs) were compared against the Gene Ontology (GO) database ( http://www.geneontology.org/ ) to annotate their functional terms, and the number of genes associated with each GO term was tabulated . Functional categories enriched with significant DEGs compared to the genomic background were identified using FDR ≤ 0.05 as the threshold for significant enrichment. Next, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database ( http://www.genome.jp/tools/kaas/ ) was utilized to annotate and classify the DEGs for pathway function . The enrichment results were analyzed using the hypergeometric test method with the ClusterProfiler package, and visualized using bubble plots and bar charts. Significant enrichment of DEGs in pathways was determined using FDR ≤ 0.05 as the threshold. 

Briefly, 12 genes were randomly selected from the differentially expressed genes (DEGs) to validate the sequencing results using RT-qPCR. Total RNA extraction was carried out using the RNA preparation kit from Tiangen Biochemical Technology Co., Ltd. (Beijing, China), followed by cDNA synthesis using the PrimeScript™ FAST RT reagent kit with gDNA Eraser kit (TaKaRa, Japan). Gene primers were designed using SnapGene software, and their specificity was confirmed using the NCBI database. The primers were synthesized by Shenggong Biotechnology Co., Ltd. (Shanghai, China), and the specific primer sequences are provided in Table . RT-qPCR analysis was performed using the 2× SYBR qPCR Mix from Jiangsu Baishimei Biotechnology Co., Ltd (Lianyungang, China) on a Bio-Rad CFX96 PCR (Bio-Rad, USA) instrument. The reaction volume was 20 µl, consisting of 10 µl 2× SYBR qPCR Mix, 0.5 µl each of forward and reverse primers (each 10 µ M), 1 µl cDNA, and 8 µl ddH 2 O. The RT-qPCR reaction followed a two-step amplification method. The reaction procedure included pre-denaturation at 95 ℃ for 30s, followed by denaturation at 95 ℃ for 10s, and annealing at 60 ℃ for 30s, for a total of 39 cycles. In addition, a dissolution curve analysis was performed. Notably, the expression data were normalized using Actin-11 as the reference gene , and the relative expression data were calculated using the 2 −∆∆Ct method . 

The qualitative and quantitative analysis of metabolites was performed based on the self-built database GB-PLANT of Beijing Biomarker Technologies Co., LTD. The metabolites in the samples were analyzed using mass spectrometry for both qualitative and quantitative determination. Characteristic ions of each substance were selected through triple quadrupole screening, and the signal intensity of these characteristic ions was recorded in the detector. After obtaining the mass spectrometry data of metabolites from different samples, the peak areas of all substance mass spectrometry peaks were integrated, and the peak integration was corrected for the same metabolite across different samples. Metabolites Extraction: The sample extracts were analyzed using an UPLC-ESI-MS/MS system (UPLC, Waters Acquity I-Class PLUS; MS, Applied Biosystems QTRAP 6500+). The analytical conditions were as follows, UPLC: column, Waters HSS-T3 (1.8 μm, 2.1 mm * 100 mm); The mobile phase was consisted of solvent A, pure water with 0.1% formic acid and 5mM Ammonium acetate, and solvent B, acetonitrile with 0.1% formic acid. Sample measurements were performed with a gradient program that employed the starting conditions of 98% A, 2% B and kept for 1.5 min. Within 5.0 min, a linear gradient to 50% A, 50% B was programmed, Within 9.0 min, a linear gradient to 2% A, 98% B was programmed, and a composition of 2% A, 98% B was kept for 1 min. Subsequently, a composition of 98% A, 2% B was adjusted within 1 min and kept for 3 min. The flow velocity was set as 0.35 mL per minute; The column oven was set to 50 °C; The injection volume was 2 uL. The effluent was alternatively connected to an ESI-triple quadrupole-linear ion trap (QTRAP)-MS. The sample extracts were analyzed using a UPLC-ESI-MS/MS system (UPLC, Waters Acquity I-Class PLUS; MS, Applied Biosystems QTRAP 6500+). Data analysis: After normalization of the original peak area information to the total peak area, subsequent analyses were performed. Principal component analysis (PCA) and Spearman correlation analysis were employed to assess the repeatability of the samples within groups and the quality control samples. The identified compounds were classified, and their pathway information was determined using KEGG, HMDB, and LIPID MAPS databases. Based on the grouping information, fold changes were calculated and compared. The significant differences for each compound were assessed using a t-test to calculate the p-value. OPLS-DA modeling was performed using the R language package “ropls”, and the reliability of the model was verified through 200 times permutation tests. The VIP value of the model was calculated using multiple cross-validation. Differential metabolites (DEMs) were screened by combining the fold changes, p-values, and VIP values of the OPLS-DA model. The screening criteria were FC > 1, p-value < 0.05, and VIP > 1. The significance of KEGG pathway enrichment for the differential metabolites was calculated using the hypergeometric distribution test. 

The DEGs and DEMs with a Pearson correlation coefficient (r) > |0.8| was selected to establish a correlation network, which was visually analyzed using Cytoscape (v.3.10.2) software . 

Physiological data were organized using Excel 2019 software. Single factor analysis (ANOVA) was performed using SPSS 27.0 software, with multiple comparisons conducted using the Duncan method. The significance threshold was set at P < 0.05. Figures were generated using Origin 2021 software. Data are represented as the mean ± standard error (SE) of three biological replicates. 

Cold stress significantly affected the phenotype and physiological indicators of two common bean varieties The phenotypes and associated physiological indices of the cold-sensitive variety ‘Bai Bu Lao’ (BBL) and the cold-tolerant variety ‘Wei Yuan’ (WY) were analyzed following exposure to 5 °C for 0, 6, and 24 h. As shown in Fig. A, under normal temperature conditions (0 h), the two bean varieties had no significant phenotypic differences. After 6 h of exposure to 5 °C, all leaves of the BBL beans wilted and drooped, whereas only some leaves of the WY beans exhibited wilting and drooping. The extent of cold damage was significantly greater in BBL beans compared to WY beans. After 24 h of exposure to 5 °C, the leaves of the BBL beans continued to droop and gradually became wrinkled, showing clear signs of dehydration. In contrast, the leaves of the WY beans ceased drooping and began to flatten, gradually resuming growth. Furthermore, the physiological indices of BBL and WY beans under low-temperature stress were measured. After 24 h of exposure to 5 °C, compared to the 0 h time point, the relative conductivity of BBL and WY beans increased by 79.08% and 108.34% ( P < 0.05) (Fig. B), the malondialdehyde content increased by 51.54% and 85.76% (Fig. C), and the SOD activity increased by 16.72% and 33.86% ( P < 0.05), respectively (Fig. D). Compared with the 0 h time point, exposure to 5 °C for 6 h and 24 h significantly increased the relative conductivity, malondialdehyde content, SOD activity, and POD activity in both BBL and WY beans (Fig. C, D, E). Notably, the increases in these physiological indices were more significant in WY beans than in BBL beans. Collectively, these results indicate a difference in cold tolerance between BBL and WY beans, with WY beans demonstrating greater cold-tolerance than BBL beans. RNA sequencing, assembly, and real-time quantitative polymerase chain reaction (RT-qPCR) validation Transcription group sequencing of BBL and WY beans under low-temperature treatment for 0 h, 6 h, and 24 h was performed using a high-flow sequencing platform. Each sample was analyzed in triplicate at each time point, resulting in 18 samples. As shown in Table , sequencing quality control yielded a total of 114.91 Gb of Clean Data, with a Q30 alkali ratio of ≥ 93.40% for each sample and a ratio between 91.43% and 96.44% for reads and the reference genome. This demonstrates that the output and quality of the archived sequencing data meet the requirements for further analysis and are suitable for subsequent bioinformatics analysis. At the same time, 12 DEGs were randomly selected for RT-qPCR analysis across the six treatments (BBL0, BBL6, BBL24, WY0, WY6, and WY24) to verify the reliability of the transcriptome data. Although some differences in gene expression were observed, as shown in Figure , the RT-qPCR results for 11 DEGs were consistent with the transcriptome expression trends, indicating that the transcriptome data from this study is reliable and suitable for further analysis. Overview of metabolic profiles in common beans under cold stress To investigate the effect of cold stress on the metabolic levels of different common bean varieties, we conducted qualitative and quantitative analyses of metabolites in the 18 samples using high-throughput and extensive targeted detection technology. A total of 923 metabolites were detected, including carboxylic acids and derivatives (12.2%), organoxygen compounds (8.5%), flavonoids (4.4%), prenol lipids (4.4%), fatty acyls (4.3%), benzene and substituted derivatives (2.3%), pure nucleosides (1.5%), and stenides (1.5%), rods and steroid derivatives (1.5%), phenols (1.4%), organitogen compounds (1.3%), and others (58.1%) (Fig. A). In the clustering heatmap analysis, BBL and WY beans, along with their respective biological replicates for each treatment, clustered together, suggesting high repeatability and correlation within the output data (Fig. B). The heat map results further reveal specific accumulation of certain metabolites in BBL beans, while others exhibit accumulation exclusively in WY beans. In addition, under cold stress, disparities in metabolite accumulation patterns between the two genotypes suggest that variations in cold tolerance within the common bean species may stem from the accumulation of DEMs. Moreover, PCA analysis was conducted on the dataset comprising 923 metabolites. The PCA score plot indicates that the first principal component (PC1) and the second principal component (PC2) account for 24.3% and 22.1% of the total variance, respectively (Fig. C). The cumulative variance explained by PC1 to PC5 is 80% (Fig. D). The close clustering of duplicate and mixed samples suggests the reproducibility and reliability of this experiment. Additionally, a Variable Importance in Projection (VIP) analysis based on the PLS-DA model was performed to identify the most informative metabolites that differentiate the two contrasting common bean varieties. The results indicated that the PLS-DA analysis effectively separated the metabolites of the two bean varieties and their three different treatment periods along the t1 and t2 axes (Figure A). The permutation test of the PLS-DA model also demonstrated a good independence between the training and testing datasets (Figure B). Analysis of the VIP values for all metabolites revealed that 428 metabolites had VIP values greater than 1 (Table ). The top twenty metabolites based on their VIP values are listed in the figure (Figure C). Joint analysis of DEGs and DEMs in common beans under cold stress This experiment established seven groups for comparison, including intra-group comparisons of BBL and WY and inter-group comparisons of BBL and WY. Initially, the study scrutinized the influence of cold stress on different common bean varieties. In the cold-sensitive common bean BBL, compared to the 0 h time point, 1923 upregulated and 1937 downregulated DEGs were identified after 6 h of low-temperature treatment (Fig. A). In addition, there were 177 upregulated and 192 downregulated DEMs in BBL (Fig. B). The KEGG joint enrichment analysis revealed significant enrichment of DEGs in pathways such as Isoflavonoid biosynthesis, Flavonoid biosynthesis, beta-Alanine metabolism, valine, leucine and isoleucine degradation, and phenylpropanoid biosynthesis pathways. Additionally, DEMs were significantly enriched in the Lysine degradation and Star and cross-metabolism pathways (Fig. C, Figure ). The GO enrichment analysis of DEGs identified annotations of 19 biological processes, three cellular components, and 12 molecular functions (Fig. G). After 24 h of low-temperature treatment, there are 3487 upregulated and 3489 downregulated DEGs in BBL (Fig. A). In addition, there were 252 upregulated and 231 downregulated DEMs in BBL (Fig. B). The KEGG joint enrichment analysis revealed significant enrichment of DEGs in pathways such as Plant hormone signal transduction, Flavonoid biosynthesis, Isoflavonoid biosynthesis, Galactose metabolism, Benzoxazinoid biosynthesis, Zeatin biosynthesis, and Phenolpropanoid biosynthesis. On the other hand, DEMs were significantly enriched in Cyanoamino acid metabolism and valine, Leucine, and isoleucine degradation (Fig. D, Figure ). The GO enrichment analysis for DEGs identified annotations for 20 biological processes, three cellular components, and 13 molecular functions (Figure A). In the cold-resistant common bean WY, compared to the 0 h time point, 1290 upregulated and 1359 downregulated DEGs were identified after 6 h of low-temperature treatment (Fig. A). Additionally, there were 180 upregulated and 244 downregulated DEMs in WY (Fig. B). The KEGG joint enrichment analysis revealed significant enrichment of DEGs in pathways such as Flavonoid biosynthesis, Plant hormone signal translation, Stilbenoid, diarylheptanoid and ginger biosynthesis, Flavone and flavonol biosynthesis, and Phenolpropanoid biosynthesis. DEMs were significantly enriched in Glyoxylate and dicarboxylate metabolism (Fig. E, Figure ). The GO enrichment analysis of DEGs identified annotations for 20 biological processes, three cellular components, and 12 molecular functions (Figure A). After 24 h of low-temperature treatment, there were 3762 upregulated and 3566 downregulated DEGs in WY (Fig. A). In addition, there were 182 upregulated and 256 downregulated DEMs in WY (Fig. B). The KEGG combined enrichment analysis revealed significant enrichment of DEGs in pathways such as isoflavone biosynthesis, Isoflavonoid biosynthesis, Plant hormone signal translation, Flavonoid biosynthesis, Galactose metabolism, Flavone and Flavonol biosynthesis, Carotenoid biosynthesis, valine, leucine and isoleucine degradation, and Glutathione metabolism (Fig. F, Figure ). Furthermore, the GO enrichment analysis of DEGs identified annotations for 20 biological processes, three cellular components, and 13 molecular functions (Figure B). Integrating the enriched metabolic pathways with differential genes and metabolites after low-temperature treatment in BBL and WY, most of the pathways are associated with sugar and acid metabolism, amino acid metabolism, flavonoid biosynthesis metabolism, and other metabolic pathways. Moreover, pathways such as plant hormone signal transduction, photosynthesis and starch, and sucrose metabolism were also identified. Notably, no significant difference was observed in GO enrichment between BBL and WY common beans after low-temperature treatment. Joint analysis of DEGs and DEMs in cold-sensitive BBL and cold-tolerant WY under cold stress Furthermore, we compared the effects of cold stress on two common bean varieties with differing cold tolerance. At 0 h, the cold-tolerant variety WY had 937 upregulated DEGs and 1170 down-regulated DEGs compared to the cold-sensitive bean variety BBL (Fig. A). There were also 229 upregulated DEMs and 196 down-regulated DEMs (Fig. B). KEGG pathway enrichment analysis revealed that DEGs were significantly enriched in pathways related to galactose metabolism, ascorbate and aldarate metabolism, valine, leucine, and isoleucine degradation, plant hormone signal transduction, isoquinoline alkaloid biosynthesis, and phenylalanine metabolism (Figure A, B). DEMs were significantly enriched in pathways such as fatty acid degradation, isoflavonoid biosynthesis, and cyanoamino acid metabolism (Figure A, B). After 6 h of low-temperature treatment, there were 1396 upregulated DEGs and 1384 downregulated DEGs in BBL compared to WY (Fig. A). Additionally, BBL had 249 upregulated DEMs and 240 downregulated DEMs compared to WY (Fig. B). The KEGG pathway enrichment analysis showed that DEGs were significantly enriched in pathways related to flavonoid biosynthesis, isoflavonoid biosynthesis, phenylpropanoid biosynthesis, amino sugar and nucleotide sugar metabolism, ascorbate and aldarate metabolism, galactose metabolism, and starch and sucrose metabolism. DEMs were significantly enriched in pathways such as isoflavonoid biosynthesis, isoquinoline alkaloid biosynthesis, and alpha-linolenic acid metabolism (Fig. A, Figure C). After 24 h of low-temperature treatment, there were 995 upregulated DEGs and 1163 downregulated DEGs in BBL compared to WY (Fig. A); There were also 194 upregulated DEMs and 282 downregulated DEMs in BBL compared to WY (Fig. B). The KEGG pathway enrichment analysis revealed significant enrichment of DEGs in pathways such as galactose metabolism, ascorbate and aldarate metabolism, amino sugar and nucleotide sugar metabolism, alpha-linolenic acid metabolism, isoflavonoid biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, flavonoid biosynthesis, and pentose and glucuronate interconversions (Fig. B, Figure D). DEMs were significantly enriched in pathways such as ascorbate and aldarate metabolism, isoflavonoid biosynthesis, and valine, leucine, and isoleucine degradation (Fig. B, Figure D). These findings indicate that at 6 h of low-temperature treatment, the cold-sensitive BBL and cold-tolerant WY common bean varieties exhibited the highest number of DEGs and DEMs. These were significantly enriched in pathways associated with flavonoid biosynthesis, sugar acid metabolism, amino acid metabolism, and plant hormone signal transduction. Therefore, in subsequent experiments, our focus was directed toward analyzing the response of common beans to low-temperature stress, specifically examining sugars and amino acids, flavonoid secondary metabolites, plant hormones, and associated genes. The DEGs petalogram results show that following low-temperature treatment, there was one upregulated DEG, Phvul.007G135400 (Fig. C), shared across all comparison groups of the cold-sensitive BBL cowpea and cold-tolerant WY common bean. There were also 25 downregulated DEGs (Fig. D). Moreover, the DEMs petalogram results reveal that following low-temperature treatment, three upregulated DEMs (Fig. D) were shared among all comparison groups of BBL and WY common bean, namely L-Tryptophan, Vasicinol, and Shizukaol D. Notably, there were eight down-regulated DEMs (Fig. E). Analysis of transcription factors in DEGs Transcription factors play a crucial regulatory role in the response mechanism to cold stress. To elucidate this regulatory network, transcription factor prediction was performed on all DEGs (11212) identified from transcriptome sequencing. The predicted results showed that the top six transcription factor families with the most annotations were RLK-Pelle-DLSV (93), AP2/ERF-ERF (76), bHLH (71), MYB (64), WRKY (56), and NAC (54), respectively (Fig. A). Most of the DEGs in transcription factors are implicated in the response to low-temperature stress, including bHLH and MYB families, which have been previously reported to regulate plant cold tolerance . The bHLH transcription factor plays a multifaceted role, regulating physiological and biochemical processes such as signal transduction. In addition, it is actively involved in modulating various stress responses, including those to drought and cold stress. This study identified 71 bHLH gene family transcription factors from the pool of DEGs (Fig. B). Among them, the cold stress treatment for 24 h significantly upregulated the expression of six bHLH transcription factors (Phvul.003G181900.v2.1, Phvul.002G018300.v2.1, Phvul.003G140800.v2.1, Phvul.002G283100.v2.1, Phvul.001G121200.v2.1, and Phvul.003G157100.v2.1) in cold-sensitive BBL beans. Cold stress for 6 h significantly upregulated the expression of six bHLH transcription factors (Phvul.009G137400.v2.1, Phvul.010G120000.v2.1, Phvul.002G017600.v2.1, Phvul.006G196600.v2.1, Phvul.001G085500.v2.1, and Phvul.002G088400.v2.1) in cold-tolerant WY beans. Simultaneously, low-temperature treatment for 24 h significantly downregulated the expression of five bHLH transcription factors (Phvul.003G067500.v2, Phvul.006G028500.v2.1, Phvul.006G184600.v2.1, Phvul.001G031300.v2.1, and Phvul.009G023500.v2.1) in cold-sensitive BBL and cold-tolerant WY beans. In addition, Phvul.001G023200.v2.1 and Phvul.003G231200.v2.1 exhibited contrasting trends between cold-sensitive BBL and cold-tolerant WY beans under cold stress conditions (Fig. B). The MYB transcription factor family is one of the largest in plants, and numerous research findings support its potential as a transcription factor for plant breeding and enhancement. In this study, 64 MYB gene family transcription factors were identified from the pool of DEGs (Fig. C). Among them, low-temperature treatment for 24 h significantly downregulated the expression of six MYB transcription factors (Phvul.002G170500.v2.1, Phvul.010G009900.v2.1, Phvul.003G176800.v2.1, Phvul.008G262900.v2.1, and Phvul.009G187700.v2.1) in cold-sensitive BBL and cold-tolerant WY beans, while significantly upregulating the expression of 25 MYB transcription factors in both varieties. At the same time, the changing trends of Phvul.007G108500-v2.1, Phvul.008G041500-v2.1, and Phvul.007G215800-v2.1 in BBL and WY beans after 24 h of low-temperature treatment are inconsistent. The effect of cold stress on sugar and amino acids in common beans Previous research has demonstrated that low-temperature stress significantly impacts sugar metabolism, amino acid synthesis, and metabolic pathways in plants . Therefore, we integrated DEGs and DEMs to construct pathway diagrams for sugar and amino acid metabolism. As shown in Fig. , low-temperature stress significantly increased the levels of galactinol, sucrose, melibiose, trehalose, and maltose in both cold-sensitive BBL common beans and cold-tolerant WY cowpeas, while it decreased the level of fructose-6P. Additionally, the expression patterns of raffinose, manninotriose, D-galactose, D-fructose, glucose, and glucose-6P differed inconsistently between cold-sensitive BBL common beans and cold-tolerant WY cowpeas (Fig. A). Concurrently, the expression levels of genes encoding myo-inositol galactosyltransferase (Phvul.001G223700.v2.1, Phvul.007G203400.v2.1), cotton galactose synthase (Phvul.004G007100.v2.1), hexokinase (Phvul.010G144900.v2.1), fructokinase (Phvul.011G091600.v2.1), pyruvate kinase (Phvul.004G000800.v2.1), phosphofructokinase/phosphatase (Phvul.003G150400.v2.1, Phvul.008G172400.v2.1, Phvul.009G054400.v2.1), α-amylase (Phvul.006G185000.v2.1), and β-amylase (Phvul.003G226900.v2.1) also showed a marked increase (Fig. B). Upon entry into the citric acid cycle, low-temperature stress significantly reduced the levels of 2-ketoglutaric acid in cold-sensitive BBL common beans and cold-tolerant WY common beans and influenced the aconitic acid content. The aconitic acid content decreased significantly with prolonged exposure to low temperature in BBL common beans, whereas in WY common beans, it initially decreased and then increased (Fig. A). Additionally, the expression levels of genes encoding isocitrate dehydrogenase (Phvul.002G114600.v2.1, Phvul.007G150400.v2.1) increased significantly with prolonged low-temperature treatment (Fig. B). Under low-temperature stress conditions, the levels of certain amino acids and their downstream metabolites, including proline, L-leucine, L-valine, tyrosine, D-aspartic acid, and L-aspartate, significantly increased in both cold-sensitive BBL common beans and cold-tolerant WY common beans. The levels of L-lysine, glutamine, and citrulline in both cold-sensitive BBL common beans and cold-tolerant WY common beans exhibited a pattern of initial increase followed by subsequent decrease with prolonged exposure to low temperature. Furthermore, the alterations in the levels of glutamate, GABA, citruline, arginine, L-isoleucine, tyramine, and L-asparagine exhibited divergent trends between cold-sensitive BBL common beans and cold-tolerant WY common beans (Fig. A). Meanwhile, the expression levels of genes encoding 2-Oxoglutarate dehydrogenase (Phvul.010G017700.v2.1, Phvul.009G156400.v2.1), isobutyryl-CoA dehydrogenase (Phvul.011G073500.v2.1), and 3-Methylcrotonyl-CoA carboxylase (Phvul.003G291600.v2.1) significantly decreased with prolonged low-temperature treatment in both cold-sensitive BBL common beans and cold-tolerant WY common beans. The expression levels of genes encoding tyrosine metabolism enzyme (Phvul.006G143500.v2.1) and aspartate aminotransferase (Phvul.001G250600.v2.1) increased with prolonged low-temperature treatment in both cold-sensitive BBL common beans and cold-tolerant WY common beans (Fig. B). The effect of cold stress on secondary metabolites of flavonoids in common beans Flavonoids constitute a crucial category of secondary metabolites in plants, exerting significant biological functions. Herein, we identified 19 types of DEMs and 15 types of DEGs associated with the biosynthesis and metabolism of flavonoids in common beans. As shown in Fig. , cold stress significantly augmented the levels of L-tyrosine, L-phenylalanine, naringenin chalcone, naringenin, genistein, apigenin, apigenin 7-Glucoside, and pratenstein, while significantly diminishing the content of coumestrol and genistein 7,4’-Di-O-β-D-glucopyranoside. These substances exhibited consistent trends in both BBL and WY beans. The levels of eriodictyol, quercetin, isoquercitrin, isosakuranetin, isoliquiritigenin, 7,4 ‘- Dihydroxyflavone, trans-5-O - (p-Coumaroyl) shikimate, luteolin, and chloroiol exhibited distinct trends between BBL and WY beans (Fig. A). At the same time, cold stress significantly upregulated the expression of genes encoding shikimate O-hydroxycinnamoyl transferase (the genes numbered 4 and marked in red in Fig. B), chalcone synthase (the genes numbered 5 and marked in red in Fig. B), and isoflavone 7-O-glucose-6 ‘’- O-malonyltransferase (Phvul.008G029400.v2.1, Phvul.008G032200.v2.1), and exhibited a consistent trend of change in BBL beans and WY beans. The gene encoding phenylalanine ammonia-lyase (Phvul.001G177700.v2.1, Phvul.001G177800.v2.1) displayed divergent trends of change between BBL and WY beans. Exposure to cold stress for 6 h significantly upregulated the expression of the PAL gene in cold-sensitive bean BBL. In addition, cold stress for 24 h significantly upregulated the expression of genes encoding 5-O - (4-coumaroyl) - D-quinoate 3 ‘- monooxygenase and chalcone synthase (the genes numbered 6–7 and marked in red in Fig. B), 2-hydroxyisoflavone synthase (Phvul.003G074000.v2.1, Phvul.003G051800.v2.1, Phvul.003G051801.v2.1, Phvul.003G051700.v2.1), and flavonoid 3’ – monooxygenase (Phvul.L001623.v2.1) in WY beans (Fig. B). These findings suggest that common beans may alleviate the effects of cold stress on plants by modulating the expression of flavonoids and genes encoding flavonoid synthesis. The effect of cold stress on hormone synthesis and transduction in common beans Plant hormones play an important role in mediating response to cold stress. In this study, we identified three DEMs and 18 DEGs implicated in the biosynthesis pathways of brassinolide (BR), abscisic acid (ABA), jasmonic acid (JA), and salicylic acid (SA) (Fig. ). In the BR pathway, the expression levels of Phvul.001G075500.v2.1, Phvul.002G318200.v2.1, Phvul.003G187200.v2.1, Phvul.003G247400.v2.1, Phvul.006G033300.v2.1, Phvul.002G047200.v2.1, and Phvul.003G143332.v2.1 significantly increased with the duration of cold stress. The expression levels of Phvul.008G186900.v2.1 and Phvul.010G056200.v2.1 significantly decreased with prolonged cold stress duration. In addition, the expression levels of Phvul.003G164800.v2.1, Phvul.003G247601.v2.1, Phvul.003G247651.v2.1, and Phvul.005G074000.v2.1 exhibited distinct trends in BBL and WY common beans (Fig. B). In the ABA pathway, metabolomic analysis revealed that cold stress for 24 h increased the ABA content in the leaves of both BBL and WY bean varieties (Fig. A). During the ABA signal transduction process, we identified five DEGs encoding PYP/PYL, 10 DEGs encoding PP2C, four DEGs encoding SnRK2, and five DEGs encoding ABF. Among these, the expression level of Phvul.002G141901.v2.1 decreases with prolonged cold stress duration, while the expression levels of most other DEGs increase with increasing cold stress duration. Furthermore, Phvul.001G246300.v2.1 exhibited differential changes only in cold-sensitive BBL beans, with no significant change observed in cold-tolerant WY beans (Fig. C). In the JA pathway, the biosynthesis of JA is achieved through the linolenic acid metabolic pathway, which is synthesized through the catalysis of 13-HPOT and OPDA enzymes. The metabolic results indicate that the content of α- linolenic acid only significantly decreased in cold-sensitive bean BBL, while the variation was minor in the cold-tolerant bean WY (Fig. A). During the JA signal transduction process, two DEGs encoding JAR1, one DEG encoding COI1, four DEGs encoding JAZ, and 11 DEGs encoding MYC2 were identified. Among them, the expression level of Phvul.003G231000.v2.1 decreased with prolonged low-temperature treatment duration. Moreover, the expression levels of Phvul.009G225300.v2.1 and Phvul.006G198400.v2.1 decreased with increasing low-temperature treatment time (Fig. C). In the SA pathway, the biosynthesis of SA is catalyzed by phenylalanine ammonia-lyase (PAL) through the phenylalanine metabolism pathway. The metabolomics results revealed a significant increase in the content of L-phenylalanine with prolonged cold stress duration, and this trend of change was similar in both cold-sensitive BBL beans and cold-tolerant WY beans. The content of SA in BBL beans increases with prolonged low-temperature treatment time, whereas in WY common beans, it exhibits a trend of initially increasing and then decreasing with prolonged low-temperature treatment time (Fig. A). During the SA signal transduction process, two DEGs encoding NPR1, two DEGs encoding TGA, and three DEGs encoding PR-1 were identified. Among them, the expression level of DEGs encoding NPR1 increased with prolonged low-temperature treatment duration. The expression level of DEGs encoding TGA decreases with prolonged low-temperature treatment duration. In addition, the expression levels of DEGs Phvul.003G148550.v2.1 and Phvul.006G197500.v2.1, encoding PR-1, increased with prolonged low-temperature treatment time, while the expression level of Phvul.006G196900.v2.1 decreased with prolonged low-temperature treatment time (Fig. C). 

The phenotypes and associated physiological indices of the cold-sensitive variety ‘Bai Bu Lao’ (BBL) and the cold-tolerant variety ‘Wei Yuan’ (WY) were analyzed following exposure to 5 °C for 0, 6, and 24 h. As shown in Fig. A, under normal temperature conditions (0 h), the two bean varieties had no significant phenotypic differences. After 6 h of exposure to 5 °C, all leaves of the BBL beans wilted and drooped, whereas only some leaves of the WY beans exhibited wilting and drooping. The extent of cold damage was significantly greater in BBL beans compared to WY beans. After 24 h of exposure to 5 °C, the leaves of the BBL beans continued to droop and gradually became wrinkled, showing clear signs of dehydration. In contrast, the leaves of the WY beans ceased drooping and began to flatten, gradually resuming growth. Furthermore, the physiological indices of BBL and WY beans under low-temperature stress were measured. After 24 h of exposure to 5 °C, compared to the 0 h time point, the relative conductivity of BBL and WY beans increased by 79.08% and 108.34% ( P < 0.05) (Fig. B), the malondialdehyde content increased by 51.54% and 85.76% (Fig. C), and the SOD activity increased by 16.72% and 33.86% ( P < 0.05), respectively (Fig. D). Compared with the 0 h time point, exposure to 5 °C for 6 h and 24 h significantly increased the relative conductivity, malondialdehyde content, SOD activity, and POD activity in both BBL and WY beans (Fig. C, D, E). Notably, the increases in these physiological indices were more significant in WY beans than in BBL beans. Collectively, these results indicate a difference in cold tolerance between BBL and WY beans, with WY beans demonstrating greater cold-tolerance than BBL beans. 

Transcription group sequencing of BBL and WY beans under low-temperature treatment for 0 h, 6 h, and 24 h was performed using a high-flow sequencing platform. Each sample was analyzed in triplicate at each time point, resulting in 18 samples. As shown in Table , sequencing quality control yielded a total of 114.91 Gb of Clean Data, with a Q30 alkali ratio of ≥ 93.40% for each sample and a ratio between 91.43% and 96.44% for reads and the reference genome. This demonstrates that the output and quality of the archived sequencing data meet the requirements for further analysis and are suitable for subsequent bioinformatics analysis. At the same time, 12 DEGs were randomly selected for RT-qPCR analysis across the six treatments (BBL0, BBL6, BBL24, WY0, WY6, and WY24) to verify the reliability of the transcriptome data. Although some differences in gene expression were observed, as shown in Figure , the RT-qPCR results for 11 DEGs were consistent with the transcriptome expression trends, indicating that the transcriptome data from this study is reliable and suitable for further analysis. 

To investigate the effect of cold stress on the metabolic levels of different common bean varieties, we conducted qualitative and quantitative analyses of metabolites in the 18 samples using high-throughput and extensive targeted detection technology. A total of 923 metabolites were detected, including carboxylic acids and derivatives (12.2%), organoxygen compounds (8.5%), flavonoids (4.4%), prenol lipids (4.4%), fatty acyls (4.3%), benzene and substituted derivatives (2.3%), pure nucleosides (1.5%), and stenides (1.5%), rods and steroid derivatives (1.5%), phenols (1.4%), organitogen compounds (1.3%), and others (58.1%) (Fig. A). In the clustering heatmap analysis, BBL and WY beans, along with their respective biological replicates for each treatment, clustered together, suggesting high repeatability and correlation within the output data (Fig. B). The heat map results further reveal specific accumulation of certain metabolites in BBL beans, while others exhibit accumulation exclusively in WY beans. In addition, under cold stress, disparities in metabolite accumulation patterns between the two genotypes suggest that variations in cold tolerance within the common bean species may stem from the accumulation of DEMs. Moreover, PCA analysis was conducted on the dataset comprising 923 metabolites. The PCA score plot indicates that the first principal component (PC1) and the second principal component (PC2) account for 24.3% and 22.1% of the total variance, respectively (Fig. C). The cumulative variance explained by PC1 to PC5 is 80% (Fig. D). The close clustering of duplicate and mixed samples suggests the reproducibility and reliability of this experiment. Additionally, a Variable Importance in Projection (VIP) analysis based on the PLS-DA model was performed to identify the most informative metabolites that differentiate the two contrasting common bean varieties. The results indicated that the PLS-DA analysis effectively separated the metabolites of the two bean varieties and their three different treatment periods along the t1 and t2 axes (Figure A). The permutation test of the PLS-DA model also demonstrated a good independence between the training and testing datasets (Figure B). Analysis of the VIP values for all metabolites revealed that 428 metabolites had VIP values greater than 1 (Table ). The top twenty metabolites based on their VIP values are listed in the figure (Figure C). 

This experiment established seven groups for comparison, including intra-group comparisons of BBL and WY and inter-group comparisons of BBL and WY. Initially, the study scrutinized the influence of cold stress on different common bean varieties. In the cold-sensitive common bean BBL, compared to the 0 h time point, 1923 upregulated and 1937 downregulated DEGs were identified after 6 h of low-temperature treatment (Fig. A). In addition, there were 177 upregulated and 192 downregulated DEMs in BBL (Fig. B). The KEGG joint enrichment analysis revealed significant enrichment of DEGs in pathways such as Isoflavonoid biosynthesis, Flavonoid biosynthesis, beta-Alanine metabolism, valine, leucine and isoleucine degradation, and phenylpropanoid biosynthesis pathways. Additionally, DEMs were significantly enriched in the Lysine degradation and Star and cross-metabolism pathways (Fig. C, Figure ). The GO enrichment analysis of DEGs identified annotations of 19 biological processes, three cellular components, and 12 molecular functions (Fig. G). After 24 h of low-temperature treatment, there are 3487 upregulated and 3489 downregulated DEGs in BBL (Fig. A). In addition, there were 252 upregulated and 231 downregulated DEMs in BBL (Fig. B). The KEGG joint enrichment analysis revealed significant enrichment of DEGs in pathways such as Plant hormone signal transduction, Flavonoid biosynthesis, Isoflavonoid biosynthesis, Galactose metabolism, Benzoxazinoid biosynthesis, Zeatin biosynthesis, and Phenolpropanoid biosynthesis. On the other hand, DEMs were significantly enriched in Cyanoamino acid metabolism and valine, Leucine, and isoleucine degradation (Fig. D, Figure ). The GO enrichment analysis for DEGs identified annotations for 20 biological processes, three cellular components, and 13 molecular functions (Figure A). In the cold-resistant common bean WY, compared to the 0 h time point, 1290 upregulated and 1359 downregulated DEGs were identified after 6 h of low-temperature treatment (Fig. A). Additionally, there were 180 upregulated and 244 downregulated DEMs in WY (Fig. B). The KEGG joint enrichment analysis revealed significant enrichment of DEGs in pathways such as Flavonoid biosynthesis, Plant hormone signal translation, Stilbenoid, diarylheptanoid and ginger biosynthesis, Flavone and flavonol biosynthesis, and Phenolpropanoid biosynthesis. DEMs were significantly enriched in Glyoxylate and dicarboxylate metabolism (Fig. E, Figure ). The GO enrichment analysis of DEGs identified annotations for 20 biological processes, three cellular components, and 12 molecular functions (Figure A). After 24 h of low-temperature treatment, there were 3762 upregulated and 3566 downregulated DEGs in WY (Fig. A). In addition, there were 182 upregulated and 256 downregulated DEMs in WY (Fig. B). The KEGG combined enrichment analysis revealed significant enrichment of DEGs in pathways such as isoflavone biosynthesis, Isoflavonoid biosynthesis, Plant hormone signal translation, Flavonoid biosynthesis, Galactose metabolism, Flavone and Flavonol biosynthesis, Carotenoid biosynthesis, valine, leucine and isoleucine degradation, and Glutathione metabolism (Fig. F, Figure ). Furthermore, the GO enrichment analysis of DEGs identified annotations for 20 biological processes, three cellular components, and 13 molecular functions (Figure B). Integrating the enriched metabolic pathways with differential genes and metabolites after low-temperature treatment in BBL and WY, most of the pathways are associated with sugar and acid metabolism, amino acid metabolism, flavonoid biosynthesis metabolism, and other metabolic pathways. Moreover, pathways such as plant hormone signal transduction, photosynthesis and starch, and sucrose metabolism were also identified. Notably, no significant difference was observed in GO enrichment between BBL and WY common beans after low-temperature treatment. 

Furthermore, we compared the effects of cold stress on two common bean varieties with differing cold tolerance. At 0 h, the cold-tolerant variety WY had 937 upregulated DEGs and 1170 down-regulated DEGs compared to the cold-sensitive bean variety BBL (Fig. A). There were also 229 upregulated DEMs and 196 down-regulated DEMs (Fig. B). KEGG pathway enrichment analysis revealed that DEGs were significantly enriched in pathways related to galactose metabolism, ascorbate and aldarate metabolism, valine, leucine, and isoleucine degradation, plant hormone signal transduction, isoquinoline alkaloid biosynthesis, and phenylalanine metabolism (Figure A, B). DEMs were significantly enriched in pathways such as fatty acid degradation, isoflavonoid biosynthesis, and cyanoamino acid metabolism (Figure A, B). After 6 h of low-temperature treatment, there were 1396 upregulated DEGs and 1384 downregulated DEGs in BBL compared to WY (Fig. A). Additionally, BBL had 249 upregulated DEMs and 240 downregulated DEMs compared to WY (Fig. B). The KEGG pathway enrichment analysis showed that DEGs were significantly enriched in pathways related to flavonoid biosynthesis, isoflavonoid biosynthesis, phenylpropanoid biosynthesis, amino sugar and nucleotide sugar metabolism, ascorbate and aldarate metabolism, galactose metabolism, and starch and sucrose metabolism. DEMs were significantly enriched in pathways such as isoflavonoid biosynthesis, isoquinoline alkaloid biosynthesis, and alpha-linolenic acid metabolism (Fig. A, Figure C). After 24 h of low-temperature treatment, there were 995 upregulated DEGs and 1163 downregulated DEGs in BBL compared to WY (Fig. A); There were also 194 upregulated DEMs and 282 downregulated DEMs in BBL compared to WY (Fig. B). The KEGG pathway enrichment analysis revealed significant enrichment of DEGs in pathways such as galactose metabolism, ascorbate and aldarate metabolism, amino sugar and nucleotide sugar metabolism, alpha-linolenic acid metabolism, isoflavonoid biosynthesis, plant hormone signal transduction, phenylpropanoid biosynthesis, flavonoid biosynthesis, and pentose and glucuronate interconversions (Fig. B, Figure D). DEMs were significantly enriched in pathways such as ascorbate and aldarate metabolism, isoflavonoid biosynthesis, and valine, leucine, and isoleucine degradation (Fig. B, Figure D). These findings indicate that at 6 h of low-temperature treatment, the cold-sensitive BBL and cold-tolerant WY common bean varieties exhibited the highest number of DEGs and DEMs. These were significantly enriched in pathways associated with flavonoid biosynthesis, sugar acid metabolism, amino acid metabolism, and plant hormone signal transduction. Therefore, in subsequent experiments, our focus was directed toward analyzing the response of common beans to low-temperature stress, specifically examining sugars and amino acids, flavonoid secondary metabolites, plant hormones, and associated genes. The DEGs petalogram results show that following low-temperature treatment, there was one upregulated DEG, Phvul.007G135400 (Fig. C), shared across all comparison groups of the cold-sensitive BBL cowpea and cold-tolerant WY common bean. There were also 25 downregulated DEGs (Fig. D). Moreover, the DEMs petalogram results reveal that following low-temperature treatment, three upregulated DEMs (Fig. D) were shared among all comparison groups of BBL and WY common bean, namely L-Tryptophan, Vasicinol, and Shizukaol D. Notably, there were eight down-regulated DEMs (Fig. E). 

Transcription factors play a crucial regulatory role in the response mechanism to cold stress. To elucidate this regulatory network, transcription factor prediction was performed on all DEGs (11212) identified from transcriptome sequencing. The predicted results showed that the top six transcription factor families with the most annotations were RLK-Pelle-DLSV (93), AP2/ERF-ERF (76), bHLH (71), MYB (64), WRKY (56), and NAC (54), respectively (Fig. A). Most of the DEGs in transcription factors are implicated in the response to low-temperature stress, including bHLH and MYB families, which have been previously reported to regulate plant cold tolerance . The bHLH transcription factor plays a multifaceted role, regulating physiological and biochemical processes such as signal transduction. In addition, it is actively involved in modulating various stress responses, including those to drought and cold stress. This study identified 71 bHLH gene family transcription factors from the pool of DEGs (Fig. B). Among them, the cold stress treatment for 24 h significantly upregulated the expression of six bHLH transcription factors (Phvul.003G181900.v2.1, Phvul.002G018300.v2.1, Phvul.003G140800.v2.1, Phvul.002G283100.v2.1, Phvul.001G121200.v2.1, and Phvul.003G157100.v2.1) in cold-sensitive BBL beans. Cold stress for 6 h significantly upregulated the expression of six bHLH transcription factors (Phvul.009G137400.v2.1, Phvul.010G120000.v2.1, Phvul.002G017600.v2.1, Phvul.006G196600.v2.1, Phvul.001G085500.v2.1, and Phvul.002G088400.v2.1) in cold-tolerant WY beans. Simultaneously, low-temperature treatment for 24 h significantly downregulated the expression of five bHLH transcription factors (Phvul.003G067500.v2, Phvul.006G028500.v2.1, Phvul.006G184600.v2.1, Phvul.001G031300.v2.1, and Phvul.009G023500.v2.1) in cold-sensitive BBL and cold-tolerant WY beans. In addition, Phvul.001G023200.v2.1 and Phvul.003G231200.v2.1 exhibited contrasting trends between cold-sensitive BBL and cold-tolerant WY beans under cold stress conditions (Fig. B). The MYB transcription factor family is one of the largest in plants, and numerous research findings support its potential as a transcription factor for plant breeding and enhancement. In this study, 64 MYB gene family transcription factors were identified from the pool of DEGs (Fig. C). Among them, low-temperature treatment for 24 h significantly downregulated the expression of six MYB transcription factors (Phvul.002G170500.v2.1, Phvul.010G009900.v2.1, Phvul.003G176800.v2.1, Phvul.008G262900.v2.1, and Phvul.009G187700.v2.1) in cold-sensitive BBL and cold-tolerant WY beans, while significantly upregulating the expression of 25 MYB transcription factors in both varieties. At the same time, the changing trends of Phvul.007G108500-v2.1, Phvul.008G041500-v2.1, and Phvul.007G215800-v2.1 in BBL and WY beans after 24 h of low-temperature treatment are inconsistent. 

Previous research has demonstrated that low-temperature stress significantly impacts sugar metabolism, amino acid synthesis, and metabolic pathways in plants . Therefore, we integrated DEGs and DEMs to construct pathway diagrams for sugar and amino acid metabolism. As shown in Fig. , low-temperature stress significantly increased the levels of galactinol, sucrose, melibiose, trehalose, and maltose in both cold-sensitive BBL common beans and cold-tolerant WY cowpeas, while it decreased the level of fructose-6P. Additionally, the expression patterns of raffinose, manninotriose, D-galactose, D-fructose, glucose, and glucose-6P differed inconsistently between cold-sensitive BBL common beans and cold-tolerant WY cowpeas (Fig. A). Concurrently, the expression levels of genes encoding myo-inositol galactosyltransferase (Phvul.001G223700.v2.1, Phvul.007G203400.v2.1), cotton galactose synthase (Phvul.004G007100.v2.1), hexokinase (Phvul.010G144900.v2.1), fructokinase (Phvul.011G091600.v2.1), pyruvate kinase (Phvul.004G000800.v2.1), phosphofructokinase/phosphatase (Phvul.003G150400.v2.1, Phvul.008G172400.v2.1, Phvul.009G054400.v2.1), α-amylase (Phvul.006G185000.v2.1), and β-amylase (Phvul.003G226900.v2.1) also showed a marked increase (Fig. B). Upon entry into the citric acid cycle, low-temperature stress significantly reduced the levels of 2-ketoglutaric acid in cold-sensitive BBL common beans and cold-tolerant WY common beans and influenced the aconitic acid content. The aconitic acid content decreased significantly with prolonged exposure to low temperature in BBL common beans, whereas in WY common beans, it initially decreased and then increased (Fig. A). Additionally, the expression levels of genes encoding isocitrate dehydrogenase (Phvul.002G114600.v2.1, Phvul.007G150400.v2.1) increased significantly with prolonged low-temperature treatment (Fig. B). Under low-temperature stress conditions, the levels of certain amino acids and their downstream metabolites, including proline, L-leucine, L-valine, tyrosine, D-aspartic acid, and L-aspartate, significantly increased in both cold-sensitive BBL common beans and cold-tolerant WY common beans. The levels of L-lysine, glutamine, and citrulline in both cold-sensitive BBL common beans and cold-tolerant WY common beans exhibited a pattern of initial increase followed by subsequent decrease with prolonged exposure to low temperature. Furthermore, the alterations in the levels of glutamate, GABA, citruline, arginine, L-isoleucine, tyramine, and L-asparagine exhibited divergent trends between cold-sensitive BBL common beans and cold-tolerant WY common beans (Fig. A). Meanwhile, the expression levels of genes encoding 2-Oxoglutarate dehydrogenase (Phvul.010G017700.v2.1, Phvul.009G156400.v2.1), isobutyryl-CoA dehydrogenase (Phvul.011G073500.v2.1), and 3-Methylcrotonyl-CoA carboxylase (Phvul.003G291600.v2.1) significantly decreased with prolonged low-temperature treatment in both cold-sensitive BBL common beans and cold-tolerant WY common beans. The expression levels of genes encoding tyrosine metabolism enzyme (Phvul.006G143500.v2.1) and aspartate aminotransferase (Phvul.001G250600.v2.1) increased with prolonged low-temperature treatment in both cold-sensitive BBL common beans and cold-tolerant WY common beans (Fig. B). 

Flavonoids constitute a crucial category of secondary metabolites in plants, exerting significant biological functions. Herein, we identified 19 types of DEMs and 15 types of DEGs associated with the biosynthesis and metabolism of flavonoids in common beans. As shown in Fig. , cold stress significantly augmented the levels of L-tyrosine, L-phenylalanine, naringenin chalcone, naringenin, genistein, apigenin, apigenin 7-Glucoside, and pratenstein, while significantly diminishing the content of coumestrol and genistein 7,4’-Di-O-β-D-glucopyranoside. These substances exhibited consistent trends in both BBL and WY beans. The levels of eriodictyol, quercetin, isoquercitrin, isosakuranetin, isoliquiritigenin, 7,4 ‘- Dihydroxyflavone, trans-5-O - (p-Coumaroyl) shikimate, luteolin, and chloroiol exhibited distinct trends between BBL and WY beans (Fig. A). At the same time, cold stress significantly upregulated the expression of genes encoding shikimate O-hydroxycinnamoyl transferase (the genes numbered 4 and marked in red in Fig. B), chalcone synthase (the genes numbered 5 and marked in red in Fig. B), and isoflavone 7-O-glucose-6 ‘’- O-malonyltransferase (Phvul.008G029400.v2.1, Phvul.008G032200.v2.1), and exhibited a consistent trend of change in BBL beans and WY beans. The gene encoding phenylalanine ammonia-lyase (Phvul.001G177700.v2.1, Phvul.001G177800.v2.1) displayed divergent trends of change between BBL and WY beans. Exposure to cold stress for 6 h significantly upregulated the expression of the PAL gene in cold-sensitive bean BBL. In addition, cold stress for 24 h significantly upregulated the expression of genes encoding 5-O - (4-coumaroyl) - D-quinoate 3 ‘- monooxygenase and chalcone synthase (the genes numbered 6–7 and marked in red in Fig. B), 2-hydroxyisoflavone synthase (Phvul.003G074000.v2.1, Phvul.003G051800.v2.1, Phvul.003G051801.v2.1, Phvul.003G051700.v2.1), and flavonoid 3’ – monooxygenase (Phvul.L001623.v2.1) in WY beans (Fig. B). These findings suggest that common beans may alleviate the effects of cold stress on plants by modulating the expression of flavonoids and genes encoding flavonoid synthesis. 

Plant hormones play an important role in mediating response to cold stress. In this study, we identified three DEMs and 18 DEGs implicated in the biosynthesis pathways of brassinolide (BR), abscisic acid (ABA), jasmonic acid (JA), and salicylic acid (SA) (Fig. ). In the BR pathway, the expression levels of Phvul.001G075500.v2.1, Phvul.002G318200.v2.1, Phvul.003G187200.v2.1, Phvul.003G247400.v2.1, Phvul.006G033300.v2.1, Phvul.002G047200.v2.1, and Phvul.003G143332.v2.1 significantly increased with the duration of cold stress. The expression levels of Phvul.008G186900.v2.1 and Phvul.010G056200.v2.1 significantly decreased with prolonged cold stress duration. In addition, the expression levels of Phvul.003G164800.v2.1, Phvul.003G247601.v2.1, Phvul.003G247651.v2.1, and Phvul.005G074000.v2.1 exhibited distinct trends in BBL and WY common beans (Fig. B). In the ABA pathway, metabolomic analysis revealed that cold stress for 24 h increased the ABA content in the leaves of both BBL and WY bean varieties (Fig. A). During the ABA signal transduction process, we identified five DEGs encoding PYP/PYL, 10 DEGs encoding PP2C, four DEGs encoding SnRK2, and five DEGs encoding ABF. Among these, the expression level of Phvul.002G141901.v2.1 decreases with prolonged cold stress duration, while the expression levels of most other DEGs increase with increasing cold stress duration. Furthermore, Phvul.001G246300.v2.1 exhibited differential changes only in cold-sensitive BBL beans, with no significant change observed in cold-tolerant WY beans (Fig. C). In the JA pathway, the biosynthesis of JA is achieved through the linolenic acid metabolic pathway, which is synthesized through the catalysis of 13-HPOT and OPDA enzymes. The metabolic results indicate that the content of α- linolenic acid only significantly decreased in cold-sensitive bean BBL, while the variation was minor in the cold-tolerant bean WY (Fig. A). During the JA signal transduction process, two DEGs encoding JAR1, one DEG encoding COI1, four DEGs encoding JAZ, and 11 DEGs encoding MYC2 were identified. Among them, the expression level of Phvul.003G231000.v2.1 decreased with prolonged low-temperature treatment duration. Moreover, the expression levels of Phvul.009G225300.v2.1 and Phvul.006G198400.v2.1 decreased with increasing low-temperature treatment time (Fig. C). In the SA pathway, the biosynthesis of SA is catalyzed by phenylalanine ammonia-lyase (PAL) through the phenylalanine metabolism pathway. The metabolomics results revealed a significant increase in the content of L-phenylalanine with prolonged cold stress duration, and this trend of change was similar in both cold-sensitive BBL beans and cold-tolerant WY beans. The content of SA in BBL beans increases with prolonged low-temperature treatment time, whereas in WY common beans, it exhibits a trend of initially increasing and then decreasing with prolonged low-temperature treatment time (Fig. A). During the SA signal transduction process, two DEGs encoding NPR1, two DEGs encoding TGA, and three DEGs encoding PR-1 were identified. Among them, the expression level of DEGs encoding NPR1 increased with prolonged low-temperature treatment duration. The expression level of DEGs encoding TGA decreases with prolonged low-temperature treatment duration. In addition, the expression levels of DEGs Phvul.003G148550.v2.1 and Phvul.006G197500.v2.1, encoding PR-1, increased with prolonged low-temperature treatment time, while the expression level of Phvul.006G196900.v2.1 decreased with prolonged low-temperature treatment time (Fig. C). 

Cold stress represents a significant abiotic stress leading to reduced yield and quality in leguminous plants. Studies have demonstrated that integrating transcriptomics and metabolomics analyses can offer comprehensive insights into the metabolic regulation and molecular mechanisms underlying cold stress . Herein, we identified 11,837 DEGs and 923 DEMs through comparative analysis of transcriptomics and widely targeted metabolomics. Our analysis unveiled their extensive involvement in the synthesis and transduction of plant metabolites and hormones. The effects of cold stress vary in terms of phenotype and physiological indicators between two different varieties of common beans Low-temperature stress has notable impacts on the osmotic regulation system and antioxidant system of plants. In Brassica napus L. (rapeseed), low-temperature stress significantly elevated the levels and relative conductivity of malondialdehyde, alongside enhancing the activity of peroxidase and superoxide dismutase . Our research findings demonstrate that low-temperature stress significantly elevated the malondialdehyde content, relative conductivity, and antioxidant enzyme activity in two distinct varieties of common beans, and the observed increments varied between the two types of common beans. These findings are consistent with those of Cai et al. in two varieties of Solanum melongena exhibiting differential cold tolerance and Wang et al. in Brassica campestris L. with varying cold tolerance levels . Together, they underscore that distinct varieties within the same species can manifest divergent responses to low-temperature stress. The role of transcription factors in the cold response of common beans Currently, a plethora of plant transcription factors have been shown to respond to low-temperature stress in plants, primarily modulating plant cold resistance through involvement in processes such as plant cell membrane fluidity, cold signal transduction, MAPK cascade, and regulation of CBF pathways . Studies have revealed that transgenic apple callus and Arabidopsis plants overexpressing MdMYB23 exhibit heightened cold tolerance . MdMYB23 has been shown to interact with the key regulator of anthocyanin biosynthesis, the MdANR promoter, thereby activating its expression, promoting anthocyanin accumulation, and clearance of reactive oxygen species (ROS) . In addition, plant alkaline helix loop helix (bHLH) transcription factors play pivotal roles in plant growth and development, secondary metabolism, and responses to abiotic stress . In rice, the upregulation of OsbHLH1 is specifically triggered by cold stress, suggesting its involvement in the cold signaling pathway of rice . In this study, a comprehensive analysis identified a total of 71 bHLH and 64 MYB gene family transcription factors (Fig. ). Interestingly, we observed a dynamic regulation of gene expression levels within these families in response to increasing duration of low-temperature treatment, while the expression levels of certain transcription factors decreased with increasing low-temperature treatment time. These findings suggest that the bHLH and MYB gene family transcription factors might exert either a positive or negative regulation on the cold tolerance of common beans. However, further validation is warranted to corroborate these observations. The effect of low-temperature stress on sugar acid metabolism and amino acid metabolism in common beans Research has demonstrated the crucial role of soluble sugars, such as galactose, starch, pyruvate, and amino acids, in augmenting plant tolerance to abiotic stress . A previous study found that the levels of soluble starch, fructose, glucose, and phenolic substances in citrus leaves undergo significant elevation after cold stress, underscoring the potential importance of sugar and secondary metabolism pathways in citrus responses to cold stress . Research by Zhao et al. suggests that overexpression of β- amylase (PbrBAM3) in Pyrus betulaefolia facilitates starch degradation under cold stress, thereby enhancing cold tolerance . In our study, we observed a significant reduction in the levels of galactinol, sucrose, melibiose, trehalose, maltose, and fructose-6P in common beans under cold stress (Fig. ). This finding suggests that common beans may enhance their cold tolerance by modulating the levels of sugar and acid substances. Furthermore, amino acids are crucial for protein synthesis, and previous studies have highlighted the pivotal role of amino acid metabolism in enhancing plant abiotic stress tolerance . Several amino acids, including proline, arginine, asparagine, glutamine, and GABA, are synthesized at high abundance under abiotic stress conditions, serving various functions such as maintaining compatible osmotic pressure, acting as precursors of secondary metabolites or serving as storage forms of organic nitrogen . Among them, proline is an osmotic regulator and is crucial in mediating plant responses to abiotic stress, especially extreme temperature stress . Our findings indicate that cold stress significantly elevates the content of Proline, L-Leucine, L-Valine, Tyrosine, D-Aspartic Acid, and L-Aspartate in bean leaves, suggesting that beans may alleviate cold stress through regulation of amino acid content. Effects of cold stress on secondary metabolite flavonoids in common beans Flavonoids play a significant role in the secondary metabolism of plants, particularly in their response to cold stress . Our study identified 19 DEMs and 15 DEGs associated with the biosynthesis and metabolism of flavonoids in common beans. Transcriptomic and metabolomic analysis of two peach trees with differing cold tolerance revealed that certain secondary metabolites, including phenolic acids and flavonoids, were exclusively upregulated in cold-sensitive peach trees . In Galega officinalis , cold stress enhances the production of phenolic compounds, including flavonoids, apigenin, coumaric acid, genistein, luteolin, trans ferulic acid, and naringenin . Combined transcriptomic and metabolomic analysis of Fagopyrum tataricum at different altitudes revealed that cold stress significantly upregulated phenylpropanoid biosynthesis and promoted the expression of anthocyanins . Plant isoflavones are naturally occurring plant estrogens belonging to the flavonoid class and cannot be synthesized in the human body . A previous study found that low-temperature treatment significantly increased the content of phenolic acids and isoflavones (genistein, daidzein, and genistein) in soybean roots, with the largest increase observed in genistein after 24 h of treatment at 10 ° C . Our findings indicate that low-temperature stress significantly increases the content of isoflavone metabolites, such as genistein, and isoflavone synthesis precursors, including L-phenylalanine, naringenin chalcone, and naringenin. The expression trends in cold-tolerant WY and cold-sensitive BBL common beans are consistent (Fig. 10), suggesting that common beans may alleviate low-temperature stress on plants by upregulating isoflavone substances. The effect of cold stress on plant hormones in common beans Furthermore, plant hormones are pivotal in co-regulating plant secondary metabolism and enhancing resistance against abiotic stress-induced damage in plants . The accumulation of anthocyanins, a secondary metabolite in plants, is influenced by various abiotic stresses, such as high light intensity, cold, drought, salinity, nutrient deficiency, and heavy metal stress, as well as the induction of endogenous plant hormones . By studying the adaptation mechanism of corn to cold stress, it was found that gibberellins play a regulatory role in the accumulation of anthocyanins induced by low temperatures . During cold stress, transgenic rice lines overexpressing OsABA8ox1 , a gene involved in abscisic acid synthesis, exhibited reduced ABA content and enhanced seedling vitality. This suggests that maintaining low levels of ABA during cold stress can promote seedling vigor . In Arabidopsis, treatment with epibrassinolide solution (EBR)enhances the tolerance of seedlings to both drought and cold stress . Research has also revealed that rice-specific microRNA miR1320 targets the ethylene-responsive transcription factor OsERF096, thereby regulating cold stress tolerance by suppressing the JA-mediated cold signaling pathway . This study identified three DEMs and 18 DEGs involved in the biosynthesis pathways of BR (brassinolide), ABA (abscisic acid), JA (jasmonic acid), and SA (salicylic acid) in common beans. We also found that the content of SA in cold-sensitive common bean BBL increased with prolonged exposure to low temperatures. Conversely, in the cold-resistant WY common bean variety, SA content exhibited an initial increase followed by a decrease with prolonged low-temperature treatment (Fig. A). These findings suggest that common beans may alleviate oxidative damage induced by low temperatures by regulating endogenous SA levels. However, variations exist among varieties with differing levels of cold tolerance. 

Low-temperature stress has notable impacts on the osmotic regulation system and antioxidant system of plants. In Brassica napus L. (rapeseed), low-temperature stress significantly elevated the levels and relative conductivity of malondialdehyde, alongside enhancing the activity of peroxidase and superoxide dismutase . Our research findings demonstrate that low-temperature stress significantly elevated the malondialdehyde content, relative conductivity, and antioxidant enzyme activity in two distinct varieties of common beans, and the observed increments varied between the two types of common beans. These findings are consistent with those of Cai et al. in two varieties of Solanum melongena exhibiting differential cold tolerance and Wang et al. in Brassica campestris L. with varying cold tolerance levels . Together, they underscore that distinct varieties within the same species can manifest divergent responses to low-temperature stress. 

Currently, a plethora of plant transcription factors have been shown to respond to low-temperature stress in plants, primarily modulating plant cold resistance through involvement in processes such as plant cell membrane fluidity, cold signal transduction, MAPK cascade, and regulation of CBF pathways . Studies have revealed that transgenic apple callus and Arabidopsis plants overexpressing MdMYB23 exhibit heightened cold tolerance . MdMYB23 has been shown to interact with the key regulator of anthocyanin biosynthesis, the MdANR promoter, thereby activating its expression, promoting anthocyanin accumulation, and clearance of reactive oxygen species (ROS) . In addition, plant alkaline helix loop helix (bHLH) transcription factors play pivotal roles in plant growth and development, secondary metabolism, and responses to abiotic stress . In rice, the upregulation of OsbHLH1 is specifically triggered by cold stress, suggesting its involvement in the cold signaling pathway of rice . In this study, a comprehensive analysis identified a total of 71 bHLH and 64 MYB gene family transcription factors (Fig. ). Interestingly, we observed a dynamic regulation of gene expression levels within these families in response to increasing duration of low-temperature treatment, while the expression levels of certain transcription factors decreased with increasing low-temperature treatment time. These findings suggest that the bHLH and MYB gene family transcription factors might exert either a positive or negative regulation on the cold tolerance of common beans. However, further validation is warranted to corroborate these observations. 

Research has demonstrated the crucial role of soluble sugars, such as galactose, starch, pyruvate, and amino acids, in augmenting plant tolerance to abiotic stress . A previous study found that the levels of soluble starch, fructose, glucose, and phenolic substances in citrus leaves undergo significant elevation after cold stress, underscoring the potential importance of sugar and secondary metabolism pathways in citrus responses to cold stress . Research by Zhao et al. suggests that overexpression of β- amylase (PbrBAM3) in Pyrus betulaefolia facilitates starch degradation under cold stress, thereby enhancing cold tolerance . In our study, we observed a significant reduction in the levels of galactinol, sucrose, melibiose, trehalose, maltose, and fructose-6P in common beans under cold stress (Fig. ). This finding suggests that common beans may enhance their cold tolerance by modulating the levels of sugar and acid substances. Furthermore, amino acids are crucial for protein synthesis, and previous studies have highlighted the pivotal role of amino acid metabolism in enhancing plant abiotic stress tolerance . Several amino acids, including proline, arginine, asparagine, glutamine, and GABA, are synthesized at high abundance under abiotic stress conditions, serving various functions such as maintaining compatible osmotic pressure, acting as precursors of secondary metabolites or serving as storage forms of organic nitrogen . Among them, proline is an osmotic regulator and is crucial in mediating plant responses to abiotic stress, especially extreme temperature stress . Our findings indicate that cold stress significantly elevates the content of Proline, L-Leucine, L-Valine, Tyrosine, D-Aspartic Acid, and L-Aspartate in bean leaves, suggesting that beans may alleviate cold stress through regulation of amino acid content. 

Flavonoids play a significant role in the secondary metabolism of plants, particularly in their response to cold stress . Our study identified 19 DEMs and 15 DEGs associated with the biosynthesis and metabolism of flavonoids in common beans. Transcriptomic and metabolomic analysis of two peach trees with differing cold tolerance revealed that certain secondary metabolites, including phenolic acids and flavonoids, were exclusively upregulated in cold-sensitive peach trees . In Galega officinalis , cold stress enhances the production of phenolic compounds, including flavonoids, apigenin, coumaric acid, genistein, luteolin, trans ferulic acid, and naringenin . Combined transcriptomic and metabolomic analysis of Fagopyrum tataricum at different altitudes revealed that cold stress significantly upregulated phenylpropanoid biosynthesis and promoted the expression of anthocyanins . Plant isoflavones are naturally occurring plant estrogens belonging to the flavonoid class and cannot be synthesized in the human body . A previous study found that low-temperature treatment significantly increased the content of phenolic acids and isoflavones (genistein, daidzein, and genistein) in soybean roots, with the largest increase observed in genistein after 24 h of treatment at 10 ° C . Our findings indicate that low-temperature stress significantly increases the content of isoflavone metabolites, such as genistein, and isoflavone synthesis precursors, including L-phenylalanine, naringenin chalcone, and naringenin. The expression trends in cold-tolerant WY and cold-sensitive BBL common beans are consistent (Fig. 10), suggesting that common beans may alleviate low-temperature stress on plants by upregulating isoflavone substances. 

Furthermore, plant hormones are pivotal in co-regulating plant secondary metabolism and enhancing resistance against abiotic stress-induced damage in plants . The accumulation of anthocyanins, a secondary metabolite in plants, is influenced by various abiotic stresses, such as high light intensity, cold, drought, salinity, nutrient deficiency, and heavy metal stress, as well as the induction of endogenous plant hormones . By studying the adaptation mechanism of corn to cold stress, it was found that gibberellins play a regulatory role in the accumulation of anthocyanins induced by low temperatures . During cold stress, transgenic rice lines overexpressing OsABA8ox1 , a gene involved in abscisic acid synthesis, exhibited reduced ABA content and enhanced seedling vitality. This suggests that maintaining low levels of ABA during cold stress can promote seedling vigor . In Arabidopsis, treatment with epibrassinolide solution (EBR)enhances the tolerance of seedlings to both drought and cold stress . Research has also revealed that rice-specific microRNA miR1320 targets the ethylene-responsive transcription factor OsERF096, thereby regulating cold stress tolerance by suppressing the JA-mediated cold signaling pathway . This study identified three DEMs and 18 DEGs involved in the biosynthesis pathways of BR (brassinolide), ABA (abscisic acid), JA (jasmonic acid), and SA (salicylic acid) in common beans. We also found that the content of SA in cold-sensitive common bean BBL increased with prolonged exposure to low temperatures. Conversely, in the cold-resistant WY common bean variety, SA content exhibited an initial increase followed by a decrease with prolonged low-temperature treatment (Fig. A). These findings suggest that common beans may alleviate oxidative damage induced by low temperatures by regulating endogenous SA levels. However, variations exist among varieties with differing levels of cold tolerance. 

This study elucidated the impact of cold stress on physiological parameters, gene expression, and metabolite profiles of common bean seedlings through an integrated analysis encompassing phenotypic physiology, transcriptomics, and metabolomics. The results suggest that under cold stress, the DEGs and DEMs of common beans are engaged in primary metabolism, secondary metabolism, and plant hormone signal transduction, especially in synthesizing secondary metabolites such as isoflavones. Furthermore, it was found that bHLH and MYB transcription factors exhibit extensive involvement in the cold response of common beans. In summary, this study offers valuable insights for a more comprehensive understanding of the cold resistance mechanism in common beans and the exploration of cold-tolerant bean germplasm resources. 

Below is the link to the electronic supplementary material. Supplementary Material 1 

