Gene Co-expression Network Analysis Identifies Stress-Specific Modules and Hub Genes Underlying Abiotic-Stress Responses in Chili Pepper (Capsicum annuum L.)

1. Introduction

Chili pepper (Capsicum annuum L.) is one of the most widely cultivated and economically important vegetable and spice crops in the world, valued for its fruit quality, pungency, colour, and rich content of vitamins, carotenoids and capsaicinoids [1]. It is a dietary staple across much of Asia, Africa and Latin America, and India is among the largest global producers, consumers and exporters of dried chili [2]. Beyond its culinary role, the crop is a significant source of income for smallholder farmers, making the stability of its yield a matter of both food and economic security.

Chili productivity is, however, repeatedly constrained by abiotic stresses, principally low and high temperature, soil salinity and water deficit, which impair seed germination, vegetative growth, flowering, pollination and fruit set, and can cause substantial yield losses [3]. As cultivation expands into marginal environments and as climatic variability intensifies, these stresses are expected to become more frequent and severe [4]. Developing stress-tolerant chili cultivars is therefore a priority for breeding programmes, and this in turn depends on a clear understanding of the molecular mechanisms underlying the plant’s response to each stress.

Plant responses to abiotic stress are highly polygenic and coordinated, involving the simultaneous reprogramming of hundreds to thousands of genes that act in signal perception, transcriptional regulation, osmotic and ionic adjustment, reactive-oxygen detoxification, hormone signalling and metabolic adaptation [5]. High-throughput RNA-sequencing has made it possible to capture these genome-wide responses, and numerous transcriptomic studies have catalogued stress-responsive genes in pepper and related Solanaceae [6]. However, analyses based on differential expression alone identify which genes change in abundance but not how those genes are organised into coordinated regulatory programmes, and they tend to highlight the most strongly responding genes rather than the most influential ones.

Weighted gene co-expression network analysis (WGCNA) provides a complementary, systems-level framework that addresses this limitation [7]. By grouping genes with highly correlated expression into modules and identifying the most interconnected “hub” genes within each module, WGCNA exploits the principle that co-expressed genes frequently participate in shared biological processes (“guilt by association”) [8]. Modules can then be related to external traits or conditions, and their hub genes nominated as candidate regulators. The approach has been applied successfully across many plant species to dissect developmental and stress-response programmes and to prioritise candidate genes for functional study and breeding [9,10]. In Capsicum, co-expression approaches have begun to be used, including a global stress co-expression network constructed from large RNA-sequencing compendia [11]; however, focused analyses that resolve the module structure and hub genes associated with individual abiotic stresses in pepper remain limited.

In the present study, we constructed a weighted gene co-expression network from a comprehensive C. annuum abiotic-stress transcriptome spanning cold, heat, salinity and osmotic treatments together with controls, with the aim of dissecting the modular organisation of the pepper abiotic-stress response. Our specific objectives were (i) to define the co-expression modules of the pepper abiotic-stress transcriptome, (ii) to identify the modules most strongly and specifically associated with each of the four stresses, (iii) to nominate the hub genes of these stress-associated modules as candidate regulators and confirm their sequence identities, and (iv) to assess the robustness of the resulting modules and hubs. By focusing on stress-specific modules and their sequence-confirmed hub genes, this work aims to provide a prioritised set of candidate genes to support the dissection and improvement of abiotic-stress tolerance in chili pepper.

2. Results

2.1. Construction of the Weighted Gene Co-Expression Network

After filtering and variance stabilisation, the 15,000 most variable genes across all 78 samples were used for network construction. Hierarchical clustering of samples confirmed coherent grouping by treatment and revealed no outlier samples (Figure 1). The scale-free topology index increased steadily with soft-thresholding power but did not exceed the conventional 0.85 threshold at biologically reasonable powers; a power of 16 was therefore adopted following the recommended default for signed networks of this sample size, at which the model fit reached R² = 0.73 and mean connectivity declined to approximately 51 (Figure 2).

