Genome-wide identification and expression profiling of cowpea HSP60 genes reveal key candidates for temperature stress responses

Cowpea (Vigna unguiculata) is an annual herbaceous plant of the Fabaceae family. This crop is cultivated around the world for its seeds and immature pods, which are rich in protein, as well as for its economic value^[1]. Cowpea cultivation contributes enormously to sustainable agriculture and food security for millions of smallholder farmers in tropical and subtropical areas^[2]. However, the escalating volatility of global weather patterns has subjected cowpea to extreme temperature fluctuations that exceed its adaptive thresholds^[3]. Heat and cold shocks now represent major abiotic barriers to phytogeography, growth, and yield and are directly threatening the stability of global food production amid climate changes^[4]. When a plant is exposed to inadequate temperatures, ionic balance is disturbed, membrane stability is compromised, proteins are denatured, and ultimately yield and seed quality are reduced^[5,6]. Therefore, increasing cowpea's resistance to temperature stress offers both strategic and practical agricultural benefits for the production of climate-resilient crops.

Plants suffer irreversible damage from extreme temperatures during their life cycle. Thus, plants have established a range of physiological and biological defense mechanisms to cope with stressful environments^[7]. At the molecular level, multiple stress proteins are induced, such as heat shock proteins (HSPs). It has been reported that HSPs act as molecular chaperones to withstand the detrimental effects of thermal stress. HSPs are categorized according to their molecular weight into five major families: HSP100, HSP90, HSP70, HSP60, and small HSPs (sHSPs)^[8]. Among these, proteins in the HSP60 family, also known as "chaperonins" (Cpn60), are crucial for protein folding, refolding, managing aggregation, and trafficking to mitochondria and chloroplasts^[9]. HSP60 chaperones are essential for organelle function, photosynthetic efficiency, development, and overall thermotolerance in plants. They are widely found in eukaryotes, eubacteria, and archaea. In plants, HSP60s are generally classified into Group I chaperonins, which resemble the bacterial Growth Relation Large (GroEL) system and act with co-chaperonins (Cpn10 or Cpn20)^[10]. HSP60 proteins are divided into α (alpha) and β (beta) subunits in plant cells. They assemble into oligomeric complexes that facilitate adenosine triphosphate (ATP)-dependent protein folding. Another form, sometimes referred to as the CCT-type chaperonins, belongs to the Group II chaperonins, mainly found in the eukaryotic cytosol and involved in folding cytoskeletal and regulatory proteins^[11]. These classes ensure adequate protein homeostasis under stress conditions.

Recent studies showed that HSP60s are a key regulator of plant thermotolerance, working through several molecular processes. Studies on their function have shown that knocking down HSP60 significantly reduces the plant's ability to withstand heat stress. The silencing of CaHSP60-6 in pepper (Capsicum annuum) led to an increase in sensitivity to heat stress, as evidenced by increased electrolyte leakage, lipid peroxidation, reactive oxygen species accumulation, reduced chlorophyll content, and diminished antioxidant enzyme activity. In contrast, heat treatment resulted in an upregulation of 93% of the HSP60 genes in both thermosensitive and thermotolerant pepper lines^[12]. In rice (Oryza sativa), OsHSP60-3B is essential for the proper formation of starch granules within pollen, thus safeguarding male fertility under elevated temperatures. The OsHSP60-3B mutant displayed typical fertility under ideal temperature conditions; however, its fertility diminished with rising temperatures, leading to impaired pollen starch accumulation and altered pollen morphology. In plastids, OsHSP60-3B interacts with FLOURY ENDOSPERM6 (FLO6) in a mechanism that stabilizes this essential enzyme for starch biosynthesis and reduces the levels of reactive oxygen species in the anthers^[13]. Comprehensive genome-wide analyses have demonstrated that HSP60 genes are subject to extensive transcriptional regulation in response to temperature stress. Multiple CsHSP60 genes in tea (Camellia sinensis) plants exhibited substantial upregulation in response to drought and heat^[14]. Chloroplast-localized HSP60 proteins, by interacting with ribulose-1, 5-bisphosphate carboxylase-oxygenase (Rubisco) activase, might mitigate its thermal inactivation, thus protecting the photosynthetic apparatus under heat stress^[15]. Levels of plastid Cpn60 increase in response to stress during seed germination and early development^[16]. Moreover, proteomic analyses of soybean (Glycine max) and azalea (Rhododendron spp.) identified HSP60 as part of wide heat-primed molecular networks that improve thermotolerance via pre-activated transcriptional programs^[17,18]. Overall, these studies demonstrate the critical function that HSP60 proteins play in various abiotic stress responses.

