Comparative Transcriptomics for Pepper ( Capsicum annuum L.) under Cold Stress and after Rewarming

: Cold stress has become one of the main abiotic stresses in pepper, which severely limits the growth and development of pepper. In this study, the physiological indicators and transcriptome of a cold-tolerance (CT) inbred line A188 and a cold-sensitive (CS) inbred line A122 under cold–rewarm treatments were studied; the aim of this study was to determine the potential of the key factors in pepper response to cold stress. Compared with CT, CS wilts more seriously after cold stress, with poor resilience, higher content of malondialdehyde, and lower content of soluble sugar and total chlorophyll. Moreover, during cold treatment, 7333 and 5953 differentially expressed genes (DEGs) were observed for CT and CS, respectively. These DEGs were signiﬁcantly enriched in pathways related to photosynthesis, plant hormone signal transduction, and DNA damage repair. Interestingly, in addition to the widely studied transcription factors related to cold, it was also found that 13 NAC transcription factors increased signiﬁcantly in the T4 group; meanwhile, the NAC8 (Capana02g003557) and NAC72 (Capana07g002219) in CT were signiﬁcantly higher than those in CS under rewarming for 1 h after 72 h cold treatment. Notably, weighted gene coexpression network analysis identiﬁed four positively correlated modules and eight hub genes, including zinc ﬁnger proteins, heat shock 70 kda protein, and cytochrome P450 family, which are related to cold tolerance. All of these pathways and genes may be responsible for the response to cold and even the cold tolerance in pepper.


Introduction
Temperature is a key factor governing the growth, geographic distribution, productivity, and seasonal behavior of plants [1,2]. With the influences of global climate change, the cold stress that occurs frequently and severely has been one of the major abiotic stresses affecting growth and development [3,4]. The interruption of photosynthetic electron transport, the carbon reduction cycle, and the decrease in stomatal conductance were induced by cold stress, which inhibited photosynthesis and reduced the accumulation of organic matter [5,6]. Abiotic stresses induce reactive oxygen species (ROS) accumulation, which results in DNA lesions-including DNA double-strand breaks, single-strand breaks, base modifications, and crosslinking-and even cause a massive loss of genetic information [7].
The stress response of plants is regulated by complex mechanisms that are helpful for reducing plant damage under cold stress [8]. The responses of plants to cold stress involve changes in photosynthetic pigmentation, protein composition, and abundance of plant cell membranes, especially thylakoid membranes, and also the significant downregulation of photosynthesis-related genes including PSI-, PSII-and LHC-polypeptides [9][10][11]. Plant hormones are also involved in response to cold stress [12]. Activating cold responsive (COR) genes, encoding cryoprotective proteins, protects plant cells from cold-induced damage [13]. Hormone components can activate COR genes through C-repeat binding factor (CBF)-dependent or nondependent pathways to regulate chilling tolerance in plants [14]. In addition, DNA damage response (DDR) is the basis of the body's defense against various kinds of damage. It has been reported that after cold-triggered DNA damage in Arabidopsis root stem cells, WEE1-mediated DDR can activate relevant signaling pathways to reduce cold damage to plants [15].
Transcriptome sequencing has developed as an important tool in genetic research to explain the molecular mechanisms of gene expression and improve the efficiency of identifying the genes of interest [16]. At present, comparative transcriptomics has been widely used to illustrate the complex genes regulatory networks of many plants under cold stress. Comparative transcriptome analysis in two peanut varieties with different cold-tolerant abilities showed that 445 TF genes were significantly differentially expressed under cold and revealed some possible mechanism in peanut cold tolerance [17]. The comparison transcriptome of the tea plant under cold found that 978 differentially expressed genes (DEGs) owned by the tolerant and sensitive plants were enriched in pathways of photosynthesis and hormone signal transduction, among others, and considered that the different expression of Lhca2 and SNF1-related protein kinase 2 (SnRK2) are correlated with cold tolerance [18]. Comparative transcriptomics of Poa pratensis systematically reveals the changes in key genes and products of glycolysis and TCA cycle in response to cold stress [19].
Pepper (Capsicum annuum L.) originating from tropical and subtropical areas is a warm season vegetable, and its optimal growth temperature ranges from 21 • C to 28 • C. Pepper plants lack mechanisms for cold acclimatization, making them very sensitive to low temperature [20]. Exposed to cold stress, pepper exhibit the reduction of nutriture and metabolic imbalances. Besides, cold stress impairs the function of the female organs and reduces the number of viable pollen grains per flower, which affects the growth, development, and productivity of pepper [21][22][23][24]. Severe stress may even lead to the death of pepper plants, resulting in serious economic losses [20]. Therefore, it is important to characterize mechanisms of tolerance to cold. The cold tolerance of plants is not only closely related to the resistance to cold stress, but also has an important relationship with the ability of plants to recover from damage after offsetting stress. Therefore, recovery after exposure to cold stress is also a relevant and exciting area of research [25]. Comparative transcriptome has also been applied to the study of stress in pepper, such as heat stress [26,27], multiple stress, hormone treatment [28], etc. However, there are few transcriptome studies of pepper related to cold. In particular, transcriptome studies focusing not only on the process of cold stress but also on the process of recovery process of pepper are scarcer. In the preliminary study, a cold-tolerance (CT) inbred line A188 and a cold-sensitive inbred line A122 (CS) were obtained through character investigation and screening. In the present study, pepper seedlings were exposed to an artificial climate chamber at 4 • C/28 • C (stress temperature/recovery temperature) to simulate cold stress and short-term recovery. Subsequently, the RNA-Seq analysis of CT and CS obtained gene expression profile and identified complex gene expression networks, which provide important information regarding the molecular basis of cold resistance adaptability in pepper, which may be useful for further improving the cold resistance of pepper.