2.2. Identification of Co-Expression Modules

Dynamic tree cutting with subsequent merging of highly similar modules partitioned the 15,000 genes into 40 co-expression modules, ranging from 54 genes (orangered4) to 2,509 genes (turquoise) (Figure 3). The modules varied considerably in size, with 10 modules containing more than 500 genes and the remainder representing smaller, more specialised gene sets. Only 165 genes (~1.1%) remained unassigned to any module (grey), indicating that the network captured the dominant transcriptional structure of the dataset. The colour band beneath the gene dendrogram (Figure 3) shows that the modules corresponded to coherent branches of the clustering tree, supporting the validity of the module partition. Each module was summarised by its module eigengene for subsequent association with the abiotic-stress traits.

2.3. Module–Trait Relationships

Correlation of module eigengenes with stress traits revealed that each of the four stresses was most strongly associated with a distinct module, all highly significant (Table 1; Figure 4). Cold stress was exceptionally tightly associated with the orange module (r = −0.91, p = 9 × 10⁻³¹³¹). Heat stress was strongly associated with the turquoise module (r = +0.83, p = 7 ×⁻²¹10⁻²¹), the largest module in the network, indicating a broad, coordinately up-regulated heat-response programme. Salinity was most associated with the steelblue module (r = +0.46, p = ⁻⁵ × 10⁻⁵) and osmotic stress with the pink module (r = −0.41, p ⁻⁴ 2 × 10⁻⁴). The negative coefficients for cold and osmotic indicate coordinated down-regulation of module genes under those treatments. A white module correlated with the general stressed-versus-mock contrast (r = +0.66, ⁻¹¹ = 5 × 10⁻¹¹), and a greenyellow module declined over the treatment time course (r = −0.6⁻¹¹, p = 4 × 10⁻¹¹).

2.4. Hub Genes Within Stress-Associated Modules

Hub genes of each stress-associated module were identified by high module membership (kME) combined with high gene significance (GS), and the identities of the top cold- and heat-module hubs were confirmed by BLASTP against the NCBI non-redundant database. The most connected hub of the cold-associated orange module (CA08g16690; kME = 0.90) was a telomere repeat-binding factor, a Myb-like telobox DNA-binding protein of a family increasingly implicated in transcriptional regulation beyond telomere maintenance. A further high-ranking hub (CA01g02040; kME = 0.88) was 1-aminocyclopropane-1-carboxylate (ACC) synthase, the rate-limiting enzyme of ethylene biosynthesis, consistent with the role of ethylene in cold-stress signaling. A notable proportion of additional cold-module hubs corresponded to conserved but functionally uncharacterised proteins, representing novel candidate regulators (Table 2).

Notably, BLASTP comparison corrected several reference-genome annotations: the top hub (CA08g16690) was mis-annotated in the CM334 reference as a myosin-related protein but is unambiguously a telomere repeat-binding factor (99% identity), and two further hubs were similarly reassigned (Table 2). This underscores the value of sequence-level identity checks when nominating candidate genes from a non-model crop with sparse functional annotation.

The top hub genes of the heat-associated turquoise module were likewise checked by BLASTP and are listed in Table 3. The most connected heat-module hub (CA06g08540; kME = 0.93) was an E3 SUMO-protein ligase SIZ1 (96.5% identity), a known positive regulator of thermotolerance, alongside photosynthesis- and metabolism-related hubs and several further conserved uncharacterised proteins (Table 3).

2.5. Functional Enrichment of Stress-Associated Modules

Because a species-specific Gene Ontology resource was not available for Capsicum annuum through the standard annotation servers, the functional composition of the stress-associated modules was examined by assigning each gene to a curated functional category on the basis of its reference protein annotation, and testing each category for over-representation within a module relative to the whole co-expression network (hypergeometric test, Benjamini–Hochberg correction).

