Two Arachis hypogaea cultivars (Silihong and Zhonghua 12) with differential cold tolerance were used in this study (Wang et al., 2021; Zhang, et al., 2021). Peanut seeds were soaked in distilled water for 6 h, potted in soil, and then placed into a growth chamber under a 16 h/8 h (light/dark) cycle at 28°C. Two-week-old seedlings were subjected to cold stress (10°C). The third unfold leaves from seedlings were collected after 0 h, 3 h, 12 h, and 24 h. All samples with three biological replicates were quickly frozen in liquid nitrogen and stored at −80°C until use.

Total RNA was extracted from peanut leaves treated with and without cold stress treatment. RNA integrity and concentrations were detected by the Bioanalyzer 2100 system (Agilent Technologies, CA, United States). One microgram of total RNA for each sample was prepared for cDNA library according to the protocol provided by Oxford Nanopore Technologies (ONT). The cDNA libraries were added to flowcells (FLO-MIN109) and run on an ONT sequencing platform (PromethION 24, P24) at Biomarker Technology Company (Beijing, China). The raw reads derived from 24 samples have been deposited in the NCBI SRA database with accession number of PRJNA1032881. Low quality (length <500 bp, Q-score <7) of raw reads and ribosomal RNA sequences were first filtered. Then, the full-length transcripts were determined by searching for clean reads that have ONT library primers at both ends. Full-length transcripts were obtained by mapping to reference genome (cv. Tifrunner, https://data.legumeinfo.org/Arachis/hypogaea/genomes/Tifrunner.gnm2.J5K5, accessed on April 4, 2023) with mimimap2 software. Mapped reads were further collapsed by cDNA_Cupcake package with minimum coverage of 85% and minimum identity of 90%. The duplicated transcripts were removed.

The obtained full-length transcripts were annotated with the following databases: NR (NCBI non-reductant protein sequences), Swiss-port database, Pfam (protein family), Gene Ontology (GO), and KEGG (Kyoto Encyclopedia of Genes and Genomes). GO enrichment analysis was performed by the GOseq R packages (Young et al., 2010). KEGG enrichment analysis was implemented by KOBAS software (Mao et al., 2005).

Transcripts were validated against available reference transcriptome annotations (https://www.peanutbase.org/download/) with gffcompare. Expression levels were estimated by reads per gene/transcript per 10,000 reads mapped (CPM). DEGs were screened using the DESeq2 R package (1.6.3) with the criteria of |log2FC (fold change)| > 1 and an adjusted p-value < 0.01 (Love et al., 2014). Alternative splicing (AS) events were analyzed by the AStalavista tool (Foissac and Sammeth, 2007), and DAS events were identified by PSI-sigma tools (Lin and Krainer, 2019). Transcription factors were detected from the transcriptome data using iTAK software (Zheng et al., 2016).

Approximately 1 µg of total RNA (DNase I-treated) was used for first-strand cDNA synthesis for each sample with MMLV reverse transcriptase kit (Thermofisher Scientific, United States). RT-PCR was used for analysis of the alternatively spliced isoforms, and real-time quantitative RT-PCR (qRT-PCR) was performed to analyze the expression level of DEGs. qRT-PCR reactions were carried out on the Bio-Rad CFX96 RT-PCR Detection system (Bio-Rad, Hercules, CA, United States) using Hieff qPCR SYBR Green Master Mix. The relative gene expression levels were quantified by 2^—△△Ct method (Livak and Schmittgen, 2001), using A. hypogaea Actin gene (accession number Aradu.W2Y55) as an internal reference gene for normalization. The expression level of 0 h in SLH was used as a control, whose value was set to 1. All primers used are listed in Supplementary Table S1.

The coding sequences of AhTPP genes (Arahy.FIK9JS, Arahy.0FY2NM) were amplified by PCR using peanut leaf cDNA as template with primers listed in Supplementary Table S1. The PCR products were purified and cloned into yeast expression vector (pYES2). The resulting pYES2-AhTPPs plasmids were transformed into yeast cells (INVSc1) individually. For cold tolerance assays, the transgenic yeast cultures harboring empty vector pYES2 (control) or pYES-TPPs were first subjected to −20°C for different times (0 h, 24 h, 48 h, 72 h, and 96 h) and serially diluted with SG-U liquid medium (synthetic complete without uracil medium + 2% galactose + 2% agar) to different concentrations (10⁰, 10⁻¹, 10⁻², and 10⁻³); then, 2 μL of each culture was spotted onto SG-U medium plates. The plates were photographed after incubating at 28°C for 72 h.

Peanut cultivar Silihong (SLH) is popularly planted in northeast China. Zhonghua 12 (ZH12) is a representative variety released by Oil Crops Research Institute, Chinese Academy of Agricultural Sciences. In a previous study, SLH was shown to have better cold tolerance than ZH12 (Zhang et al., 2021). First, both the two cultivars showed similar seed germination rates (>90%) under normal condition. The germination rate was severely down to approximately 50% in ZH12 (cold susceptible) after 72 h cold treatment, but it was not affected in SLH (cold tolerant). Second, the levels of reactive oxygen species (H₂O₂ and O₂⁻) and contents of malondialdehyde (MDA, an important indicator for cell membrane damage caused by oxidative stress) were significantly lower in the seedlings of SLH than those in ZH12 under cold stress. Thus, the differential performance and physiological response to low temperature indicated that SLH was relatively more cold-resistant than ZH12 (Wang et al., 2021).

The long-read ONT sequencing platform could theoretically obtain the full-length transcript information without requiring fragmented RNA. To investigate the molecular basis of peanut in response to cold response, we use ONT sequencing to analyze the transcriptome of SLH and ZH12 after cold treatment for different time (0 h, 3 h, 12 h, and 24 h). As a result, a total of 168 million clean reads derived from the 24 samples (three biological replicates per time points) were obtained, with an average length of 1,040–1,393 nucleotides per each sample (Supplementary Table S2). The number of full-length transcripts varied from 4.5 to 8.6 million per library, accounting for more than 84% of clean reads (Supplementary Table S3), indicating that the data were good enough for following bioinformatic analysis.

Full-length transcriptome based on ONT sequencing could facilitate to optimize the structure of transcripts by alignment of long-read transcripts with the reference genome and/or transcriptome of peanut cv. Tifrunner (Bertioli et al., 2019). When mapped to the genome sequences, a total of 38,511 genes were discovered in our transcriptome data, including 8,949 novel gene loci (Supplementary Table S4). Some novel genes with high expression levels were found to be associated with stress response. For example, ONT.16282 was annotated as a homolog of abscisic stress-ripening protein 2 (ASR2), which confers tolerance to drought in plants (Srivastava et al., 2022). In addition, 125,159 redundant-removed consensus transcript sequences were obtained through minimap2 software. Among of them, only 29,958 transcripts (23.93% of total transcripts) were known isoforms, while 81,364 (65.0%) were new isoforms mapped to the known genes, and 13,837 (11.05%) were considered as novel isoforms that mapped to unannotated genomic regions (Supplementary Table S5). The results indicated that ONT sequencing could be used as an efficient and robust tool to provide significant additional annotations to the available reference genome.

The cold-responsive genes were identified by screening DEGs between cold-treated samples (3 h, 12 h, and 24 h) and normal control (0 h). In both cultivars, the number of DEGs in 12 h or 24 h was greater than that in 3 h (Figure 1A), suggesting that prolonged duration of cold treatment led to more profound alternations of transcriptome profiles in peanut. There was a total of 10,677 and 16,363 DEGs in the three comparison sets (24 h vs. 0 h, 12 h vs. 0 h, 3 h vs. 0 h) of SLH and ZH12 respectively (Supplementary Table S6), and 7,935 common DEGs were shared between the two cultivars (Figure 1B). We performed GO analysis of these common DEGs to explore their possible functional processes. The most enriched GO terms associated with “biological process” categories was “response to stress” (GO:0006950) (Figure 1C; Supplementary Table S7). Most stress-associated genes were significantly upregulated by low temperature in this GO term (Figure 1E; Supplementary Table S8), including glutathione S-transferases (Arahy.F8LB2J, Arahy.67E0CQ), suggesting their vital roles in peanut cold resistance. The top frequent GO terms belonging to “cellular component” and “molecular function” categories were “chloroplast part” (GO:0044434) and “transcription factor activity” (GO:0003700), respectively, indicating that photosynthesis process and transcription factors play vital roles in cold stress response in peanut. KEGG enrichment analysis showed that “circadian rhythm-plant” (ko04712) and “glutathione metabolism” (ko00480) were most abundantly enriched (Figure 1D). Moreover, 2,742 DEGs were found specifically in SLH. “Transketolase activity” (GO:0004802) was the most enriched GO term, and biosynthesis of amino acids (ko01230) was the top enriched pathway (Supplementary Figure S1). This result suggested that the pentose phosphate pathway and protein metabolisms were more associated with cold response in resistant peanut variety (SLH). To verify the reliability of the expression profiles of DEGs from ONT sequencing data, 10 DEGs were selected for qRT-PCR analysis. As a result, the transcript levels obtained by qRT-PCR showed similar expression trends to the corresponding genes in RNA-seq data, confirming the reliability of transcriptome analysis (Figure 2).

Since transcription factors (TFs) played a crucial role in responding to cold stress, we performed genome-wide analysis of genes containing functional domains of TFs family in the two peanut cultivars. In total, 4,493 TFs were identified in our transcriptome data (Supplementary Table S9). The largest number of TF family was the bHLH member, followed by AP2/ethylene-responsive transcription (AP2/ERF) and MYB (Supplementary Figure S2). A total of 1,350 cold-responsive TFs were discovered in the two peanut cultivars (Supplementary Table S10). Among of them, 526 specific TFs had differential expression levels after cold treatment in ZH12, while only 193 specific TFs in SLH (Figure 3A). These 1,350 cold-responsive TFs were most significantly enriched in “plant hormone signal transduction” (ko04075) and “circadian rhythm-plant” (ko04712) KEGG pathways (Figure 3B). The AP2/ERF members accounted for the largest proportion of cold-responsive TFs (Figure 3C). CBFs/DREBs, belonging to AP2/ERF family, were extensively studied TFs that have important roles in plant cold acclimation (Kim et al., 2015; Tang et al., 2020). We identified 28 AhDREBs in the transcriptome (Supplementary Table S11), including AhDREB3 (Arahy.SE8WTS), AhDREB27 (Arahy.KQL0FY), AhDREB35 (Arahy.5KI8FH), and AhDREB39 (Arahy.L5HMHT). Most of them were significantly upregulated under cold stress in both SLH and ZH12 (Figure 3D). In addition, the 193 SLH-specific TFs were most enriched in “plant hormone signal transduction” KEGG pathway, with AP2/ERF members being the largest (Supplementary Figure S3). These results indicated that AP2/ERF play vital roles in signal transduction in adapting peanut to cold stress.

As mentioned above, we identified 95,201 novel isoforms in our transcriptome (Supplementary Table S5). Most of these novel isoforms were probably derived from different forms of alternative splicing (AS) and extensions or truncations at 5′ or 3′ ends of mRNA. Polyadenylated (PolyA) tail, a long chain of adenine nucleotides, is added to the 3′ untranslated region (3′ UTR) of mRNA molecules for regulating their stability and translation (Jia et al., 2021). Due to presence of two or more poly A signals within the transcript, alternative polyadenylation (APA) events often occur to generate more isoforms and increase the diversity of transcriptome (Tu et al., 2020). We performed analysis of the characteristics of APA sites in the two peanut cultivars (Figure 4A; Supplementary Table S12). It was found that more than SLH (21,568) had more genes with APA sites than ZH12 (20,956) under normal condition (0 h). After cold treatment, the number of gene that possessed more than five APA sites was significantly increased in SLH, while other types of APA events were decreased in this cultivar. In contrast, the number of gene with one APA site was raised in ZH12, and the numbers of the remaining types of APA events were down. This result indicated that the cold stress promoted more APA sites that occurred in SLH (cold resistant) but that was inhibited in ZH12 (cold sensitive). In addition, the top 10 putative polyadenylation motifs were identified in the transcriptome data. The APA sites were most enriched in the motif WTGKA, accounting for 19.48% in the top 10 motifs (Figure 4B; Supplementary Table S13).

AS is an important posttranscriptional process of gene expression that affects the stability of transcripts and increases the diversity of proteins. A total of 23,673 (60.7%) genes were found to possess two or more alternative splicing variants in this transcriptome (Figure 5A; Supplementary Table S14). The number of AS events was significantly increased in SLH, while it was decreased in ZH12 under cold conditions (Figure 5B). This result indicated the differential impact of cold stress on AS events in the two cultivars. Moreover, five major types of AS events were discovered in the samples, including exon skipping (ES), intron retention (IR), mutually exclusive exon (MXE), alternative 5′ splice sites (A5SS), and alternative 3′ splice sites (A3SS). IR was the most abundant type of AS events and accounted for more than 40% of all AS events in SLH and ZH12 (Figure 5C; Supplementary Table S15). In both cultivars, after cold treatment, the proportion of A3SS and IR events decreased while ES increased compared to their normal controls (0 h). We performed RT-PCR analysis of six selected genes, including DREBs and splicing factors, to validate the cold-responsive alternative splicing events. The results showed that the alternative splicing patterns were consistent with their profiles revealed by ONT sequencing. As shown in Figure 5D, AhDREB3 had two AS transcripts. The larger size of amplicon was AhDREB3.1 (706 bp), which possessed an IR event in the coding region of the primary transcript, while the smaller one corresponded to AhDREB3.2 (196 bp). Similar results were also observed in the other four genes (AhDREB35, AhSP41, AhSP35, and AhRNApol). On the other hand, AhERF underwent an A5SS event in the second exon of the primary mRNA, which resulted in two AS transcripts. The larger size of amplicon was AhERF.1 (524 bp), and the smaller one corresponded to AhERF.2 (120 bp). There was no positive signal in the second lane for AhRNApol gene, which was probably due to the fact that the transcript level of AhRNApol was substantially decreased after 24 h cold treatment in SLH.

DAS events were detected between cold-treated samples (3 h, 12 h, and 24 h) and normal control (0 h). There was a total of 402 differential alternative splicing genes (DASGs) in the three comparison sets (24 h vs. 0 h, 12 h vs. 0 h, 3 h vs. 0 h) of SLH and 344 DASGs in ZH12, respectively (Supplementary Table S16). Unlike DEGs, the counts of gene with DAS events greatly decreased in ZH12 but slightly declined in SLH after 24 h of cold treatment (Figure 6A). To investigate the effect of cold-induced alternative splicing on cellular processes, we performed functional GO category and KEGG enrichment analysis of DASGs in the two cultivars. The DASGs in SLH had different enriched GO categories and KEGG pathways with those in ZH12 (Figures 6B–F). For instance, “nitrogen metabolism” and “ubiquitin mediated proteolysis” were the most enriched pathways in SLH, while “sulfur metabolism” and “ubiquinone and other terpenoid–quinone biosynthesis” were top enriched in ZH12. Next, a comparison was made between DASGs and DEGs to investigate how many DASGs were also changed their transcriptional levels under cold stress (Figure 6G). A total of 295 DASGs were differentially expressed in response to cold stress, while 427 DASGs did not have altered transcript levels. Moreover, we identified 142 DEGs with DAS events in SLH and 165 DEGs with DAS events in ZH12 (Supplementary Tables S17, S18). It has been shown that cold-induced DAS events often led to the formation of premature termination codons (PTCs) (Shen et al., 2016; Li et al., 2020a). Most of DEGs with DAS events were downregulated or degraded under cold stress in both of the two cultivars (Figures 6H, I).

In our previous study, we have developed a highly efficient yeast functional system to characterize potential cold tolerance genes in peanut (Wang et al., 2023). Based on our ONT transcriptome data, we selected two candidate genes (Arahy.FIK9JS and Arahy.0FY2NM) to verify their functions using the yeast platform. Both of the genes were significantly upregulated by low temperature in the two peanut cultivars and encoded putative trehalose-6-phosphate phosphatases (AhTPPs) (Supplementary Table S6). The two AhTPPs were transformed into yeast cells individually, and their cold tolerance were evaluated after low temperature treatment. As shown in Figure 7, the yeast transformants overexpressing AhTPPs exhibited better growth than the cells harboring pYES2 empty vector (control). Specifically, after cold treatment for 96 h, the control yeast cannot survive on the solid medium when the cultures were diluted to 10⁻³, whereas the transgenic yeast cells containing AhTPP were still able to grow under the same stressful conditions. These results indicated that overexpression of AhTPP improves cold tolerance in transgenic yeast.