Plant Materials and Treatments
Two pepper cultivars, a cold-tolerance inbred line A188 and a cold-sensitive inbred line A122, were provided by the Vegetable Institution of Hunan Academy of Agricultural Science. Seeds were sown in plastic pots with the nutrient substrate. In the greenhouse, plants were received unified routine management (~28 • C with 16 h light and~20 • C with 8 h dark) until the plants grew to the five-leaf stage. Seedlings were placed in an artificial climate chamber (CU-41L4, Percival Scientific Inc., Dallas County, LA, USA) for cold treatment (~4 • C for 72 h) and short-term recovery (~28 • C). Phenotypes were observed and recorded at cold treatment for 0 h (CK), 1 h (T1), 2 h (T2), 12 h (T3), and rewarming for 1 h after 72 h cold treatment (T4). Leaf samples were collected with liquid nitrogen and stored at −80 • C. The leaf samples, obtained from each group at cold for 0 h (CK), 1 h (T1), 2 h (T2), 12 h (T3), and rewarming for 1 h (T4), were used for the measurements of the content of total soluble sugars (SS), malondialdehyde (MDA), and total chlorophyll (Chl) using assay kits (Solarbio, Beijing, China). According to the instruction, the OD values under the spectrophotometer were operated, measured, and then calculated.

RNA Extraction, cDNA Library Construction and Sequencing
Total RNA was obtained from each group at CK, T1, T2, T3, T4, with a total of 30 samples, and used for RNA-Seq and qRT-PCR analyses. The total RNA sample quality control (QC), library construction, and sequencing on BGISEQ-500 were performed at Beijing Genomics Institute (BGI). The total RNA samples were used for subsequent tests after qualified QC. The Agilent 2100 Bioanalyzer was used to evaluate the integrity of RNA. The mRNA was enriched by magnetic beads with Oligo(dT). The obtained RNA was fragmented by breaking the buffer. Then, the RNA fragments were reverse-transcribed to double-strand cDNA (ds-cDNA) by N6 random primer. The synthesized ds-cDNA was subjected to end-repair and phosphorylated at the 5 end. Then, the sticky ends with an "A" protruding from the 3' end ligate to bulbous adaptors projecting a "T". The PCR products of the ligation products were denatured by heat and the single-strand DNA was cyclized by splint oligo and DNA ligase. Finally, the BGISEQ-500 sequencing platform was used to sequence the DNA library to obtain "raw reads". The raw data were uploaded to NCBI with the accession number PRJNA646356 (https://www.ncbi.nlm.nih.gov/sra/PRJNA646356, Release Date, 15 July 2020).

Sequencing Reads Mapping and Identification of DEGs
Raw reads were statistically analyzed using SOAPnuke and filtered using Trimmomatic [29]. Clean data were obtained by removing low-quality reads, reads containing adapter sequences, and reads with excessive unknown base N content. Then, clean data were aligned to the pepper reference genome sequence (Zunla-1) using HISAT [30]. The gene expression levels were quantified with fragments per kilobase million (FPKM) method using RSEM [31]. Genes with an adjusted p-value (Q-value) ≤ 0.01 and a |log2fold change (FC)| ≥ 1.5 were identified as DEGs by using DESeq2 [32,33].