Two modules showed statistically significant functional over-representation. The heat-associated turquoise module was enriched for genes of the “stress/defence” category (89 genes; 3.5% of the module versus 2.3% of the network; fold = 1.52; FDR = 3.0 × 10⁻⁴⁴), consistent with its strong correlation with heat treatment and with the heat-responsive hub genes identified above. The osmotic-associated pink module was enriched for “RNA / translation”-related genes (26 genes; 4.8% versus 2.4%; fold = 1.99; FDR = 0.011), in agreement with the established sensitivity of translation and ribosome biogenesis to osmotic and water-deficit stress (Table 4).

Several further associations were suggestive but did not survive multiple-testing correction and are reported only as trends: the salinity-associated steelblue module showed nominal enrichment of lipid-metabolism, cell-wall and transcription-factor categories (fold = 2.2–6.8; FDR ≈ 0.22), and the cold-associated orange module showed nominal excess of hormone-related and transcription-factor genes (fold = 2.2 and 1.6; FDR ≈ 0.84) (Table 4).

Notably, the cold-associated orange module despite its exceptionally strong correlation with cold treatment (r = −0.91) showed no functional category that was significantly over-represented after correction. Rather than indicating an absence of biological coherence, this reflects that the module is composed of functionally diverse genes spanning many categories, together with a substantial proportion (collectively more than two-thirds across the stress modules) of genes annotated only as uncharacterised, hypothetical, or of unassigned function. This high fraction of poorly annotated genes is characteristic of the Capsicum genome and underscores that a considerable part of the stress-responsive co-expression programme is mediated by genes of currently unknown function themselves a valuable set of candidates for future characterisation.

Taken together, the functional-composition analysis provides supporting evidence that the heat and osmotic modules are enriched for biologically appropriate processes, and indicates that the stress-associated modules are coordinated by functionally heterogeneous gene sets rather than by single dominant pathways. The most informative functional insights at the gene level came from the individual hub genes with confirmed identities, with the module-level composition providing their broader context.

2.6. Robustness of the Modules and Hub Genes

To confirm that the co-expression modules represented robust transcriptional structure rather than features of a particular sample partition, module preservation was assessed by comparing the full network with a random half of the samples (modulePreservation, 30 permutations). All assigned modules were preserved, with composite Zsummary values exceeding 10 the conventional threshold for strong preservation throughout, the only exception being the unassigned grey module (Zsummary = 8.0), as expected. The four stress-associated modules were all preserved: the osmotic-associated pink (Zsummary = 49.6) and heat-associated turquoise (Zsummary = 44.1) modules most strongly, the cold-associated orange module robustly (Zsummary = 24.6), and the salinity-associated steelblue module the smallest stress module above the strong-preservation threshold (Zsummary = 13.3). This indicates that the modular organisation recovered by the network was reproducible and stable.

The stability of the leading hub genes was further examined by bootstrap resampling. The most central hub of each stress-associated module was consistently recovered, remaining among the most highly connected genes of its module across resamples, whereas the exact rank ordering among lower-ranked hubs of closely similar connectivity was less stable, as expected. Together, these analyses indicate that both the stress-associated modules and their leading hub genes are robust features of the pepper abiotic-stress transcriptome.

3. Discussion

This study used weighted gene co-expression network analysis to resolve the modular organisation of the Capsicum annuum abiotic-stress transcriptome and to nominate candidate regulators of stress responses. The network partitioned cleanly into 40 modules with very few unassigned genes (~1.1%), and, most strikingly, each of the four abiotic stresses mapped onto a distinct, statistically robust module, demonstrating that cold, heat, salinity and osmotic stress drive separable co-expression programmes in pepper. Module-preservation analysis confirmed that all modules, including the four stress-associated modules, were strongly preserved under resampling (Zsummary > 10), indicating that this modular organisation reflects reproducible transcriptional structure rather than features of a particular sample partition.

