In the previous study, total RNA was extracted from cold-treated leaves and roots of pepper (Liu et al., 2017). The elite breeding pepper (Capsicum annuum) line 6421 was selected from a long red pepper landrace grown widely in the West of the Xiangjiang River, in Hunan Province, China. It is resistant to anthracnose, bacterial spot, and bacterial wilt. Control plants were mock-treated with nutrient solution alone. Leaf and root tissues were collected from both treated and control plants at 0.5, 1, 3, 6, 12, and 24 h post-treatment (HPT).

We performed quality analysis using fastp (Chen et al., 2018) and FastQC v0.11.7 (Andrews, 2017), and aligned sequences with the reference genome of pepper (Zunla genome) (Qin et al., 2014). Default mapping parameters (10 mismatches/read; nine multi mapping locations/read) were analyzed using HISAT 2.2.1 (Kim D. et al., 2015). We used v1.20.0 of DESeq2, an R-based software package provided by Bioconductor, for differential gene expression analysis (Love et al., 2014).

The abundance of transcripts was measured as the average normalized count of the reads mapped to the transcript, and the difference in their abundance was examined under two conditions to identify transcripts that are differentially expressed under both conditions (Love et al., 2014). The difference in expression was quantified based on the logarithm (change in the logarithmic multiple) of the average normalized count ratio between the two conditions. The differentially expressed transcripts in our experiment were defined as those with adjusted P-values < 0.01; cut-off threshold, |log2fold change (FC)| ≥ 2 (negative binomial Wald test, followed by the Benjamini-Hochberg correction). Differentially expressed genes (DEGs) were classified as up-regulated or down-regulated genes based on their significant positive or negative logarithmic changes in value.

Venn diagrams were constructed (Lin et al., 2016). Heatmaps were generated using the seaborn heatmap available in python. Statistical sequencing reads have been included in Supplementary Table 1 and all DEGs genes information in different stages in Supplementary Table 2.

Co-expression analysis was performed on samples from 12 control tissues and 12 cold stress-treated tissues using the k-means (Gasch and Eisen, 2002) method in python. The normalized expression value of genes was calculated by dividing their expression levels in all samples by their maximum observed transcripts per million (TPM); the cluster information for each expressed gene is displayed in Supplementary Table 3. Principal component analysis (PCA) and hierarchical clustering (HCL) were performed using the Kernel principal component analysis method in python; it was convenient to graphically explain the correlations between all samples using default settings. Transformed normalized gene expression values with z-scores were used for PCA and hierarchical clustering.

The identified clusters that responded to cold stress were followed up via gene enrichment analysis with GOATOOLS (Klopfenstein et al., 2018), a python package used for gene ontology (GO) enrichment analysis and determination of false discovery rate (FDR) values of statistically significant GO terms (FDR < 0.01). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed by KofamKOALA (Aramaki et al., 2020) to protein sequences by homology search, and the enrichment analysis of the KEGG pathway was carried out R package clusterProfiler (Yu et al., 2012; P value Cutoff = 0.05).

The Arabidopsis genome website contains Arabidopsis gene sequence information¹ and enables us to determine the Arabidopsis CBF pathway through a comparison of the protein sequence generated by the Arabidopsis gene and the protein sequence generated by the pepper genome gene. Related genes corresponded to orthologous genes in pepper. For important genes in the ICE-CBF-COR-signaling module, the pepper genes with the closest corresponding homology relationship that contain specific structural domains were also identified through homology comparison.

Fourteen key genes were identified in the ICE-CBF-COR signaling module in pepper, and their expression under control and cold stress-treated conditions was studied. The TPM value of gene expression was normalized using the z-score (Supplementary Table 4), and the expression was plotted using a heat map. Then, we used cold tolerance (CT) inbred line A188 and a cold sensitive (CS) inbred line A122 under cold-rewarm treatments and used the RNA-seq data (accession number: PRJNA646356), as described by Miao et al. (2021), to determine the expression levels of the 14 identified genes; then, the Z-score was normalized and a heat map was drawn.

We performed Pearson correlation analysis of a total of 29,053 genes that were expressed in at least one sample, obtained the correlation matrix between gene pairs, and extracted 14 CBF pathway-related genes and their correlation information. Python NetworkX enabled us to determine the network relationship diagram between 14 genes.