Functional Annotation and Enrichment Pathway Analyses of DEGs
Mfuzz was used for gene clustering analysis based on loose clustering algorithm [34]. The function of DEGs were annotated using the gene ontology (GO) and KEGG databases, and the enrichment analyses of DEGs were performed using the Phyper function in R software. GO analysis using the Blast2GO software revealed biological significance of the above DEGs.

Gene Coexpression Network Analysis
The weighted gene coexpression network analysis (WGCNA) package was used to perform coexpression network analysis [35]. Gene expression values were imported into WGCNA to construct coexpression modules and to explore the relationship between modules and samples. The networks were visualized using Cytoscape v3.7.1.

Synthesis of cDNA and qRT-PCR
cDNA was synthesized using the all-in-one first-strand cDNA synthesis supermix for qPCR kit (TransGen, Beijing, China). qRT-PCR was completed using the green qPCR supermix UDG kit (TransGen, China) and qPCR instrument Qtower 3 G (Analytik Jena AG, Jena, Germany). To verify the results of RNA-Seq, eight DEGs were selected, including genes related to cold response. The primers for qRT-PCR were listed in Supplementary  Table S1. The relative expression levels of target genes were calculated using the 2 −∆∆Ct method. Pepper gene β-Actin was selected as the housekeeping gene.

Data Analysis
The experimental results were expressed as mean ± standard error and analyzed using Excel 2010 and SPSS 22.0. Significance of differences among different data sets were analyzed using Duncan's multiple range test.