The cold response stood out for the exceptional strength of its module association (r = −0.91), among the stronger module–trait correlations reported in plant WGCNA studies, and the negative sign indicates that the orange-module programme is predominantly repressed under cold. Its hub genes offer a coherent biological narrative. The most connected hub was a telomere repeat-binding factor, which remained the leading hub of the module across bootstrap resamples; while historically associated with telomere maintenance, this Myb-like telobox class of DNA-binding proteins has increasingly been linked to transcriptional regulation and stress responses [14], making its central position in a cold module a noteworthy and testable hypothesis. A second hub, ACC synthase, is the rate-limiting enzyme of ethylene biosynthesis, and the connection between ethylene and cold/chilling responses is documented across plant systems [15], lending biological plausibility to the module. Notably, a substantial fraction of cold-module hubs was conserved but functionally uncharacterised proteins; rather than a limitation, these represent the kind of novel candidate regulators that network analysis is uniquely able to surface.

The heat response was governed by the largest module in the network, consistent with the broad, system-wide transcriptional reprogramming that high temperature elicits, and its top hub genes were checked individually by BLASTP. The most compelling hub was an E3 SUMO-protein ligase SIZ1 (CA06g08540; 96.5% identity), a master regulator of stress-responsive protein SUMOylation with a well-documented role in heat-stress signalling and acquired thermotolerance [16], lending strong biological coherence to the module. Other hubs spanned processes central to the heat response, including photosynthetic machinery (a Photosystem II CP43 reaction-centre protein), membrane lipid metabolism (phosphatidate cytidylyltransferase), aromatic-amino-acid and phenylpropanoid precursor supply (arogenate dehydratase), auxin-responsive growth regulation (a SAUR-family protein) and specialised oxidative metabolism (an SRG1-like 2-oxoglutarate-dependent dioxygenase). Consistent with this hub composition, the heat module was significantly over-represented for genes of the stress/defence functional category (FDR = 3.0 × 10⁻⁴⁴), providing module-level support for its heat-responsive character independently of individual hub annotations. Intriguingly, a BAHD acyltransferase related to the capsaicinoid-biosynthesis enzyme Pun1 [17] was also among the hubs, hinting at a possible connection between pungency-related specialised metabolism and the heat response that merits further investigation. As in the cold module, several high-connectivity hubs corresponded to conserved but functionally uncharacterised proteins, representing novel candidate regulators.

The salinity- and osmotic-associated modules, although significant, showed more moderate correlations, which may reflect the partially overlapping signalling of ionic and osmotic stress and the more diffuse transcriptional signatures these treatments produce. Nonetheless, the osmotic-associated module was significantly enriched for genes of the RNA/translation category (FDR = 0.011), in line with the established sensitivity of translation and ribosome biogenesis to osmotic and water-deficit stress, while the salinity-associated module showed suggestive but non-significant trends towards lipid-metabolism, cell-wall and transcription-factor categories. By contrast, and despite its very strong eigengene correlation, the cold-associated module showed no single significantly over-represented functional category, indicating that it is composed of functionally heterogeneous genes including many of unknown function rather than being dominated by a single pathway.

Several limitations should be acknowledged. The analysis is based on a single, albeit comprehensive, publicly available dataset, and co-expression relationships are correlative and do not establish causation. Functional annotation of pepper remains incomplete: a substantial proportion of the genes in the stress-associated modules, including several high-connectivity hubs, correspond to uncharacterised or hypothetical proteins, and the absence of a species-specific Gene Ontology resource limited the functional analysis to curated functional categories rather than formal GO/KEGG enrichment. In addition, the mapping between annotation versions, although supported by sequence homology, introduces a degree of uncertainty for individual genes. These constraints notwithstanding, the convergence of strong and robustly preserved module–trait associations, bootstrap-stable hub genes, sequence-confirmed hub identities and biologically coherent functional enrichment provides multiple, mutually supporting lines of evidence for the stress-associated modules and their candidate regulators.