Cowpea, an important warm-season crop, is sensitive to temperature stress, particularly during germination and flowering^[19]. In subtropical zones, lower temperatures (below 15 °C) can inhibit or reduce germination rates by triggering imbibitional chilling injury and severe membrane fracturing. Following emergence, seedlings suffer from a rapid accumulation of toxic reactive oxygen species, which disrupt cell metabolism and lead to phenotypic wilting, chlorophyll degradation, and leaf chlorosis. Conversely, high daytime temperatures (above 35 °C) coupled with long photoperiods abort early floral buds, whereas elevated nocturnal temperatures 7–9 d before anthesis induce male sterility by impairing pollen viability, thereby preventing successful fertilization. This heat stress ultimately suppresses pod development and causes near-total fruit-setting failure, severely limiting subtropical yields^[20,21]. Despite the well-known importance of these genes for alleviating plant environmental stress, a comprehensive assessment of the HSP60 family in cowpea has not been reported to date. Identifying these genes and their expression profiles is therefore essential to pinpointing the candidates responsible for thermal adaptation. In this work, we present a thorough analysis of the HSP60 gene family in cowpea under heat and cold stress. Gene structure, conserved domains and motifs, chromosomal distributions, phylogenetic analysis, cis-acting elements in the 5' upstream region, and expression patterns were examined. The findings provide a vital baseline for functional validation and genetic enhancement of cowpea cultivars that are resistant to temperature stress.