Accession Numbers
RNA-Seq data generated in this study are available from the SRA-Archive (http://www. ncbi.nlm.nih.gov/sra, Release Date, 15 July 2020) with accession number: PRJNA646356.

Different Cold Tolerance in Two Peppers during Cold Treatments
To study the mechanism of pepper resistance to cold stress, two cultivars with different tolerance to cold stress were used for the studies of morphological, physiological, and comparative transcriptomics. Results showed that the leaves of CT having slight water loss and wilting under cold and were slightly damaged after recovery ( Figure 1a). However, the leaves and stalks of CS seriously wilted under cold stress. After rewarming, the stalks stood upright and the leaves died except for the growth point ( Figure 1b). Physiological and biochemical indexes were measured, whose results showed that the change in MDA content is a trend of short-term increase and then followed by a slower decline, and the content in CT is significantly lower than that in CS ( Figure 1c); with the extension of stress time, the content of SS was first increased and then decreased after recovery ( Figure 1c); the change in Chl content was opposite to that of SS content, but the content of SS and Chl of CT were higher than those of CS ( Figure 1c).

Identification of DEGs in CT and CS by RNA-Seq
To investigate genes associated with the response to cold stress and rewarming stage in pepper, the transcriptome profiles in pepper leaves were examined. Using the BGISEQ-500 platform to sequence the RNA of 30 samples, the clean reads produced from each sample varied from 42.31 to 46.71 million. In this project, more than 83.03% reads were aligned to the pepper reference genome, and at least 77.7% reads could be aligned to the unique locations of the reference genome. Additionally, Q30 of each sample exceeded 86.25%, and all GC contents of each sample exceeded 40.07% (Supplementary Table S2). Moreover, total of 30,490 known genes were detected, among which the DEGs were selected on the basis of DEG-seq method with the following parameters: |log2FC| ≥ 1.5 and Q-value ≤ 0.01.
The DEGs specific to each cold-rewarm treatment were studied by comparing the stage with the previous one. The DEGs detected from cold stress treatment (T1, T2, and T3 groups) and rewarming treatment (T4 group) in CT were 1425, 2315, 547, and 6339, respectively ( Figure 2a). In CS, the total DEGs were 258, 322, 921, and 5231, respectively (Figure 2a), which increased following cold exposure. More DEGs were found in T1 and T2 groups and T4 treatment, while fewer DEGs were observed in T3 group were observed. Of particular note, more genes are upregulated than downregulated in T1 group, while more genes are downregulated in T2, T3, and T4 groups ( Figure 2a).  Furthermore, the DEGs of T1, T2, and T3 groups using the starting stage treated for 0 h as a control were also analyzed. Results indicated that the T1 group led to the most DEGs in CT, while the T3 group led to the most DEGs in CS, suggesting that most DEGs regulated preferentially in early stage were associated with greater cold resistance in CT compared with CS. Analysis of differentially expressed genes in the T4 group versus T3 group identified key genes responsible for recovering from cold stress. Venn showed that 34 DEGs were shared by the cold stress, while the number was reduced by only 3 DEGs after recovering treatment in CT (Figure 2c). Hence, these 31 common DEGs could be crucial for the cold-rewarm response, and the 3 DEGs not exhibited in the recovering treatment were responsible for cold stress only. In CS, 42 DEGs were shared in the process of cold stress, and only 14 DEGs of them acted together during the rewarming treatment (Figure 2d), which is less than CT. In total, 7333 and 5953 DEGs responded to one or more of the four treatments in CT and CS, respectively. The 4391 DEGs were shared between the two cultivars (Figure 2b), which were the basic response of peppers to cold treatment. In addition, the 2942 and 1562 DEGs were regulated in only one cultivar, which may be the main reason for the difference in cold resistance of the two cultivars.

Functional Annotation of DEGs under Cold-Rewarming Treatments
The DEGs were significantly (Wallenius test, Q-value ≤ 0.001) enriched to 51 GO terms in comparisons between one or more groups (Supplementary Table S3), which were classified into the "biological process", "cellular component", and "molecular function" categories. The top 10 GO terms in each category, ranked by the number of DEGs, varied with treatment stages, suggesting that the functional responses were different at each stage. Amongst the biological process, the most represented terms including "protein phosphorylation" (GO:0006468) and "regulation of transcription, DNA-templated" (GO:0006355) followed by "transcription, DNA-templated" (GO:0006351) and "defense response" (GO:0006952) were reprogrammed in both CT and CS (Figure 3a). Concerning the cellular component, the major GO terms related to cell part and organelle were enriched upon cold treatment (Figure 3a). For molecular function, the terms with more DEGs were "DNA binding" (GO:0003677), "protein serine/threonine kinase activity" (GO:0004674), and "transcription factor activity" (GO:0003700), showing the importance of these terms in the progress of cold stress response (Figure 3a).
According to Q-value, it was found that the differences in GO terms were significantly enriched between different groups. Compared with CS, the GO terms-such as "regulation of transcription, DNA-templated" (GO:0006355), "photosystem II" (GO:0009523), "transcription factor activity" (GO:0003700)-showed significant enrichment in the early cold treatment of CT. In addition, the GO terms that only enriched significantly in CT included "transcription, DNA-templated" (GO:0006351), "cell redox homeostasis" (GO:0045454), "chloroplast thylakoid membrane" (GO:0009535), etc. The above GO terms may be responsible for the difference in cold tolerance between CT and CS. In the rewarming treatment, it was found that "flavonoid biosynthetic process" (GO:0009813) and "translation" (GO:0006412) showed significant enrichment simultaneously in CT and CS.

Response of Key Unigenes to Photosynthesis, Plant Hormone Signal Transduction, and DDR Processes
To further understand the response of photosynthesis, plant hormone signal transduction, and DDR processes to cold and recovery processes, DEGs related to seven pathways (including carbon fixation in photosynthetic organisms, photosynthesis, photosynthesisantenna proteins, plant hormone signal transduction, DNA replication, HR, and MMR) were annotated. Under cold-rewarm treatments, 52, 21, and 6 DEGs were classified into photosynthesis, photosynthesis-antenna proteins, and carbon fixation in photosynthetic organisms pathways, respectively (Figure 4). In CT and CS, the changes in the regulation of DEGs related to photosynthetic electron transport, photosynthesis-antenna proteins, and carbon fixation in photosynthetic organisms were basically consistent, in which light-harvesting complex II chlorophyll a/b binding protein (LHC) (Capana10g000141 and Capana09g001520), plastocyanin (petE) (Capana04g000129), ferredoxin (petF) (Capana08g001896 and Capana10g001106), and ribulose bisphosphate carboxylase (Rubisco) (Capana02g002781 and Capana02g002784) were firstly upregulated in T1 group; then subsequently downregulated; and finally, significantly decreased in T4 group (Figure 4a,b). Attractively, compared with CS, DEGs involved in photosystem I, photosystem II, Cytochrome b6/f complex, and F-type ATPase were significantly upregulated in T1 and T2 groups of CT (Figure 4a). These DEGs mainly included encoding photosystem I P700 apoprotein (psaA) (Capana00g002019 and Capana07g000048), photosystem II P680 reaction center D2/D1 protein (psbA/psbD) (Capana11g002217/Capana03g001477), photosystem II CP43 chlorophyll apoprotein (psbC) (Capana03g001476 and Capana00g001698), photosystem II CP47 chlorophyll apoprotein (psbB) (Capana00g001154), cytochrome b6 (petB) (Capana01g003214), ATP synthase cf1 (CF1) (Capana00g000237 and Capana03g003995), and ATP synthase cf0 (CF0) (Capana04g000951). Plant hormone signal transduction plays an important role in plant response to cold stress. The DEGs were mainly enriched in auxin, abscisic acid (ABA), and jasmonic acid (JA) pathways ( Figure 5). These might be important hormones related to pepper resistance to cold stress and rewarming recovery. In the two pepper cultivars, most of the genes such as the auxin transporter-like protein (AUX1), transport inhibitor response 1 (TIR1), and auxin-responsive protein IAA (IAA), Auxin response factor (ARF), Gretchen Hagen3 (GH3), and small auxin up RNA (SAUR) had little changes in T2 and T3 group, but their expression levels decreased significantly in T4 group. However, the expression levels of IAA8 (Capana00g000845) and SAUR (Capana03g001923) were significantly increased in T4 group. The expression levels of abscisic acid signal transduction related genes did not change significantly under T2 and T3 group, but the expression levels of most genes includ-ing PYL, PP2C, SnRK2, and ABF increased significantly under T4 group, and the expression levels of PYL (Capana00g003406 and Capana03g002049), ABF (Capana04g000551), and SnRK2 (Capana02g003135) in CT were significantly higher than those in CS. In addition, jasmonate ZIM domain-containing protein (JAZ) (Capana03g000117, Capana09g000144, and Capana07g000750) was upregulated at T3 and T4 groups. Interestingly, under T3 and T4 groups, the expression level of JAZ in CS was significantly higher than that in CT. This suggests that the expression of JAZ may be more sensitive to cold stress. These pathways-including DNA replication, HR, and MMR-were significantly enriched, which belonged to DDR. Replication protein a (RPA) (Capana03g000868 and Capana00g004882) and DNA polymerase delta (POL δ) (Capana06g002135) were involved in the above three pathways ( Figure 6). Compared with CS, most DEGs were significantly upregulated in T2 and T3 groups of CT; these DEGs showed lower expressions in T1 and T4 groups of CT, while the expressions in CS were higher than those in CT. It indicated that a large number of genes were upregulated during DNA damage repair after CT was subjected to cold stress, thus conducting DNA repair of cold stress injuries.

Classification of Transcription Factors (TFs) Associated with Cold Treatments
The TFs of 57 families were included in the DEGs of CT and CS during cold-rewarming (Supplementary Table S5). The largest group of TFs was the ethylene-responsive factor (ERF) family, followed by the basic helix-loop-helix (bHLH), whereas other TFs belonged to the MYB, WRKY, NAC, C2H2, and GRAS families (Figure 7a). The number of TFs such as NAC, B3, and HD-Zip in CT were much more than those in CS, but they were the opposite for C2H2 and Dof. In the NAC families, the expression levels of 13 NAC transcription factors remained low in the T1, T2, and T3 groups, but increased significantly in T4 group (Figure 7b). Among the 13 NAC transcription factors, the expression levels of NAC8 (Capana02g003557) and NAC72 (Capana07g002219) in CT were significantly higher than those in CS under T4 group.

Gene Coexpression Network Analysis
The gene coexpression profiles of pepper in response to cold stress and rewarming conditions were analyzed using WGCNA. Seventeen coexpression modules and correlation coefficients were identified (Figure 8). The resistance of cold was positively correlated with the MEturquoise, and the correlation coefficient was 0.85. The MEyellow (0.76), MEmidnightblue (0.97), and MEblue (0.91) modules had higher correlation coefficients at T1, T3, and T4, respectively. Further, it was found that CT-T1 was more positively correlated with the MEtan (0.97). The genes related to cold stress and recovery of pepper were mainly concentrated in the MEturquoise, MEmidnightblue, MEblue, and MEtan modules. Therefore, the hub genes in each module were selected for further analysis and drawing the network (Figure 9). The genes, including Capana08g000846, Capana10g001816, Capana10g000688, Capana07g001882, Capana03g003720, Capana04g000027, Capana04g002349, and Capana11g000347, which encoded translation initiation factor 3 subunit F (EIF3F), zinc finger protein, two-component response regulator-like APRR9, cytochrome P450 family 4 subfamily X (CYP), CCR4-associated factor 1 (CAF1), calmodulin-lysine n-methyltransferase (CAMKMT), bzip transcription factor, and heat shock 70 kda protein 8 (HSP70), respectively, had higher indegrees in four modules (Figure 9 and Supplementary Table S6).

qRT-PCR Validation of Gene Expression Patterns
In order to validate the results of RNA-Seq, qRT-PCR analyses were performed on 8 selected genes using gene-specific primers. qRT-PCR analysis showed similar trends of transcript abundance when assessed by real-time RT-PCR ( Figure 10), thus confirming that the FPKM valve of each transcript determined by RNA-Seq was reliable in this study.

Cold Stress Affects the Physiology and Phenotype of Pepper Cultivars
As plants suffer cold stress, a series of physiology and biochemistry processes can be induced to reduce the cold damage to them [36]. Under cold stress, excessive ROS will be generated among plant cells. Some osmotic adjustment substances such as SS and proline (Pro), acting as cryoprotectants, can against freezing-dehydration damage and remove ROS [37,38]. Furthermore, the accumulation of ROS will degrade polyunsaturated lipids to form MDA, which causes cellular dysfunction and tissue damage and is the most suitable index to determine the degree of oxidative damage of plant biofilms [39,40]. In this study, the content of SS in CT was significantly higher than that in CS, while the content of MDA in CT was significantly lower than that in CS in T1, T2, and T3 treatments, indicating that the cold damage of CS was more serious. Moreover, the chlorophyll content significantly decreased in T3 treatments both in CT and CS cultivars. This is because cold stress could slow down the plant metabolism and reduce the raw materials used for chlorophyll synthesis, resulting in the reduction of chlorophyll content [41].

Response of Photosynthesis and Related Processes to Cold and Rewarming Conditions
In addition to changes in chlorophyll content, transcriptome data showed that the DEGs related to photosynthesis were significantly enriched in pepper during cold stress and recovery stage. Photosynthesis is one of the most temperature-sensitive processes in plants. Plants could adjust their photosynthetic characteristics to adapt to cold [42]. Previous studies have shown that the reduction of photosynthesis is parallel to a strong decline of PSII photochemistry in various plant species exposed to cold, and suggested PSII as one of the primary targets of cold stress [43,44]. Arabidopsis were exposed to a relatively moderate cold of 12 • C, which caused the reduction of PSII reaction center polypeptide psbA [45]. In this study, genes in the PSII system including psbA (Capana11g002217), psbB (Capana00g001154), psbC (Capana03g001476 and Capana00g001698), and psbD (Capana03g001477) were downregulated as cold stress treatment continued, which could affect the photochemical reaction of PSII. In addition to PSII, exposure of various plant species to cold stress can also cause photoinhibition damage in the PSI reaction center [43,46]. The lower amount of psaB might be also due to downregulation of the genes and/or degradation of PSI to response to the cold stress [47]. Besides, psaA is a core protein and crucial for the functional assembly of PSI [48]. Except for its significant increase in the T1 group, the downregulation of psaA (Capana00g002019 and Capana07g000048) under cold stress were also found in this study. However, it seems to be more meaningful to maintain the stability of PSI reaction center complex rather than PSII, which has a much faster repair cycle, as PSI may be the limiting factor for the recovery of photosynthesis in rewarming stage [45]. In T4 group, the above genes were all downregulated. These results indicate that the gene expression of PSI and PSII is difficult to return to normal level after suffered cold stress.
Furthermore, as an electron transfer enzyme, the Cytochrome b6/f complex can be involved in electron transport to create proton gradient, which drives ATP synthesis in chloroplasts [49,50]; petA, petB, and petD are the major chloroplast pet genes that encode the cytochrome b6/f complex [51]. In this study, significant changes of petB (Capana01g003214) occurred only in CT cultivars. Previous studies have shown various LHC compounds involved in several chloroplast functions, including photosynthetic pigment metabolism and photosynthetic electron transport chain, to adapt cold stress [52]. The increase of LHC is beneficial to adjust the organization of photosystem II and photoprotection. In this study, the expression levels of LHCB4 (Capana09g001520) and LHCA3 (Capana10g000141) in CT were higher than those in CS under T1 and T2 groups, and the LHCB1 (Capana08g000679) was only expressed in CT cultivar. The specific expression of these genes under low temperature stress may be an important cause of different cold tolerance. Rubisco also continued to increase after short-term cold stimulus (T1 group); cold increases the affinity of Rubisco for CO 2 relative to O 2 and the solubility of CO 2 , which is crucial for C3 photosynthesis plants [53].

Plant Hormones Reflect Cold and Rewarming Condition of Pepper
Plants can quickly sense and respond to stress through hormone-mediated signal transduction pathways. Plant hormone signal transduction plays a central role in response to cold stress, especially the role of ABA in sensing and receiving cold signals, which triggers the corresponding responses of plants to adapt to abiotic stress [54]. Under cold stress, ABA content in plants increased, and ABA entered the central hydrophobic site of PYL and induced the closure of related structures of PYL, creating conditions for the binding of PP2Cs. Then, the tryptophan residues in PP2C were inserted into the ABA structure and tightly locked ABA in the appropriate position [55,56]. In the PYL-ABA-PP2C complex, the binding of ABA-PYL with PP2Cs blocked the binding of SnRK2 with PP2Cs to prevent their interaction. The released SnRK2s were self-activated through phosphorylation effect proteins including TF (ABF, etc.) and downstream stress response genes, thus improving the tolerance of plants to the low-temperature environment [57]. It has been reported that SnRK2.6 (OST1) activated by cold stress phosphorylates ICE1 to activate CBF-COR gene expression cascade and cold tolerance [58]. In this study, PYL (Capana00g003406 and Capana03g002049), ABF (Capana04g000551), and SnRK2 (Capana02g003135) were all upregulated in T4 treatments and the expression levels of these genes in CT were significantly higher than those in CS. This suggests that the high expression of these genes may contribute to the recovery of pepper seedlings after cold stress.
JAZ also interacts with ICE1 to inhibit the ICE-CBF signaling pathway. Cold induces the expression of JA synthesis-related genes and the synthesis of bioactive JA-Ile, which activates JA receptor COI1 to bind to JAZ1, resulting in the degradation of JAZ1 through 26S proteome after ubiquitination. Then, the ICE-CBF transcriptional regulation cascade signaling pathway is activated, and the expressions of cold-regulated genes are induced to improve plant cold tolerance [59]. The downregulation of JAZ can also reduce the binding with ICE to regulate cold stress. In this study, under T3 and T4 groups, the expression level of JAZ in CS was significantly higher than that in CT. This may be an important reason why cold-resistant cultivars recover more easily than cold-sensitive cultivars. Further, there are increasing evidences that SAUR genes, a family of auxin-responsive genes, are involved in stress/defense responses [60]. Under cold treatment, 8 SAUR genes were upregulated, which indicated they may play an important role when the poplar is subjected to cold stimulation [61]. In this study, three SAUR genes (Capana03g001923, Capana09g001771 and Capana01g003581) were upregulated in T4 group, this suggested that SAUR may play an important regulatory role in recovery stage after cold stress in peppers.

Response of Defence Mechanism to Cold and Rewarming Conditions
Maintaining genome integrity is essential for an organism to reach its full life span. In order to maintain genome integrity and prevent DNA lesion transmission to daughter cells, DNA damage is recognized by a variety of sensors, including MutS, the RPA complex and the MRE11-RAD50-NBS1 complex [62]. These sensors activate DDR mechanisms, specifically evolved by plants, including direct repair, nucleotide excision repair, base excision repair, nonhomologous end joining, MMR, and HR [63]. This study focused on MMR-and HR-related genes; most DEGs in MMR and HR were significantly upregulated in T2 and T3 groups of CT compared with CS, suggesting that cold tolerance cultivars may have a strong DNA damage repair ability when suffering cold stress, while cold sensitive cultivars have a weak DNA repair ability. Mismatch repair corrects mismatches caused by misincorporation of nucleotides during DNA replication and recombination. After identifying the lesion, MutS recruits MutL [64], and PCNA loaded by replication factor C interacts with MutLa to enable PMS1 to exert its endonuclease activity. The endonuclease activity of MutLa generates new entry sites for the 5 -3 ExoI, which removes the nicked strand [65,66]. Finally, POL δ and DNA ligase I complete the repair process [67]; HR is an important repair mechanism used to eliminate DNA double-strand breaks. The MRE11-RAD50-NBS1 complex recognizes DNA damage. Following resection by MRE11, the resulting single stranded DNA is bound by RPA and then replaced by RAD51 [68].
The formed RAD51 nucleoprotein filament facilitates DNA strand invasion and exchange, resulting in the formation of a D-loop [69]. Then, DNA synthesis is conducted by three mechanisms, including non-crossover double-strand break repair, synthesis-dependent strand annealing, and break-induced replication. The DRR mechanism of plants was significantly enriched in this study and showed great differences between the two cultivars, but few studies have been conducted on it at cold stress, so the specific mechanism of action in DDR remains to be further studied.

TFs Regulate the Cold Resistance of Pepper
Transcription regulation is a key step in plant response to abiotic stress by regulating their target genes of transcription [70]. TFs can improve the stress tolerance of plants by regulating downstream genes. Most TFs related to cold tolerance in pepper have been reported, including ERF, bHLH, MYB, WRKY, NAC, GRAS, etc. [71][72][73][74]. Li et al.'s [75] results suggested that SlNAC1 functions as a stress-responsive NAC protein involved in the abscisic acid-dependent signaling pathway and may have potential applications in transgenic breeding to enhance crops' cold stress tolerances. In Sweet Osmanthus, the expression patterns of most branch '5B' members, especially OfNAC49 and 59, were significantly induced by cold stress, indicating that branch '5B' members could have important roles in cold tolerance in O. fragrans [76]. CaNAC064-silenced pepper plants exhibited more serious wilting, higher MDA contents, and chilling injury index in the leaves after cold stress. The CaNAC064-overexpressing Arabidopsis plants exhibited lower MDA content and chilling injury index compared with WT plants under cold stress [77]. In this study, the expression levels of 13 NAC transcription factors increased significantly in T4 group; meanwhile, the NAC8 (Capana02g003557) and NAC72 (Capana07g002219) in CT were significantly higher than that in CS under T4 group. These results may indicate that NAC plays an important role in regulating the recovery of pepper after cold stress, and the expression levels of NAC8 and NAC72 determine its recovery ability.

Network Analysis Revealed the Key Genes of Cold Resistance of Pepper
In this study, WGCNA analysis revealed eight hub genes in modules after the four key modules related to cold resistance and different stages of treatment of pepper were identified. These genes were divided into three groups. The first group was encoding and regulating genes, related to protein synthesis and modification, including the bzip transcription factor, CAF1, EIF3F, CAMKMT, and APRR9. This kind of research is relatively clear and is also the most direct response to stress. The second is coding effector proteins, such as zinc finger protein and HSP, which can directly protect plant cells from stress. Overexpression of SlCZFP1 from tomato in transgenic Arabidopsis and rice induced constitutive expression of COR genes and conferred enhanced tolerance to cold treatments for transgenic plants [78]. HSPs, as a molecular chaperone, can regulate the correct folding of proteins, which is a key protein to maintain cell homeostasis under stress. Cold stress resulted in upregulation of Hsp70 and downregulation of Hsp90 in winter wheat, which contributed to plant phenotypic variation in plant response to stress [79]. The third encodes enzymes, such as CYP, which are involved in catalysis of variety of reactions that include growth, development, and secondary metabolite biosynthetic pathways. CBF/DREB1 transcription factor regulates the expression of SLCYP72A184, SLCYP85A1, and SLCYP96A48 genes, which can be used as candidate genes for cold tolerance in tomato [80].

Conclusions
In this study, CT and CS with different cold tolerances were used to study the coldrewarming-related mechanisms of pepper. The physiological indicators, transcriptome, and WGCNA analysis provide a great deal of information for important pathways and hub genes related to cold-rewarm treatments. These pathways related to photosynthesis, plant hormone signal transduction, and DNA damage repair were significantly enriched. Importantly, the high accumulation capacity of common hub genes and important pathways-related genes may contribute to the cold tolerance of CT. Our results provide a better understanding of the extremely complex regulatory mechanisms in the process of cold-rewarm responses in pepper.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/app112110204/s1, Table S1: The primers used in this study, Table S2: Summary of RNA-seq data quality, Table S3: The GO merge of DEGs, Table S4: The KEGG merge of DEGs, Table S5: The numbers of TFs family, Table S6: Cytoscapelnput-edges.