Future work should prioritise experimental validation of the nominated hub genes, for example by quantitative RT-PCR across the stress time courses and by functional perturbation (overexpression or gene editing) of the highest-priority candidates, such as the cold-module telomere repeat-binding factor and the heat-module SUMO ligase SIZ1. Obtaining comprehensive Gene Ontology and pathway annotation for Capsicum annuum, and integrating additional pepper stress datasets, would further test the generality of the modules identified here and convert the candidate regulators reported in this study into experimentally validated targets for the breeding of abiotic-stress-tolerant chili pepper.

4. Materials and Methods

4.1. Transcriptome Data Source

Publicly available RNA-sequencing data for Capsicum annuum cv. Criollo de Morelos 334 (CM334) were retrieved from the NCBI Gene Expression Omnibus (accession GSE132824; [12]). The dataset comprises 78 leaf samples representing four abiotic stresses cold (10 °C), heat (40 °C), salinity (400 mM NaCl) and osmotic stress (400 mM mannitol) together with mock-treated controls. Stress samples were collected at 3, 6, 12, 24 and 72 h after treatment and controls at 0, 3, 6, 12, 24 and 72 h, each with three biological replicates. Gene-level read counts (35,884 gene models) were obtained from the processed expression matrix provided with the series.

4.2. Data Preprocessing and Normalisation

Analyses were performed in R (v4.0.5). Genes expressed at low levels were removed by retaining only gene models with at least 10 read counts in at least 25% of samples. The filtered count matrix was variance-stabilised using the vst function of DESeq2 (v1.30.1) [13]. To focus the network on informative variation and reduce computational load, the 15,000 genes with the highest variance across samples were retained for co-expression analysis. Sample quality was assessed with the goodSamplesGenes function of the WGCNA package (v1.70.3), and potential outliers were examined by average-linkage hierarchical clustering of samples; no outlier samples were detected, and all 78 samples were retained.

4.3. Co-Expression Network Construction and Module Detection

A signed weighted gene co-expression network was constructed using the WGCNA package (v1.70.3) in R. The soft-thresholding power was evaluated over the range 1–30 using pickSoftThreshold with a signed network type. The scale-free topology fit index did not reach the conventional 0.85 threshold until a power of 24, at which network connectivity was severely reduced (mean k ≈ 14), a behaviour commonly observed for datasets dominated by several strong treatment effects. A soft-thresholding power of 16 was therefore adopted as a compromise that retained adequate connectivity (mean k ≈ 51) while approaching scale-free topology (R² = 0.73), consistent with the WGCNA authors’ recommended default for signed networks of comparable sample size. The network and modules were generated in a single block with the blockwiseModules function using the following parameters: soft power 16, signed network and signed topological overlap matrix (TOM), minimum module size 30, merge cut height 0.25, reassignment threshold 0, and maximum block size 20,000. This procedure yielded 40 modules; 165 genes (~1.1%) that were not assigned to any module were placed in the grey module. Modules were summarised by their module eigengenes (the first principal component of each module’s expression), and similar modules were merged at a dissimilarity threshold of 0.25.

4.4. Module–Trait Association Analysis

Sample traits were encoded numerically for correlation analysis: each of the four stresses was represented as a binary indicator (treated = 1, otherwise = 0), treatment duration was encoded as a continuous variable (hours), and a binary stressed-versus-mock indicator was included; mock-control samples therefore served as the common reference baseline. Pearson correlation coefficients between module eigengenes and traits were calculated, and their significance was assessed using Student asymptotic p-values (corPvalueStudent). Modules showing the strongest and most significant correlation with each stress trait were selected for downstream analysis.

4.5. Identification of Hub Genes