Genomes of cowpeas (V. unguiculata v1.2) and Arabidopsis thaliana (Araport11) were obtained from the Phytozome13 database (https://phytozome-next.jgi.doe.gov). A. thaliana HSP60 proteins were taken from the NCBI database (www.ncbi.nlm.nih.gov) and used as queries in TBtools' BLAST to find their protein homologs in the V. unguiculata genome. The Pfam database (http://pfam.xfam.org) provided the hidden Markov model (HMM) profile of the domain of each family in A. thaliana. We identified the domains in the cowpea proteome using HAMMER 3.0. After duplicates were removed, sequences lacking the Cpn60_TCP1 domain were also excluded using NCBI's Conserved Domain Database (CDD) (www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART (http://smart.embl-heidelberg.de)^[22,23]. In the HSP60 family, the genes encoding a protein with the required domain were taken as candidates for further analysis. Hydropathy, instability index, molecular weight (MW), sequence size, and isoelectric point (pI) were among the physical and chemical properties of the chosen candidate genes, as determined by TBtools II. To determine the subcellular location of proteins inside eukaryotic sequences, we used the CELLO subcellular localization system (https://cello.life.nctu.edu.tw)^[24]. According to their chromosomal positions and distributions, these genes were systematically designated using the prefix 'Vu' for V. unguiculata.

Chromosome arrangement and gene duplication analysis

Using TBtools II software, we evaluated HSP60 proteins from the Phytozome v13 database to establish their distribution on different chromosomes and to calculate the ratio of nonsynonymous to synonymous substitutions (Ka/Ks) for the gene duplication events. The Nei–Gojobori approach was used to remove any blank or missing values from the analysis^[25].

Evolutionary relationship and sequence alignment

The phylogenetic tree was constructed using sequences of HSP60 proteins from three different species (V. unguiculata, A. thaliana, G. max, and O. sativa). A neighbor-joining (NJ) tree with a bootstrap value of 1,000 repetitions was constructed using the Molecular Evolutionary Genetic Analysis (MEGA12) program. To show the evolutionary tree, the Interactive Tree of Life (ITOL; https://itol.embl.de) was used^[26]. ClustalW was used to align the HSP60 protein sequences in MEGA12.

Analysis of conserved motifs, gene structure, and miRNA prediction

In this study, the VuHSP60 protein sequences were used to assess the exon/intron architecture using the TBtools II software. NCBI's CDD (www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) was used to identify the conserved domains. MEME in TBtools was used to identify conserved protein motifs^[27]. In addition, we used the Plant miRNA Encyclopedia (PmiREN; https://pmiren.com) and Plant Small RNA Target Analysis Server (psRNATarget; www.zhaolab.org/psRNATarget) databases to predict miRNAs associated with VuHSP60 proteins. The miRNA–target gene network diagram was generated by www.bioinformatics.com.cn, an online platform for data analysis and visualization^[28].

Examination of synteny, collinearity, and cis-acting elements

For synteny and collinearity studies, orthologs of the cowpea HSP60 genes were identified across four species (namely, A. thaliana, G. max, Zea mays, and Medicago truncatula) using genomic data retrieved from the Phytozome v13 database. TBtools II's one-step MCScanX method was used to examine the orthologs, and the interspecific collinearity map was displayed using the Multiple Synteny Plot option^[29]. The 2,000-bp promoter sequences of the cowpea VuHSP60 gene family were investigated using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html) to predict cis-acting regulatory elements^[30]. Data visualization was facilitated with TBtools II software.

Plant materials and temperature stress treatments

Low- and high-temperature-sensitive cowpea ('BBWQ') seeds were provided by Nanjing Agricultural University's Brassica System Biology Laboratory. These seeds were planted in pots with regular potting mix and grown in a plant growth chamber with the ideal controlled conditions: 25 °C, 70% relative humidity, and a 16-h light and 8-h dark photoperiod. At exactly 2 weeks after the sowing date, when the seedlings reached the early trifoliate stage, uniform plants were selected for stress treatments and while some seedlings were maintained under ideal conditions to serve as controls. For cold stress, the seedlings were kept at 10 °C for a full day with a cycle of 16 h of light and 8 h of darkness, whereas heat stress was induced by the exposure of seedlings to 40 °C during 24 h. Young leaf samples were harvested at 1, 3, 6, 12, and 24 h after the onset of treatment, with three distinct biological replicates at each time point. Leaves of the control plants were also collected. In preparation for the subsequent RNA extraction and analysis, all samples were promptly frozen in liquid nitrogen and kept at −80 °C.

Quantitative real-time polymerase chain reaction

Total RNA was extracted from the samples according to the manufacturer's instructions using the Total RNA Extraction Reagent (Vazyme Biotech Co., Ltd., Nanjing, China). A NanoDrop™ One spectrophotometer measured the RNA's purity and concentration at 260 nm. RNA integrity was further confirmed using 1.0% agarose gel electrophoresis. We used the HiScript® II 1st Strand cDNA Synthesis Kit (Vazyme Biotech Co., Ltd.) to generate the first strand of cDNA. The study used a Bio-Rad CFX96 Touch Real-Time PCR Detection System and 2× ChamQ Blue Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China). Each 20 μL reaction mixture included 10 μL of 2× SYBR Master Mix, double-distilled water (ddH₂O), 2.0 μL of the diluted cDNA template (1:10), and 0.4 μL of the primers. The thermal cycling conditions included an initial denaturation at 95 °C for 4 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. The 2^−ΔΔCᴛ technique was used to calculate the relative gene expression levels. The internal reference was the cowpea actin gene (Accession No. XM_028073549). Primer specificity was confirmed with NCBI Primer-BLAST, and primer pairs were generated using the Primer3 online tool based on the coding sequences (CDS) of cowpea HSP60 genes (Supplementary Table S1). Supplementary Table S2 shows the primer sequences used in this study. Each quantitative real-time polymerase chain reaction (qRT-PCR) test consisted of three biological duplicates and three technical repetitions each^[31,32].

Expression data analysis and statistical methods

Analysis of variance (ANOVA) in SPSS version 23.0 was used for statistical analysis. Significant differences between the control and treatment groups, as well as at different time points, were established using Tukey's test at the 5% significance level (p < 0.05). Three biological replicates were used in each experiment, and each treatment's standard error of the mean (SE) was determined. The data were visualized using GraphPad Prism 10^[33].

In total, 25 HSP60 genes were identified in the current study through genome-wide scanning of V. unguiculata. These proteins were dispersed randomly over nine distinct chromosomes. At minimum, three genes were found on chromosomes 3, 7, and 11, among other chromosomes. Other chromosomes contain one or two genes of HSP60 (Fig. 1). However, VuHSP60-18 was not assigned to a specific chromosome and remains on an unplaced contig (Contig_339). According to their chromosomal positions and distributions, these V. unguiculata HSP60 genes were named with the prefix 'Vu'. The VuHsp60 genes are fully described in Table 1 below. The HSP60s ranged in molecular weight from 44.91 kDa (VuHSP60-15) to 71.76 kDa (VuHSP60-18) and in amino acid content from 414 (VuHSP60-15) to 680 (VuHSP60-18). The pIs ranged from 5.15 (VuHSP60-17) to 9.11 (VuHSP60-9.1). Of the 25 HSP60s, five were identified as basic (pI > 7), while the remaining 20 were classified as acidic (pI < 7). According to the grand average of hydropathy (GRAVY), the majority of VuHSP60s have negative GRAVY values, which suggest that they are hydrophilic, with the exception of VuHSP60-3, VuHSP60-5.1, VuHSP60-5.2, VuHSP60-5.3, and VuHSP60-10, which are hydrophobic, with GRAVY > 0. Moreover, 24 of the 25 VuHSP60s are stable, with instability indices ranging from 23.64 (VuHSP60-14.1) to 39.99 (VuHSP60-2), but only one protein is unstable with a value of 46.12 (VuHSP60-11.1). The aliphatic index of the majority of HSPs in cowpea ranges from 96.04 (VuHSP60-13) to 115.78 (VuHSP60-3). Subcellular localization predicted that VuHSP60s are mainly identified in the cytoplasm, mitochondria, and chloroplasts. The highest number of HSP60s are localized in the cytoplasmic compartment, with 13 genes, followed by the mitochondria, where 5 HSP60 genes were counted (Table 1).

Figure 1. 

VuHSP60 genes' chromosomal positions.

Gene duplication analysis showed that there were 93 paralogous pairs in the VuHSP60 gene family (Supplementary Table S3). When genes are duplicated and positioned on distinct chromosomes, this scenario is known as segmental duplication. Segmental duplication was the most prevalent kind in VuHSP60, with 91 gene pairs. In VuHSP60, dispersed duplication was less common, occurring in two gene pairs. These results suggest that the evolution of this gene family in cowpea was significantly influenced by large-scale chromosomal duplication events. The predicted Ka, Ks, and Ka/Ks values for these gene pairs are shown in Supplementary Table S3. Since all Ka/Ks ratios are less than 1, these genes likely underwent strong purifying selection during the process of evolution. This suggests functional conservation and selective pressure to maintain protein stability and essential biological functions.

Figure 1 indicates the localization of HSP60 gene family members. Megabases (Mb) indicate chromosomal size and the genes' locations. For loci exhibiting alternative splicing events, the base gene name represents the physical locus position, and only the primary transcript model coordinates were utilized for visualization of the map layout. The chromosome numbers of V. unguiculata are numbered from Vu01 to Vu10 and are indicated at the left of each bar. The arrangement of genes in cowpea chromosomes was developed by the TBtools II program.

Evolutionary relationships of HSP60s in V. unguiculata

We used the neighbor-joining method in MEGA 12 software to generate a rooted phylogenetic tree (Fig. 2). The full-length HSP60 protein sequences from three dicots (cowpea, A. thaliana, and soybean) and one monocot (rice) were used (Supplementary Table S4). According to their localization, this study categorized 98 HSP60 proteins into four groups (Group I, II, III, and IV). Group I, which is the biggest, contains 30 cytoplasmic proteins. Of these, 13 (52% of the total VuHSP60s) come from cowpea, nine come from Arabidopsis, eight come from rice, and 18 from soybean. Groups II and III comprise mitochondrial and chloroplastic proteins, respectively; cowpea genes represent 28% of VuHSP60s in Group II and 20% in Group III and contain representatives from four plant species. Group IV is completely devoid of both cowpea (V. unguiculata) and soybean genes, containing exclusively orthologs from A. thaliana and rice. Furthermore, each group includes a mix of monocot and dicot proteins.

Figure 2. 

Phylogenetic trees of the VuHSP60 gene family.

MEGA12 software was used to generate the phylogenetic trees using the HSP60 protein sequences from cowpea, A. thaliana, soybean, and rice. The figure shows the HSP60 genes' phylogenetic tree, which is divided into four groups (I, II, III, and IV) of 98 genes. The Interactive Tree of Life (iTOL) was used to illustrate the figure, with bootstrap values derived from 1,000 replicates.

Gene structure and conserved motif analysis

We investigated the structural variation of the VuHSP60 genes by analyzing the exon–intron organization, domains, and motifs (Fig. 3). Closely related members within the same phylogenetic group exhibited similar gene architectures. Intron counts varied widely across the family, ranging from 0 to 17. The highest intron numbers (16–17) were observed in VuHSP60-4, VuHSP60-14.1, VuHSP60-14.2, VuHSP60-10, and VuHSP60-12, which encode proteins predicted by CELLO to be localized to the cytoplasm and mitochondria. VuHSP60-3 differs from the other VuHSP60 members in that it lacks introns. According to the results, the number of introns of the remaining VuHSP60 genes range from 5 to 14. To further characterize the proteins, we used the MEME suite (integrated in TBtools) to identify the conserved motifs shared by related proteins in each subfamily. We found 10 potential motifs of 21–50 amino acids across the HSP60 family (Supplementary Table S5). The arrangement and composition of motifs vary throughout the family and are correlated with phylogenetic clustering. Motifs 1 and 2, which are associated with the Cpn60_TCP1 domain, are present in all members of the VuHSP60 family. Conserved domain annotation (Pfam, CDD, SMART) databases revealed a highly conserved domain architecture, confirming their classification as molecular chaperones.

Figure 3. 

Conserved motifs, gene structure, conserved domains, and phylogenetic relationship of VuHSP60 proteins.

The left panel (a) displays the phylogenetic tree of the three VuHSP60 protein family. The full-length amino acid sequences of VuHSP60 proteins were used to construct the rooted neighbor-joining phylogenetic tree with MEGA 12 (Supplementary Table S6). The bootstrap test replicates were set to 1,000. The second panel (b) illustrates the distributions of conserved motif proteins. The 10 possible motifs are marked by numbered frames in various colors, ranging from 1 to 10. Panel (c) represents the conserved domains that are present in VuHSP60 proteins. The right panel (d) illustrates the exon–intron layout of the VuHSP60 genes. Gray lines symbolize introns, and the yellow rectangles symbolize exons (CDS). Green boxes represent the upstream and downstream regions (untranslated regions, UTRs) of VuHSP60 genes. The graphic also shows a description of each colored box.

In silico prediction of microRNAs potentially targeting the VuHSP60 gene family

In silico prediction identified a total of 20 miRNAs predicted to target and regulate VuHSP60 genes (Fig. 4). These miRNAs were classified into eight families involved in RNA cleavage and translation inhibition. Important miRNAs known to regulate VuHSP60s include Vun-miR3309, Vun-miR165, Vun-miR2111, Vun-miR394, and Vun-miR408. Vun-miR167 is involved in plant growth and targets VuHSP60-5.1, VuHSP60-5.2, and VuHSP60-5.3. The stress-responsive Vun-miR164 targets VuHSP60-9.1 and VuHSP60-9.2. A multigene targeting network was observed, particularly for the plant growth regulator Vun-miR396, which participates in stress responses to drought, high salinity, and extreme temperatures and targets five specific genes: VuHSP60-8, VuHSP60-15, VuHSP60-16, VuHSP60-14.1, and VuHSP60-14.2.

Figure 4. 

Predicted miRNA–VuHSP60 gene interaction network. Light blue circles represent VuHSP60 genes, yellow triangles represent miRNAs, and gray lines show predicted targeting interactions.

Synteny and collinearity investigation of the HSP60 genes from V. unguiculata and other species

To infer evolutionary conservation and genomic divergence, we used TBtools to perform synteny and collinearity analyses of the HSP60 gene family in V. unguiculata with A. thaliana, G. max, Z. mays, and M. truncatula (Fig. 5). Syntenic blocks between V. unguiculata and three of the four comparison species were detected. Specifically, G. max and V. unguiculata showed 27 collinear gene pairs within the syntenic blocks, indicating strong microsynteny conservation despite legumes' diversification. In contrast, V. unguiculata and A. thaliana showed nine collinear gene pairs, whereas M. truncatula showed only eight. These results demonstrate that genomic conservation is strictly higher within the tribe Phaseoleae and substantially reduced in both distant legumes and nonlegume dicots. Surprisingly, no collinear gene pairs were found between Z. mays and V. unguiculata, reflecting genomic alterations in specific lineages, along with the deep evolutionary divergence between dicots and monocots.

Figure 5. 

Synteny and collinearity analysis. Blue lines indicate the collinear HSP60 gene pairs within conserved syntenic blocks, while the faint gray lines in the background represent the genome-wide synteny (all other collinear gene pairs) in five species.

Cis-acting regulatory element analysis

To predict potential cis-regulatory elements, we uploaded 2,000-bp upstream promoter regions from VuHSP60 genes to the PlantCARE database. The results revealed 37 distinct types of cis-regulatory elements (Fig. 6). These components participate in multiple biological processes, such as development, stress response, hormone signaling, light perception, and promoter activity. Of the 25 genes, 12 (48%) contained low-temperature responsive elements (LTRs), and 9 (36%) contained defense and stress responsive elements (TC-rich repeats). In terms of hormone-responsive elements, 17 (68%), 7 (28%), 11 (44%), 6 (24%), and 18 (72%) genes contained one or more abscisic acid-responsive elements (ABRE), auxin (indole-3-acetic acid; IAA)-responsive elements, gibberellin (GA)-responsive elements, salicylic acid (SA)-responsive elements, and methyl jasmonate (MeJA)-responsive elements. Other elements detected in 20 and two VuHSP60 genes were anaerobic response elements (AREs) and wound-response elements (WUN motifs), respectively. The existence of LTR elements and hormones responding to environmental stresses points to the predicted role for the VuHSP60 genes in mediating cowpea's adaptation to extreme temperatures.

Figure 6. 

Cis-elements predicted in the promoter regions of VuHsp60 genes.

Cis-regulatory elements were predicted by PlantCARE in the promoter sequences (2,000 bp) of the 25 VuHSP60 genes that were analyzed. The various types of cis-acting elements are represented by distinct colored shapes. We used the TBtools software to visualize all components.

Relative expression levels of VuHSP60s under temperature stress

We investigated the expression levels of 25 members of the HSP60 gene family in response to heat and cold stressors (Fig. 7). The expression trends of nine genes, including VuHSP60-1, VuHSP60-3, VuHSP60-5.1, VuHSP60-5.2, VuHSP60-5.3, VuHSP60-11.2, VuHSP60-12, VuHSP60-14.1, and VuHSP60-17, exhibited a significant increase under high-temperature treatment in comparison with the control. VuHSP60-12 and VuHSP60-14.1 exhibited the highest expression levels among these genes, peaking at 1 h (6.64-fold) and 12 h (4-fold) after treatment, respectively, before gradually declining. The expression trends of VuHSP60-5.1, VuHSP60-5.2, and VuHSP60-5.3 increased gradually and reached their highest point at 24 h post-treatment. VuHSP60-1, VuHSP60-3, VuHSP60-11.2, and VuHSP60-17 showed different variations. Conversely, the remaining 16 family members displayed either nonsignificant changes or significant downregulation relative to the control. Under cold treatment, 10 genes (VuHSP60-1, VuHSP60-4, VuHSP60-7, VuHSP60-8, VuHSP60-9.1, VuHSP60-9.2, VuHSP60-12, VuHSP60-13, VuHSP60-15, and VuHSP60-16) were significantly upregulated. Notably, VuHSP60-4, VuHSP60-7, and VuHSP60-15 exhibited the highest expression levels, reaching 8-, 7.5-, and 7.2-fold increases over the control, respectively. The remaining 15 VuHSP60 genes showed differential expression levels with nonsignificant increases. In general, these distinct temporal expression profiles demonstrate that specific VuHSP60 genes actively participate in heat, cold, or dual thermal stress responses, indicating functional diversification within the gene family.

Figure 7. 

Expression levels of 25 VuHSP60 genes.

Figure 7 illustrates the qRT-PCR results of the leaf tissues of the 'BBWQ' cowpea cultivar under heat stress (40 °C) and cold stress (10 °C). Panels are categorized into Group I, II, and III according to the phylogenetic classification of the VuHSP60 gene family. The 2^−ΔΔCᴛ method served to determine the expression levels; the bars represent the mean ± standard error. The Tukey test in SPSS 23.0 demonstrated statistically significant differences between the control and treatments. Asterisks indicate statistically significant differences compared with the control (CK) (* p < 0.05, ** p < 0.01, *** p < 0.001, and **** p < 0.0001), whereas the absence of asterisks shows no significant difference (p ≥ 0.05). The graphs were generated using GraphPad Prism 10.

Cowpea (V. unguiculata) is a globally important crop^[34], but its yield is severely constrained by specific phenological vulnerabilities, including cold-induced chlorosis during early vegetative establishment and acute heat-induced pollen sterility during blooming in subtropical regions^[35,36]. These thermal shocks cause protein denaturation and metabolic disruption of important chloroplast and mitochondrial bioenergetic pathways^[37]. According to earlier studies, HSP60s are essential for how plants respond to temperature stress^[38]. Therefore, this gene family has been studied in the genomes of several species. However, very little is known about HSP60 genes in cowpea. In this study, 25 HSP60 genes were detected in a systematic genome-wide analysis of V. unguiculata (Table 1), a number similar to those found in multiple plant species such as O. sativa (22)^[39], mung bean (Vigna radiata) (23)^[40], Sorghum bicolor (30)^[41], A. thaliana (29)^[42], C. sinensis (19)^[14], and Jerusalem artichoke (Helianthus tuberosus) (22)^[43]. The VuHSP60 genes in this study have mostly experienced segmental duplication (97.8%) and minor dispersed duplication (2.2%) (Supplementary Table S3), indicating that the expansion of this gene family in V. unguiculata was primarily driven by large-scale genomic events. These findings are consistent with those of other crops where segmental duplication has played a significant role in the spread of the HSP60 gene family, including tea^[14], S bicolor^[41], and soybean^[44]. Additionally, the Ka/Ks ratios of all duplicated VuHSP60 gene pairs were less than 1, indicating that this family has evolved under purifying selection that maintains a functional conservation among the paralogs. The physicochemical characteristics of proteins, such as their MW and pI, provide important information about their functions^[45]. In the present study, the VuHSP60 proteins exhibited diverse MW and pI values (Table 1). Most of them possessed high aliphatic indices, indicating enhanced thermostability, alongside hydrophilic and acidic properties^[46]. In molecular chaperones, these biochemical properties are typically associated with increased protein solubility and the prevention of thermal aggregation within aqueous cellular environments^[47]. Therefore, various profiles of VuHSP60 proteins may reflect an evolutionary adaptation to preserve organellar structural integrity during the severe thermal variation characteristic of cowpea's subtropical environment.

Phylogenetic analysis revealed the closeness of evolution among HSP60 sequences, suggesting a common ancestor and underscoring possible functional conservation within some clades^[48]. In previous studies, HSP60 has been classified according to the predicted subcellular localizations. For example, the evolutionary study of HSP60s in the model plant A. thaliana partitioned them into the mitochondria and chloroplasts^[42]. Similarly, three groups of VuHSP60 genes were identified; these genes are primarily predicted to localize in the cytoplasm, mitochondria, and chloroplasts (Fig. 2, Table 1), with a mix of monocot and dicot species, which suggests that these HSP60 genes diverged before the monocot–dicot split. Group IV completely lacks both cowpea and soybean representatives (Fig. 2). This selective absence suggests a lineage-specific gene loss event or a distinct evolutionary divergence within the Fabaceae family that warrants further genomic investigation. According to the domain architecture analysis, all identified proteins contained the highly conserved Cpn60_TCP1 domain, which is essential for ATP-dependent protein folding activity. In addition, significant structural diversity in gene architecture was observed, with intron counts ranging from 0 to 17 (Fig. 3). This architectural variability aligns with recent findings reported in other crops, including tea^[14] and cassava (Manihot esculenta)^[49]. This structural variation in VuHSP60 genes reflects a history of intron loss and gain events. Low-intron variants are structurally optimized for rapid transcription by minimizing RNA splicing requirements, whereas multi-intron genes provide the architectural framework that supports alternative splicing.

In the current study, the predicted multigene network showed a centralized post-transcriptional regulatory framework for cowpea chaperonins. Screening based on psRNATarget categorized these interactions into RNA cleavage and translation inhibition mechanisms across eight miRNA families (Fig. 4). The network topology demonstrated that growth- and stress-associated regulators, particularly Vun-miR396 and Vun-miR164, exhibit high sequence complementarity with multiple VuHSP60 paralogs simultaneously. This multitarget architecture suggests an evolutionary mechanism where single small RNA nodes synchronize the fine-tuning of multiple chaperonin loci^[50]. Ultimately, these in silico interactions provide a robust candidate pool for future experimental validation of cowpea's stress response networks.

Characterization of the upstream promoter regions revealed a diverse architectural layout of cis-regulatory elements across the VuHSP60 gene family (Fig. 6). The diverse promoter profiles predicted by the PlantCARE database suggest a highly integrated transcriptional network for the VuHSP60 family. The high abundance of predicted ABREs) and MeJA-responsive elements suggest that these chaperonin loci are computationally located within core phytohormone-mediated stress cascades^[51,52]. Rather than operating in isolation, the predicted coexistence of these hormone-responsive components alongside specialized motifs, such as LTRs, defense-related TC-rich repeats, and AREs, points to potential functional diversification at the promoter level. This specialized architectural layout reveals a structural framework that supports the potential for differential or selective transcript activation among individual VuHSP60 paralogs. Finally, the present in silico multi-stress framework provides a mechanistic hypothesis of cowpea's molecular chaperone system to coordinate its transcription under environmental fluctuations, establishing a defined target pool for future functional validation under low and high temperature and hormone treatments.