In order to thoroughly observe the response network of genes in the leaves and roots of pepper plants before and after their exposure to cold stress, we analyzed the complete gene expression in the roots and leaves of control and stress-treated seedlings at six time points, and established a root-based and spatiotemporal dynamic expression data set of whole genes in the leaf tissue. PCA results for the transcriptome showed that the sample exhibited significant tissue specificity (Figure 1A). PCA also showed that different samples can be separated by both PC1 (41.76%) and PC2 (10.58%). In the early stage, when leaves were subjected to cold stress, at 0.5–1 h, the untreated leaves were closer together, and were clustered together at 3–24 h, indicating that the CORs in the leaves of peppers were induced to express proteins 3 h after cold stress treatment. Between 0.5 and 24 h of the root tissue being subjected to cold stress, three main clustering modules could be observed according to the time ranges of 0.5–1, 3–6, and 12–24 h. We circled them with the same colors. The difference in the expression at the root tissues of the three of these and the control showed an increasing trend, and a clear distinction was observed between different treatment times (Figure 1A). It can also be seen that after pepper was subjected to cold stress, the genes in its leaves and roots had strong tissue specificity with regard to the speed and degree of response to cold stress. The Pearson correlation analysis of datasets displayed a gene expression pattern and significant tissue specificity. The same tissues are clustered together (Figure 1B), which is consistent with the phenomenon shown in Figure 1A.

In order to further analyze whether the dynamic expression of genes in the roots and leaves of peppers after cold stress treatment exhibits a certain trend of differentiation, we used the k-means clustering algorithm to divide them into at least one sample according to the expression pattern of 29,053 genes. After 15 clusters, the clustering information for each gene is shown in Supplementary Table 3. Upon analyzing these 15 clusters (Figure 2), we found that the expression levels of 1,373 genes in cluster 2 were much higher in the leaf tissues than in the control, 3–24 h after cold stress treatment. The genes in the cluster responded to cold stress. However, the expression levels of genes in cluster 3 were much higher in the leaf tissues under stress at 24 h than those observed in the untreated samples at the same time point when a peak was observed. The expression levels were almost the same at other time points. The 15 clusters formed according to the trend of dynamic expression for the entire gene set had different expression patterns, and gene expression patterns in the same cluster were almost the same. We can classify the COR gene response by time and organization, according to the clustering information. However, we simultaneously found that after cold stress treatment, the expression of some genes in clusters did not change significantly at each sampling time point in the roots and leaves, such as in cluster 5, cluster 10, and cluster 12. Genes may not be sensitive to cold stress, and their expression is always stable.

In order to analyze the functional situation of gene clusters with increased expression after cold stress, the three clusters of the cluster 3, cluster 9, and cluster 15 sets were subjected to a GO term enrichment analysis (Figures 3A–C). Those genes in the cluster 3, 9, and 15 were significantly enriched in the “protein-chromophore linkage,” “photosynthesis,” “peptide biosynthetic process,” “amide biosynthetic process,” and “phosphatidylinositol phosphate biosynthetic process” biological processes. The top enriched GO terms in “molecular functions” were related to “protein serine/threonine kinase activity,” “phosphotransferase activity,” “chlorophyll binding,” “structural constituent of ribosome,” “transferase activity,” and “DNA-binding transcription factor activity.” In addition, the top enriched GO terms in “cellular component” were “chloroplast thylakoid membrane,” “plastid thylakoid membrane,” “intracellular non-membrane-bounded organelle,” and “ribosome.”

In order to gain insight into the response mechanism of peppers to cold stress, we conducted a comparative analysis of the transcriptional differences between the leaves and roots treated at 10°C at six time points (Figure 4). DEGs among samples were defined using fold change values, assessed using the expression of assembled transcripts. We have identified 2,306 up-regulated genes and 1,543 down-regulated genes in the leaves of peppers at different processing times. We have also identified 1,551 up-regulated genes and 933 down-regulated genes in the roots; among these, 169 genes exhibited two states of up-regulation and down-regulation in six stress-treated and control leaves. There were 79 genes in the roots in two states, and the remaining DEGs were up-regulated or down-regulated at all periods. Removed genes were identified in both leaves and roots, and the analyzes collectively yielded 4,872 DEGs, constituting ∼16.8% of the expressed genes in the dataset.

In the dynamic situation of gene expression within 0.5–24 h of pepper plants being subjected to cold stress, both tissues had the largest number of DEGs at 24 h. This indicates that a longer duration of exposure to cold stress leads to more profound changes in the pepper transcriptome profile. Currently, there were 1,123 up-regulated genes and 597 down-regulated genes in the roots, and 1,624 up-regulated genes and 1,107 down-regulated genes in the leaves. As shown in Figure 4, most of the differentially expressed genes were unique at 24 h. There were 784 unique up-regulated genes and 490 unique down-regulated genes in the roots, and 1,178 unique up-regulated genes and 964 unique up-regulated genes in the leaves. Among the unique down-regulated genes, the longer the pepper plants were subjected to treatment at 10°C, the greater the change in gene expression at the genome-wide level, and the number of differentially expressed genes in the leaves at 24 h was greater than that in the roots, indicating that pepper leaves have a stronger response to cold stress.