For each stress-associated module, intramodular connectivity was quantified as module membership (kME), defined as the correlation between each gene’s expression and the module eigengene (signedKME). Gene significance (GS) was defined as the absolute correlation between each gene’s expression and the trait of interest. Hub genes were defined as those exhibiting both high module membership (|kME| ≥ 0.8) and high gene significance for the corresponding stress.

4.6. Gene Identifier Mapping and Hub-Gene Identity Confirmation

Because the expression dataset used Pepper Genome Annotation scaffold-based identifiers (e.g., TC.CA.PGAv.1.6.scaffold606.26), whereas the reference proteome used chromosome-level identifiers (e.g., CA08g16690), gene models were mapped between the two systems using the CM334 annotation assignment table (Pepper v.1.55, Sol Genomics Network). For each scaffold, genes were ordered by chromosomal coordinate, and the scaffold-internal index of each expression-dataset identifier was matched to the correspondingly ordered chromosome-level identifier. The reliability of this mapping was confirmed by independent BLASTP searches of the mapped protein sequences against the NCBI non-redundant (nr) protein database, which returned high-identity matches to Capsicum proteins. The identities of the top hub genes were assigned from their BLASTP best hits.

4.7. Functional Composition Analysis of Modules

Because a species-specific Gene Ontology resource was not available for Capsicum annuum through standard annotation servers, the functional composition of the stress-associated modules was characterised using the functional descriptions of the CM334 reference proteome (Pepper v.1.55, Sol Genomics Network). Each gene was assigned to one of a set of curated functional categories (e.g., transcription factor/DNA-binding, protein kinase, heat shock/chaperone, stress/defence, hormone-related, photosynthesis, redox, transport, lipid metabolism, carbohydrate metabolism, RNA/translation) on the basis of keyword matching to its protein annotation. For each stress-associated module, over-representation of each functional category was tested against the whole co-expression network as background using the hypergeometric test, and p-values were corrected for multiple testing by the Benjamini–Hochberg method; categories with an adjusted p-value below 0.05 were considered significantly over-represented.

4.8. Robustness Assessment

The robustness of the co-expression modules was evaluated by module-preservation analysis using the module Preservation function of WGCNA, with the full dataset as the reference network and a random 50% subsample of the samples as the test network (30 permutations, signed network). Preservation was summarised by the composite Zsummary statistic, for which values above 10 indicate strong preservation and values between 2 and 10 indicate moderate preservation. The stability of hub genes was assessed by bootstrap resampling: samples were resampled with replacement 100 times, module membership (kME) was recomputed in each resample, and the proportion of resamples in which each of the ten leading hub genes of a stress-associated module remained among the ten most highly connected genes of that module was recorded.

5. Conclusions

This study provides a systems-level view of the Capsicum annuum abiotic-stress transcriptome, resolving it into 40 co-expression modules in which each of the four stresses cold, heat, salinity and osmotic is captured by a distinct, robustly preserved module. By combining module–trait association, module preservation and hub bootstrapping with sequence-level (BLASTP) confirmation of hub identities, the analysis nominates a prioritised and reliability-graded set of candidate regulators rather than an unfiltered gene list. The cold module, the most strongly associated of the four (r=−0.91), is centered on a telomere repeat-binding factor and an ethylene-biosynthesis enzyme (ACC synthase), while the heat module is anchored by the SUMO E3 ligase SIZ1, a known regulator of thermotolerance; alongside these, a substantial number of conserved but uncharacterised hub proteins emerge as novel candidates that would be difficult to surface by single-gene approaches. Because the findings derive from a single public dataset and from correlative relationships, they are best regarded as well-supported, testable hypotheses rather than established mechanisms. Their value lies in focusing future experimental effort: the hub genes identified here, in particular the cold-module telomere repeat-binding factor and the heat-module SIZ1, represent immediate priorities for functional validation and, ultimately, for the molecular breeding of abiotic-stress-tolerant chili pepper.