According to the qRT-PCR expression profiling, VuHSP60 gene family members exhibit differential transcriptional regulation under heat and cold stress (Fig. 7), consistent with findings reported in S. bicolor^[41] and C. sinensis^[14]. The significant upregulation of nine genes under heat stress and 10 under cold stress highlights distinct transcriptional variations within this gene family. Although the majority of these genes displayed strict stimulus-specific expression (seven exclusively under heat and eight exclusively under cold), VuHSP60-1 and VuHSP60-12 exhibited dual responsiveness to both treatments. This overlapping accumulation suggests that these two specific family members are transcriptionally activated by upstream signaling networks common to both heat and cold stress. Similar transcriptionally distinct profiles within chaperone families have been documented in other legumes such as V. radiata, where specific VrHSP60 genes displayed pronounced, stress-specific expression patterns during thermal stress. For example, VrHSP-60.22 showed differential transcription depending on the stress type, exhibiting responsiveness to both heat and cold stress, whereas other chaperones in mungbean demonstrated strict single-stimulus induction^[40]. Similarly, cowpea displays specific upregulation of VuHSP60 genes exclusively under high- and low-temperature conditions. Collectively, these expression variations identify specific VuHSP60 candidates that are responsive to temperature stress; however, further functional studies are required to determine whether these transcript-level changes translate to altered protein abundance or functional differences.

This study establishes the first genome-wide characterization of the HSP60 gene family in cowpea and thus provides a detailed baseline for chaperone architecture in legumes. The data also reveal structural conservation across paralogs by showing that family expansion was driven by segmentally duplicated events under purifying selection. The different transcriptional profiles obtained by qRT-PCR showed that VuHSP60-1 and VuHSP60-12 are candidate genes that are responsive to high and low temperatures. These in silico and expression datasets provide precise, high-priority targets to guide future functional validation and molecular breeding strategies for climate-resilient cultivars.

This work presents the complete genomic analysis of the HSP60/Cpn gene family in cowpea, characterizing a structurally diverse group of 25 members shaped by segmental duplication and purifying selection. The differential and overlapping expression profiles of these genes at severe temperatures suggest their specific potential functions in maintaining cellular integrity. Instead of a mere roster of stress-responsive candidates, this work identified VuHSP60-1 and VuHSP60-12 as candidate genes for future functional studies. In addition, protein interactions, and multi-omics approaches should be integrated to further investigate the specific regulatory networks and translate these findings into robust vegetable cultivars suited to increasing climatic variability.