In order to further analyze the possible functions of the identified cold-responsive genes shown in Figure 4, the up-regulated and down-regulated differentially expressed genes (DEGs) identified in the leaf (Figures 5A,B) and root (Figures 5C,D) were used to perform GO enrichment analysis, and Venn diagrams were used to show the aggregation of up-regulated and down-regulated DEGs in the leaf and root, respectively (Figure 5C). In the leaf and root (Figures 5A,B,D,E), among the “biological process” category, the most significantly enriched terms were “olefinic compound metabolic process,” “phototropism,” “response to water,” and “cell wall organization.” The top enriched GO terms in “molecular functions” were related to “DNA-binding transcription factor activity,” “monooxygenase activity,” “oxidoreductase activity,” and “protein serine/threonine kinase activity.” In addition, the top enriched GO terms in “cellular component” were “integral component of membrane,” “intrinsic component of membrane,” “extracellular region,” and “nucleus.” Low temperatures can slow down the fluidity of cell membranes and stimulate plants to respond to the low temperature.

KEGG pathway enrichment was carried out in order to better understand the biological functions of cold stress DEGs in leaf and root. The enrichment analysis indicated that “MAPK signaling pathway,” “Phenylpropanoid biosynthesis,” “Cytochrome P450,” “Glutathione metabolism,” “Zeatin biosynthesis,” and “Diterpenoid biosynthesis” were common in leaf and root (Figure 6). However, “Plant-pathogen interaction”, “Protein kinases,” “Photosynthesis proteins,” “alpha-Linolenic acid metabolism,” “Cutin, suberine and wax biosynthesis” and “Brassinosteroid biosynthesis” were unique in leaf (Figure 6A) and “Pentose and glucuronate interconversions” and “Sesquiterpenoid and triterpenoid biosynthesis” were unique in root (Figure 6B).

Several studies have focused on ICE-CBF-COR in Arabidopsis and some model plants (Thomashow, 1999; Chinnusamy et al., 2003; Kim Y. S. et al., 2015; Barrero-Gil and Salinas, 2017). The ICE-CBF-COR-related genes in peppers have not been systematically identified and functionally analyzed. The expression of ICE-CBF-COR-related genes in peppers under low-temperature stress and their regulatory relationships needed to be examined. Therefore, based on the rich research background associated with ICE-CBF-COR in Arabidopsis, and based on the conservation of ICE-CBF-COR in the plant and protein sequence alignment, we identified that pepper contains 14 CBF1/2/3-related genes, including ICE1/2 and ICE1/2, and MYB15, which have been thoroughly studied in Arabidopsis (Agarwal et al., 2006), SOC1, positive regulators of EIN3, Brassinazole-resistant1 (BZR1) (Li et al., 2017b), CCA1, CESTA (CES) (Eremina et al., 2017), and positive regulators of LHY (Seo et al., 2009; Dong et al., 2011; Lee and Thomashow, 2012; Eremina et al., 2016). Notably, the normalized data for the control tissues obtained from the leaves and roots of pepper included that for the 14 orthologous genes in pepper and six time points at which transcript levels observed after cold stress enabled us to distinguish between negative regulators and positive regulators (Figure 7A). In Arabidopsis, the expression of CBF1 and CBF3 in response to a low temperature precedes that of CBF2, and CBF2 negatively regulates CBF1 and CBF3 expression (Novillo et al., 2004). In our pepper dataset, we can also see that CBF1 and CBF3 were clustered together, but the distance from CBF2 was high, indicating that the gene functions of CBF1/2/3 were conserved in pepper, implying that there is a certain degree of conservation in A. thaliana. Genes that are negatively or positively regulated by CBF also have similar effects in pepper, which is of great value for the study of the ICE-CBF-COR signaling module in pepper.

We used this data set to establish a co-expression regulatory network of these 14 genes in pepper and visualized the mutual regulatory relationship between these 14 genes. The gray lines represent negative regulatory relationships between genes, and the red lines represent the positive regulatory relationship between genes. Positive regulatory relationships have also been shown. The increased thickness of the line represents a stronger predicted relationship, and the dotted line represents the Pearson coefficient absolute value being less than 0.5. We found that the relationship predicted in pepper is the same as that observed in A. thaliana (Figure 7B). For example, MYB15, SOC1, and EIN3 negatively regulate ICE1 expression in Arabidopsis (Agarwal et al., 2006). The same result is also shown in our network diagram. The co-expression relationship between these 14 genes can provide certain insights into the correlation between related genes in pepper and the ICE-CBF-COR signaling module at low temperatures. Furthermore, we used predecessors in the cold-resistant (A122) and cold-sensitive (A188) varieties of pepper under cold stress (4°C for 0, 1, 2, and 24 h) conditions and after rewarming (28°C for 1 h), and obtained the RNA-seq data (Miao et al., 2021). We analyzed the expression of the 14 genes identified in the dataset and found that after the rewarming process, the expression of EIN3.1 and EIN3.2 was significantly up-regulated, and the expression of ICE1/2 in A188 was significantly higher than that in A122 (Figure 8).