Integration of chromatin accessibility and gene expression reveals new regulators of cold hardening to enhance freezing tolerance in Prunus mume

Abstract

Low temperature is one of the most important abiotic factors limiting the growth, development and geographical distribution of plants. Prunus mume is an attractive woody ornamental plant that blooms in early spring in Beijing. However, the molecular mechanisms underlying cold hardening to enhance freezing tolerance in Prunus genus remains elusive. This study examined the dynamic physiological responses induced by cold hardening, and identified freezing-tolerance genes by RNA-seq and ATAC-seq analyses. Cold hardening elevated the content of soluble substances and enhanced freezing resistance in P. mume. Transcriptome analysis indicated that the candidate differentially expressed genes (DEGs) were those enriched in Ca²⁺ signalling, mitogen-activated protein kinase (MAPK) cascade, abscisic acid signalling, and inducer of CBF expression 1 (ICE)-C-repeat binding factor (CBF) signalling pathways. The openness of gene chromatin positively correlated with the expression level of these genes. Thirteen motifs were identified in the open chromatin regions in the treatment group subjected to freezing after cold hardening. The chromatin opening of transcription start site at the proximal –177 region of cold-shock protein CS120-like (PmCSL) was markedly increased, while the expression level of PmCSL was significantly up-regulated. Overexpression of PmCSL in Arabidopsis significantly improved the freezing tolerance of transgenic plants. These findings provide new insights into the regulatory mechanism of freezing tolerance to improve breeding of cold-hardy P. mume plants.

Introduction

Low temperature stress not only affects plant growth and development, but also restricts the geographical distribution of plants (Zhu, 2016). Cold hardening to enhance freezing tolerance is an important strategy for plants to cope with low temperature stress (Miki et al., 2019; Ambroise et al., 2020). In the process of cold hardening, plants synthesize a large number of protective substances to enhance their cold tolerance (Rinne et al., 1994; Cook et al., 2004; Kaplan et al., 2007). Accumulation of soluble substances protects plant cells by regulating the osmotic potential of cells and reducing water loss under freezing conditions (Strand et al., 1999; Takagi et al., 2003; Kaplan and Guy, 2004; Rohde et al., 2004). Low temperature stress increases ROS (reactive oxygen species) content in plants, and excessive ROS can be toxic to cells (Lyons and Raison, 1970). Cold hardening can significantly improve the activities of peroxidase, ascorbate peroxidase, and glutathione reductase, promote the synthesis of antioxidants, and alleviate the oxidative stress damage caused by freezing stress (Lyons and Raison, 1970; Einset et al., 2007). Various antioxidant enzymes play synergistic functions in different biological processes to improve the low temperature tolerance of plants (Dai et al., 2009).

Plant response to low temperature involves a precise and complex regulatory mechanism of gene expression, through long-term evolution (Guy et al., 1985; Thomashow, 1999; Kosmala et al., 2009; Winfield et al., 2010). Many studies have shown that plants respond to low temperature through multi-gene expression (Guy et al., 1985; Thomashow, 1999; Winfield et al., 2010). In A. thaliana, the differentially expressed genes (DEGs) after 24 h of cold treatment account for 4% of all genes in the genome, while the DEGs after 14 d of cold treatment account for 20% of the genome (Hannah et al., 2005; Lee et al., 2005). In nature, freezing-tolerant species or varieties undergo cold hardening earlier in response to cold winters, compared with freezing-sensitive plants (Winfield et al., 2010). Meanwhile, more DEGs are expressed in freezing-tolerant, rather than freezing-sensitive varieties, under freezing stress (Thomashow, 1999; Liang et al., 2020). Most plants are exposed to unfreezing low temperature environments, and gradually develop freezing resistance through the activity of cold response genes (Örvar et al., 2000; Guo et al., 2018). The expression of cold response genes is not only related to the time of low temperature induction, but also closely related to the temperature threshold (Fowler, 2008).

Genomic DNA strands in eukaryotes are usually tightly intertwined with histones to form nucleosomes, which are folded and compressed to form higher chromosomal structures (Fyodorov et al., 2018). Normally, the chromatin in the closed state does not replicate and transcribe. The transformation of chromatin from inhibited state to active state is crucial for plants to survive a constantly changing environment (Park et al., 2018). Open chromatin allows the binding of regulatory factors at transcriptional levels. Therefore, the degree of open chromatin is closely related to the gene expression level. The assay for transposase-accessible chromatin (ATAC)-seq is a high-throughput sequencing technology that uses transposase to study chromatin accessibility (Buenrostro et al., 2015). The application of ATAC-seq to delineate open chromatin regions across each genome has been gradually developed in plant species, such as Arabidopsis thaliana, Oryza sativa, Solanum lycopersicum, and Medicago truncatula (Maher et al., 2018). The chromatin accessibility and bivalent histone modification of active genes are achieved in Solanum tuberosum under cold conditions. Chromatin accessibility and translational landscapes of Camellia sinensis under chilling stress have been examined by ATAC-seq (Wang et al., 2021). Under specific conditions, the chromatin accessibility landscape can help us investigate the transcriptional regulatory mechanisms in the whole genome by detecting the range of protein binding sites (Zeng et al., 2019).

Prunus mume, one of the most traditional flowers in China, has important economic, ornamental, and cultural values. The natural distribution centre of P. mume lies in the Hengduan mountains of Sichuan, Yunnan, and Tibet (Zhang et al., 2012, 2018). The Yangtze River Basin is the main birthplace of the diversity of P. mume (Zhang et al., 2012, 2018). After more than 60 years of research on freezing resistance breeding, a few P. mume varieties have successfully achieved winter hardiness and blossomed in Beijing. Compared with annual plants, perennial woody plants have more complex and diverse cold resistance mechanisms, including cold hardening and winter dormancy (Hincha and Zuther, 2014; Tylewicz et al., 2018). Studies on the molecular mechanism of improving freezing tolerance by cold hardening are helpful to breed excellent freezing-resistant varieties. Here, we examined the changes in physiological responses induced by cold hardening of P. mume varieties with varying cold resistance, and determined the key stages for cold hardening. Integrated analysis of RNA-seq and ATAC-seq was performed to explore the genes related to freezing-tolerance after cold hardening, and a potential molecular regulatory mechanism was proposed for P. mume. Furthermore, the functions and molecular regulation mechanisms of hub genes were also verified in Arabidopsis.

Materials and methods

Plant materials and growth conditions

Five P. mume cultivars (‘Songchun’, ‘Danfenghou’, ‘Zao Lve’, ‘Xiangxue Gongfen’ and ‘Fenhong Zhusha’) were obtained from the nursery of Beijing Forestry University. ‘Songchun’ and ‘Danfenghou’, which were developed by crossing Eumume and Apricot through one round of artificial hybridization, belonged to Apricot Mei group with strong cold resistance. ‘Zao Lve’, ‘Xiangxue Gongfen’ and ‘Fenhong Zhusha’ belonged to Eumume group with poor cold resistance. The Apricot rootstock seedlings of P. mume were cultivated in the greenhouse for 2 years. Plants with good growth status and consistent growth rates were selected in August 2019 and used in subsequent analysis (Supplementary Fig. S1).

Plants were grown for 1 week in a phytotron (16 h light/8 h dark cycle, 24 °C, 70% relative humidity and ~100 μmol·m^–2·s^–1 light intensity) prior to low temperature treatment. The cold hardening treatment was applied at 4 °C for 0, 1, 3, 5, 7, 9, and 11 d, and those at 24 °C were used as control. After 7 d, ‘Zao Lve’, belonging to Eumume group, was selected for freezing treatment test, while those at 24 °C were used as control. The selected plants were frozen at 0, –4, –8, –12 and –20 °C for 4 h using a horizontal chest freezer with a temperature controller. We randomly selected annual branches with the same growth potential and removed the bud points, which were used as materials for subsequent experiments.

Arabidopsis thaliana Col-0 plants were grown in an artificial climate chamber with an average temperature of 22 ± 2 °C, 16 h of light and 8 h of darkness, 100 μmol·m^–2·s^–1 light intensity, and 60–75% relative humidity.

Measurement of physiological indices

Physiological indices including soluble sugar, soluble protein, free proline, superoxide dismutase (SOD) and peroxidase (POD) activity were measured from five cultivars treated at 4 °C for different time lengths. The content of soluble sugar was measured colorimetrically by anthrone test using glucose as a standard (Moon and Sohn, 1988). The content of soluble protein was measured colorimetrically by Coomassie Brilliant Blue G-250 test using bovine serum albumin (BSA) as a standard. The amount of free proline was determined by sulfosalicylic acid spectrophotometry (Bates et al., 1973). POD and SOD activity were ascertained according to a previous method (Kaouthar et al., 2016). We used dimensionless function to integrate five physiological indices to evaluate the effect of cold hardening. Correlations between varieties were calculated using the Corrplot package with the ‘hclust’ method in the R project.

Determination of relative electrolyte leakage and recovery culture of plants

The relative electrolyte leakage (REL) of the freeze-treated plants was measured as described previously, using a conductivity meter (Li et al., 2017). We also performed freeze-thaw plant recovery culture. At first, the temperature was controlled to rise slowly to ~ 20 °C. Next, these plants were transferred to a phytotron (16 h light/8 h dark cycle, 24 °C, 70 % relative humidity and ~100 μmol·m^–2·s^–1 light intensity).

RNA extraction and library preparation

We prepared RNA-seq libraries from ‘Zao Lve’ branches at control (24 °C, 7 d, named as ‘CG’), cold hardening (4 °C, 7 d, named as ‘CE’), freezing (24 °C, 6 d + 20 h; –4 °C, 4 h, named as ‘NF’), and freezing after cold hardening (4 °C, 6 d + 20 h; –4 °C, 4 h, named as ‘CF’) treatments. Three individuals (obtained through asexual reproduction from ‘Zao Lve’) with consistent growth rate and phenotype were used for an independent test treatment, and three biological replicates were performed for each treatment. For sampling, after removing stem tip and leaf buds, a large number of ~1 cm long semi-lignified stem segments of 1-year-old branches were randomly collected from control and treatment groups. Total RNA of samples was extracted using Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). The degradation and contamination of total RNA were detected by 1% agarose gel electrophoresis. The purity and integrity of RNA were detected by the NanoPhotometer N80 (Implen, Germany) and Agilent 2100 Bioanalyzer (Agilent, USA), respectively. The libraries were sequenced using standard methods on an Illumina HiSeq2500 (Illumina, USA). All protocols were conducted as per the manufacturers’ instructions.

RNA sequencing and analysis

The splice sequence of raw reads, more than 10% unknown sequence (N) and 50% low-quality base reads were removed using Fastp software (Chen et al., 2018). To maximize the exclusion of rRNA interference, clean reads were compared with the ribosomal database using Bowtie2 software (Langmead and Salzberg, 2012), and reads compared with the ribosomal database were removed. The filtered reads were mapped to the P. mume reference genome (Zhang et al., 2012), and the distribution position of reads in the reference genome was counted using HISAT2 software (Kim et al., 2015). The transcripts were reconstructed, and the sequencing depth and transcript length were corrected using StringTie v1.3.1 software (Pertea et al., 2015). Differentially expressed genes (DEGs) were identified between samples using DESeq2 software with false discovery rate <0.05, |fold change| >2 and hypothesis test P-value <0.05 (Love et al., 2014). The DEGs were annotated and enriched using Genome-wide Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Ashburner et al., 2000; Kanehisa and Goto, 2000).

ATAC-seq library preparation

The materials used for ATAC-seq libraries were the same as that of RNA-seq libraries. Three biological replicates were prepared for each treatment. The nuclei were isolated from plant materials by the sucrose precipitation method (Bajic et al., 2018). Chromatin was fragmented and tagged following the standard ATAC-seq protocol (Buenrostro et al., 2015). PCR amplification and library construction of DNA were performed according to a previously tested method (Bajic et al., 2018). DNA fragments of ~100–800 bp were recovered by PCR amplification. The libraries were sequenced using an Illumina HiSeqTM 4000 platform (Buenrostro et al., 2013). All protocols were conducted as per the manufacturers’ instructions.

ATAC sequencing and analysis

The quality of the raw reads was tested using MultiQC software (Ewels et al., 2016). The flanking sequence of raw reads, more than 10% unknown sequence (N) and 50% low-quality base reads were removed using Trimmomatic software (Bolger et al., 2014). The reads were compared with the reference genome, and those compared to mitochondria or chloroplasts were removed using Bowtie2 software (Langmead and Salzberg, 2012). The distribution of reads relative to TSS position was analysed using deepTools software (Ramírez et al., 2016). Genome-wide open chromatin enrichment regions were retrieved using MACS v.2.1.2 software with parameters of --nomodel --shift -100 --extsize 200 -B -q 0.05 (Zhang et al., 2008). Peaks were filtered between repeated samples within the group, and the peaks with an overlap greater than 50% were retained when at least two samples co-existed. ChIPseeker software was used to confirm the genes related to transposase hypersensitive sites (THSs) and to analyse the distribution of THSs in different genomic regions (Yu et al., 2015). MEME suit (http://meme-suite.org/) was used to detect the motifs. MEME-ChIP was used to scan motifs with high reliability through peak regions, and MEME-AME was used to confirm the existence of any specific known motifs (Bailey et al., 2009; Machanick and Bailey, 2011).

In order to study the correlation between the level of chromatin opening and gene expression, the reads quantity of the gene promoter region was converted to reads per kilobase per million (RPKM) value to obtain the opening value of each gene in the sample. According to previous methods, chromatin openness could be divided into four levels, including high (RPKM>25), medium (5<RPKM≤25), low (1<RPKM≤5) and closed (RPKM≤1) (Guo et al., 2014; Meredith et al., 2015). Finally, the gene expression level and standard deviation of each sample in the four grades were calculated. Differential peaks were identified between samples using DiffBind software with |log₂(fold change)|≥1 and false discovery rate <0.05 (Ross-Innes et al., 2012). The peaks detected in the 2 kb upstream region of TSS were defined as the proximal peaks, and those detected in the upstream 2 kb to 100 kb region of TSS were defined as the distal peaks. The peak-associated genes were obtained on a genome-wide scale based on the regional information of the peak on the genome. All the differential proximal peak associated genes (DPGs) and distal peak associated genes (DDGs) were retrieved. The effect of the peak at the proximal or distal end on downstream gene expression were then analysed and compared.

Gene cloning and gene expression analysis

Plant genomic DNA and total RNA were extracted using DNAsecure Plant Kit (DP320, TIANGEN, Bejing, China) and RNAprep Pure Plant Plus Kit (DP441, TIANGEN), respectively. The cDNA library was constructed using PrimeScript™ II 1st Strand cDNA Synthesis Kit (6210A, TaKaRa, Bejing, China). The target genes were cloned by PrimeSTAR® HS (Premix) Kit (R040A, TaKaRa). Primers for gene cloning were designed using the Primer Premier 5.0 software (Supplementary Table S1). The reaction system and conditions of target gene cloning were maintained according to our previously published methods (Zhang et al., 2022; Zheng et al., 2022). qRT–PCR was performed with SYBR Premix EX Taq II (TaKaRa) on a CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The level of gene expression was detected by qRT–PCR using the reaction system and conditions described previously (Zhang et al., 2022; Zheng et al., 2022). The relative expression levels of genes were calculated using the 2-delta-delta cycle threshold (Ct) method, using the protein phosphatase 2A gene (PmPP2A) from P. mume as a reference (Zhang et al., 2012). Primers for gene expression are shown in Supplementary Table S1. The cis-acting elements of gene promoter sequences were analysed using PlantPAN 3.0 online tool (Chow et al., 2019). The codon preference in the coding region was analysed using NUCLEIC CODON USAGE included in the EMBOSS tool (Rice et al., 2000). The hydrophilicity of proteins was analysed using ProtParamand ProtScale tools (Wilkins et al., 1999).

Plant transformation and low temperature treatment

To verify the biological functions of candidate genes, the coding sequence (CDS) of PmCSL was obtained via PCR and inserted into the plant expression vector pCAMBIA1304 (stored in our laboratory), which contained a hygromycin resistance gene. The primers for overexpression vectors are shown in Supplementary Table S1. The recombinant plasmids were transformed into Agrobacterium tumefaciens (GV3101) by using the repetitive freeze thaw method (Weigel and Glazebrook, 2006). The DNA containing the target gene was integrated into the A. thaliana genome via the floral-dip method using Agrobacterium (Clough and Bent, 1998). To screen for transformants, the seeds were plated on solid Murashige and Skoog (MS) media supplemented with 50 mg l^–1 hygromycin. Viable and healthy Arabidopsis seedlings were transplanted into the soil and continued to be cultivated until seeds were obtained. Then, the positive Arabidopsis lines were screened for hygromycin resistance again until T₂ generation transgenic plants were obtained. At the same time, we used PCR detection technology for gene identification of transgenic plant lines. The genomic DNA of young leaves was extracted by DNAsecure Plant Kit (TIANGEN, Beijing, China). The hygromycin B 7'' -O-kinase (HYG) gene was used as a molecular marker gene, and the primers are shown in Supplementary Table S1. The reaction system and conditions of target gene cloning were maintained according to our previously published methods (Zhang et al., 2022; Zheng et al., 2022).

The disinfected wild-type and transgenic seeds were evenly planted on demand in solid medium (half-strength MS+ 3% sucrose+ 0.7% agar), vernalized (4 °C) for 2 d, and then transferred to a phytotron (16 h light/8 h dark cycle, 24 °C, 70% relative humidity and ~100 μmol·m^–2·s^–1 light intensity) for 14 d. The freezing treatment experiment was divided into two parts, direct freezing treatment (–5 °C, 1 h) and freezing after cold hardening treatment (4 °C, 3 d; –8 °C, 1 h). The survival rate and electrolyte leakage of the plants were subsequently measured. After freezing treatment, the seedlings were cultured at 4 °C for 1 d in the dark, and then transferred to a phytotron (16 h light/8 h dark cycle, 24 °C, 70% relative humidity and ~100 μmol·m^–2·s^–1 radiation) for 3 d (Li et al., 2017). We counted plant mortality and observed phenotypes.

Results

Physiological responses to cold hardening and freezing tolerance

To investigate the effects of cold hardening on the physiological and biochemical characteristics of P. mume, we selected five varieties with normal growth and flowering in the Beijing region, and measured the physiological and biochemical indices under cold hardening treatment (Fig. 1A). Cold hardening increased the content of soluble substances in P. mume. The soluble sugar content of the five varieties increased 2.1 times on average; for example, ‘Zao Lve’ increased from 6.73 mg g^–1 to 12.29 mg g^–1. The soluble protein content of five varieties showed an increasing trend during 0–5 d of cold hardening. The highest soluble protein content of 54.7 mg g^–1 and 52.89 mg g^–1 were observed in ‘Songchun’ and ‘Xiangxue Gongfen’, respectively. The soluble protein content of ‘Zao Lve’ was increased by 1.7-fold after 7 d. Meanwhile, the change in free proline content of ‘Zao Lve’ was the most prominent, reaching 365.11 μg g^–1 after 11 d of cold hardening, which was four times that of the control group (Fig. 1B). Cold hardening also affected POD and SOD activities of P. mume. We found that the optimal cold hardening effect was achieved after 9 d of cold treatment (Supplementary Table S2). For example, POD activity of ‘Danfenghou’ increased from 181.7 U g^–1 to 411.1 U g^–1, and that of ‘Zao Lve’ increased from 167 U g^–1 to 329.5 U g^–1. Collectively, the SOD activity first showed an increasing trend and then decreased gradually (Fig. 1B). In addition, ‘Zao Lve’ had the most prominent cold hardening effect among the Eumume branch varieties, and its freezing resistance was markedly enhanced (Fig. 1C; Supplementary Figs S2, S3).

Differentially expressed genes in response to low temperature

A total of 101.33 Gb of raw data were obtained from 12 samples using Illumina HiSeq 2500. The average clean dataset for each filtered sample was 8.36 Gb (Supplementary Table S3). According to principal component analysis (PCA) analysis, the dataset was divided into 12 dimensions, among which the variance contribution rate of PC1 dimension was 83.4%, and that of PC2 dimension was 8.6% (Supplementary Fig. S4A). The intra-group evaluation of PC1 and PC2 dimensions ranged from 320.9–1990.7, and 725.1–1559.2, respectively, and the Pearson correlation coefficients were all above 95% (Supplementary Fig. S4B).

The DEGs included 1281 up-regulated and 1227 down-regulated genes in the cold hardening treatment. A total of 809 genes were up-regulated, and 782 genes were down-regulated in the freezing treatment, and the number of DEGs was reduced compared with cold treatment. The freezing after cold hardening treatment group showed the highest number of DEGs, including 2202 up-regulated and 4182 down-regulated genes (Fig. 2A). Here, the number of down-regulated genes was approximately twice that of up-regulated genes. Among the 506 shared DEGs in three treatment groups, 223 were up-regulated and 283 were down-regulated (Fig. 2B; Supplementary Fig. S5A, B), suggesting that these genes might be broad spectrum low temperature-response genes.

The shared DEGs were divided into eight clusters (Supplementary Fig. S5C). In addition to their broad spectrum expression in different treatment groups, these DEGs were either induced or inhibited by different low temperature treatments. The top 30 DEGs with significant differences were selected for prediction of interactions. The results showed that the proteins encoded by these genes had strong interactions in each treatment condition (Fig. 2C-E; Supplementary Figs S6-S8). A total of 97, 57, and 122 transcription factors (TFs) were identified among DEGs in the cold hardening, freezing treatment, and freezing after cold hardening treatment groups, respectively (Supplementary Figs S9-S11). MYB transcription factors, Helix-Loop-Helix factors (bHLH), NAC transcription factors, ethylene-responsive element-binding factors (ERF), WRKY transcription factors, heat shock transcription factors (HSFs), and dehydration responsive element binding (DREB) protein family members were shared in the three treatment groups. Meanwhile, 91 DEGs were involved in calcium signalling, MAPK signalling, plant hormone signal transduction, and ICE-CBF signalling pathways (Fig. 2F; Supplementary Table S4). A total of 1365 DEGs had similar expression patterns in the cold hardening treatment and freezing after cold hardening treatment groups, of which 535 genes were up-regulated, and 830 genes were down-regulated. Of these DEGs, we found that 21 up-regulated and six down-regulated genes were enriched into the pathways containing 91 low temperature-response DEGs. Compared with the cold hardening treatment group, four CBFs, three basic region/leucine zipper transcription factors (bZIPs), five cold-responsive genes (CORs), calcineurin B-like protein 1 (CBL1), CBL-interacting serine/threonine-protein kinase 10 (CIPK10), mitogen-activated protein kinase kinase kinase 18 (MAPKKK18), mitogen-activated protein kinase kinase 9 (MKK9), protein phosphatase 2C 51 (PP2C51), abscisic acid receptor PYL2 (PYL2), and elongation of fatty acids protein 3-like (HOS3) genes were up-regulated in the freezing after cold hardening treatment group (Supplementary Table S4).

ATAC-seq analysis under low temperature treatment

We obtained a total of 212.3 Gb of raw reads and 106.4 Gb of clean reads from 12 samples through ATAC-seq. Base quality assessment showed that Q20 was above 90% and the average GC content was 40.2%, which was higher than the genomic GC content (38.33%; Supplementary Table S5). The average proportion of clean reads mapped to the nuclear, chloroplast, and mitochondrial genomes was 80.24%, 17.95% and 1.81%, respectively (Supplementary Fig. S12). After deduplication, the reads were unevenly distributed on chromosomes, with consistent distribution on positive and negative strands (Supplementary Fig. S13). Reads were mainly located upstream of the TSS of genes in all samples (Supplementary Figs S14-S19). More than 50% of the peaks were located in the 2 kb upstream of the TSS and the 300 bp downstream of the transcription termination site (TTS). The proportion of the promoter region of each sample was in the range of 46.72–51.14%, and the average proportion of the 300 bp downstream was 9.6% (Supplementary Fig. S15).

The chromatin region of the nucleosome was identified by THS, which potentially contained cis-acting elements that regulate the expression of neighbouring genes. In order to identify the genes related to the specificity of THS, each THS was matched to the TSS region of the closest gene. In CG, CE, NF, and CF groups, 6789, 5890, 12014 and 3187 THS were enriched, which were mapped to 5150, 4518, 9348, and 2359 neighbouring genes, respectively (Supplementary Table S6). These results indicated that a single gene was associated with one or more enriched THS, and most genes with enriched THS were associated with a single locus. Based on the comparison of THS mapping adjacent genes, 1617 shared THS-associated genes were retrieved in these groups (Fig. 3A). The quantitative data showed that the chromatin accessibility of the whole gene maintained a similar trend, but there were many different genes among the groups (Fig. 3B). Chromatin open regions potentially contained cis-acting elements that recruited transcription factors and regulated the expression of neighbouring genes. In CG, CE, NF and CF groups, 105, 108, 124, and 104 motifs were detected, corresponding to 104, 107, 123, and 103 TFs, respectively, of which 78 TFs were shared among the four groups (Supplementary Fig. S20). Based on the STRING database, we constructed an interaction network of 147 motifs corresponding to TFs, among which 124 TFs had interaction relationships (Fig. 3C). In the TF interaction network, bHLH family members were the most abundant, followed by bZIP and TCP family members.

Impact of chromatin accessibility on gene expression

The reads of the 2 kb upstream region of the TSS of genes were converted to RPKM values, and the open gene activity value of each gene was obtained. The average FPKM of RNA-seq was concentrated in the range of 5–15, and the average RPKM of ATAC-seq was concentrated in the range of 5–7.5 (Fig. 4A; Supplementary Fig. S21). In order to better understand the relationship between RNA-seq and ATAC-seq, the gene activity level was divided into four grades. A positive correlation was found between chromatin openness and gene expression level (Fig. 4B; Supplementary Fig. S22). The chromatin openness of 41 genes was highly open, with a wide range of gene activity and high expression levels (Supplementary Table S7). The functional annotation showed that most of the genes were related to chloroplast and mitochondria (Supplementary Fig. S22; Supplementary Table S7). We divided the open chromatin region into TSS proximal chromatin open region and distal chromatin open region. Correlation analysis of ATAC-seq and RNA-seq differential genes showed that the number of down-regulated associated genes was more than up-regulated genes, and the number of proximal TSS associated genes was greater than that of distal TSS (Supplementary Fig. S23). Here, 16 DPGs were detected in the CF down-regulated group (Supplementary Table S8). We found that the proximal chromatin openness of these genes was minimal, and their expression level was low in the CF group (Fig. 4C, D; Supplementary Table S9). These results further justified that chromatin openness was positively correlated with gene expression. Interestingly, three DPGs were detected in the CF up-regulated group, including cold-shock protein CS120-like (CSL), β-cyanoalanine synthase (CAS1), and calcium uniporter protein (MFO20.0; Fig. 4E, F). Multiple motifs that could bind TFs were found in the open chromatin region of the CSL gene. Meanwhile, the CSL gene was annotated to be closely related to plant response to low-temperature stress (Fig. 4F).

Chromatin accessibility, expression level and function analysis of PmCSL

PmCSL encodes a highly conserved protein with a dehydrin domain. The effective number of codons is 41.925, whose value is in the expected range of 20–61. The analysis results showed that the codon selection bias was moderate, and the contents of GC and GC3s were 52.38% and 54.66%, respectively, indicating that the coding region of PmCSL preferred G+C for base selection, and the codon preferred to end with G/C (Supplementary Fig. S24). PmCSL contained 78 phosphorylation sites phosphorylated by 14 specific phosphatases, and the main phosphorylation site was for Ca²⁺/calmodulin-dependent protein kinases (PKC; Supplementary Fig. S25). PmCSL was a hydrophilic protein, in which glycine (Gly) at positions 127 and 162 was the most hydrophilic amino acid (Supplementary Fig. S26). We identified 17 CSL sequences ranging from 154–331 amino acids, all of which belonged to Rosaceae in the NCBI database (Blastp, E-value<1e^–10). Phylogenetic results showed that CSL were highly conserved and closely related to each other in the same genus (Supplementary Fig. S27).

To further analyse the chromatin structure, the chromatin opening degree of PmCSL in response to low temperature was analysed. After cold or freezing treatment, significant enrichment of PmCSL peaks was observed on chromosome 8, especially in the upstream promoter region. The chromatin opening of TSS proximal –177 region of PmCSL was remarkably increased, and the gene expression level was notably up-regulated (Figs 4F, 5A, C). The bZIP binding site of PmCSL was found in the –215 to –236 region, and the ERF, WRKY, and TCP binding sites were also found in the chromatin open region (Fig. 5B). The expression level of PmCSL was up-regulated in response to cold hardening, and prolonged cold hardening promoted gene expression (Supplementary Fig. S28). The expression level of PmCSL further increased under freezing after cold hardening; nevertheless PmCSL was suddenly down-regulated under the –8 °C treatment of NF (Fig. 5D). Furthermore, PmCSL was ectopically overexpressed in Arabidopsis (Supplementary Fig. S29). The survival rate of PmCSL transgenic Arabidopsis plants was notably higher, and the electrolyte leakage was markedly lower, than that of wild-type plants (Fig. 5EH).

Discussion

As a major abiotic stress, low temperature not only affects plant growth and development, but also restricts the geographical distribution of plants (Zhu, 2016). For most plants, low temperature (0–15 ºC) does not cause fatal injury. For perennial plants thriving in cold and north temperate zones, cold acclimation is one of the most important factors to ensure winter protection. Cold acclimation/hardening has been widely used in freezing tolerance breeding of P. mume (Zhang, 1992; Chen et al., 2008). Soluble sugars and proteins are accumulated during cold acclimation/hardening, thus effectively preventing damage to plant cells caused by low temperature stress (Shahryar and Maali-Amiri, 2016). The freezing-resistance of P. mume ‘Zao Lve’ was significantly enhanced after cold hardening (Fig. 1C; Supplementary Figs S2, S3). As the cold hardening duration was prolonged, P. mume accumulated more soluble sugars and proteins. This upward trend was slightly decreased after 7 d (Fig. 1B). Our results are in line with previous studies documenting the increased freezing resistance of plants after cold hardening, such as A. thaliana, wheat, and peach (Warren et al., 1996; Hincha and Zuther, 2014; Muthuramalingam et al., 2022). Gene regulatory networks play pivotal roles at micro levels to instigate such responses.

The change of macromolecular substances in plants is a complex process, which involves changes in gene expression, signal transduction, and substance synthesis (Guo et al., 2018; Pu et al., 2019; Ambroise et al., 2020). Exposure to a non-freezing low temperature environment develops freezing in plants, which is regulated by the expression of genes related to the cold response (Örvar et al., 2000; Guo et al., 2018). In rapeseed, 690, 2538, and 2403 DEGs were identified by freezing treatment for 1, 3, and 24 h, respectively (Pu et al., 2019). In this study, the DEGs at 7 d of cold treatment accounted for 10.1% of all genes in the P. mume genome, and the DEGs at freezing after cold hardening treatment accounted for 25.7% of all genes in the genome (Fig. 2A). Thus, the freezing after cold hardening significantly increased the number of DEGs, including MYB, bHLH, NAC, ERF, WRKY, HSF, bZIP, and DREB transcription factor family members (Supplementary Figs S9-S11). Based on this, we hypothesize that compared with cold hardening and freezing treatment, freezing after cold hardening treatment is a more complex process in the cold response pathways. The induced expression of cold response genes not only requires a certain period of low temperature induction but it is also closely related to the temperature threshold (Fowler, 2008). Our study showed that cold hardening had a positive effect on the gene expression of P. mume in response to freezing stress.

The low temperature signal transduction pathways in plants include the calcium signalling pathway, MAPK pathway, ABA signal transduction, and ICE-CBF signalling pathways (Ma et al., 2015; Zhu, 2016; Guo et al., 2018; Shi et al., 2018). The spatial and temporal Ca²⁺ signal in plant cells, induced by the different stresses and decoded by many Ca²⁺ sensors, activates downstream cascades such as phosphorylation, leading to expression of stress-related genes (Yuan et al., 2018). In this study, three calcineurin B-like protein (CBL) and 13 CBL-interacting protein kinase (CIPK) genes were detected in the DEGs, including nine CIPK and three CBL genes in the freezing after cold hardening treatment (Supplementary Table S4). The CIPK-CBL complexes may connect low temperature-induced calcium signalling with the ICE or SnRK2.6/OST1 cascades, because CIPK-CBL binds to calcium and calmodulin (Zhu, 2016). The MAPK signalling pathway transmits Ca²⁺ signals to MEKK, MKK, and MPK kinases through calcium-activated protein kinase (CPK) transduction, and then activates the ICE-CBF signalling pathway (Shi et al., 2018). A total of eight MAPK-related kinases and one downstream gene of MAPK cascade (calmodulin-binding transcription activator, CAMTA) was detected in P. mume in response to low temperature treatment (Supplementary Table S4). We found that 27 DEGs were involved in the ABA signal transduction pathway. Serine/threonine-protein kinase (SnRK2.6/OST1) was up-regulated in response to cold signalling and further regulated ICE gene expression (Supplementary Table S4). However, some studies have shown that the process of SnRK2.6/OST1 regulating ICE does not belong to the ABA signal transduction pathway (Zhu, 2016). Forty-six genes related to the ICE-CBF-COR pathway were identified in this study (Supplementary Table S4). In conclusion, the Ca²⁺ signalling, MAPK cascade, ABA signalling, and ICE-CBF-COR signalling pathways also play an important role in P. mume in response to low temperature regulation.

In recent years, the development and improvement of ATAC-seq technology has greatly improved the detection of chromatin open regions in the whole genome of plants (Buenrostro et al., 2013; Lu et al., 2017). The open chromatin regions of A. thaliana, alfalfa, rice, and tomato exhibit similar distribution patterns in genomes (Maher et al., 2018). Meanwhile, the chromatin open regions show a similar distribution in the stem cells and mesophyll cells of A. thaliana and rice under drought and heat stresses (Wilkins et al., 2016; Sijacic et al., 2018). We found that more than 50% of the chromatin open regions in all samples were located 2 kb upstream and 300 bp downstream of the TTS. Chromatin open regions are mainly distributed upstream of the TSS, and the number of chromatin open regions in the intergenic region positively correlates with the size of the genome (Supplementary Figs S14-S19). The chromatin open region allows regulatory proteins such as transcription factors and co-factors to bind to cis-acting elements and regulate the expression of target genes. In rice, 445 TFs, including bZIP, ERF/AP2, and WRKY family members are involved in drought and heat stress regulation, and 115 TFs are involved in environmental gene regulatory networks (EGRINs; Wilkins et al., 2016). The AP2/ERF, WRKY, TCP, bZIP, and bHLH family members are enriched in the chromatin open area of tea tree under low temperature stress (Wang et al., 2021). Based on the ATAC-seq, we detected 107, 123, and 103 TFs in cold hardening, freezing, and freezing after cold hardening treatments, respectively, of P. mume (Supplementary Fig. S20). These results indicate that the chromatin structure of plant changes a lot under low temperature stress, and specific transcription factors are combined with cis-acting elements of loose DNA regulatory regions to regulate low temperature stress. Studies on humans and animals have shown that the integrated use of ATAC-seq and RNA-seq strategies has been quite helpful to address scientific problems and introduce innovative analytical methods (Ranzoni et al., 2021; Siahpirani et al., 2022).

Chromatin openness not only directly affects gene expression levels but also recruits transcription factors to regulate gene expression. In rice, under drought and heat stress, OsHSFA2a combines with cis-acting elements in the chromatin open region to regulate 53 target genes (Wilkins et al., 2016). Specific TFs bind to cis-acting elements in the open chromatin region to regulate the hypothermia inducing gene HD.04G0004270 under low temperature stress in tea tree (Wang et al., 2021). The chromatin opening of TSS proximal –177 region of PmCSL was significantly increased, and the gene expression level was significantly up-regulated (Fig. 5A). In this study, the bZIP, ERF, WRKY, and TCP binding sites were also found in the PmCSL promoter region, suggesting a positive correlation between chromatin accessibility and expression levels of target genes (Fig. 5B). The de-repression of genes at the chromatin level is thought to be the starting point of gene transcription (Park et al., 2018).

Cold shock proteins (CSPs) play an important role in the process of adaptation to low temperature environment and freezing tolerance. During cold acclimation, wheat WCSP is accumulated in the crown tissues of winter cultivars, but there is no significant difference in spring cultivars, inferring that WCSP plays a vital role in cold acclimation (Karlson and Imai, 2003; Radkova et al., 2014). Many studies have shown that DREB/CBFs transcription factors directly regulate the expression of COR genes, but CBFs usually have the highest expression level after 3–8 h of cold treatment (Cook et al., 2004; Hincha and Zuther, 2014; Guo et al., 2018; Shi et al., 2018). The expression level of PmCSL was up-regulated in response to cold hardening, and prolonged cold hardening promoted gene expression (Fig. 5C, D). Compared with freezing treatment, the expression level of PmCSL was significantly increased in freezing after cold hardening treatment. In A. thaliana, CSP gene members play different roles in the process of cold temperature regulation. AtCSP3 positively regulates plant freezing tolerance, while AtCSP2 negatively regulates cold acclimation (Sasaki et al., 2007; Kim et al., 2009). Overexpression of PmCSL in A. thaliana improves the freezing tolerance of transgenic plants, suggesting that PmCSL plays important roles in cold hardening and freezing tolerance of P. mume (Fig. 5E, H). Due to disadvantages of heterologous expression, the regulation mechanism and gene regulatory network of PmCSL in cold hardening still need further research.

Supplementary data

The following supplementary data are available at JXB online.

Fig. S1. Growth state of plant materials.

Fig. S2. Comprehensive evaluation of physiological and biochemical indices of different P. mume varieties.

Fig. S3. Freezing treatment and recovery culture of P. mume.

Fig. S4. Correlation analysis for RNA-seq data.

Fig. S5. Analysis of differential genes between treatment groups.

Fig. S6. The expression patterns of top 30 DEGs in P. mume following cold hardening treatment.

Fig. S7. The expression patterns of top 30 DEGs in P. mume following freezing treatment.

Fig. S8. The expression patterns of top 30 DEGs in P. mume following freezing after cold hardening treatment.

Fig. S9. Expression patterns of differentially expressed transcription factor genes following cold hardening.

Fig. S10. Expression patterns of differentially expressed transcription factor genes following freezing treatment.

Fig. S11. Expression patterns of differentially expressed transcription factor genes following freezing after cold hardening treatment.

Fig. S12. Reads distribution of chromatin open regions in the genome.

Fig. S13. The comparative chromosome analysis in the genome.

Fig. S14. Distribution of chromatin open regions.

Fig. S15. The distribution of reads relative to transcription start site locations of genes.

Fig. S16. The normalized reading signal of ATAC-seq samples in P. mume CG group enriched THS region.

Fig. S17. The normalized reading signal of ATAC-seq samples in P. mume CE group enriched THS region.

Fig. S18. The normalized reading signal of ATAC-seq samples in P. mume NF group enriched THS region.

Fig. S19. The normalized reading signal of ATAC-seq samples in P. mume CF group enriched THS region.

Fig. S20. Identification of motifs with open chromatin enrichment and comparison between groups.

Fig. S21. Scatter plot of the correlation between chromatin accessibility and gene expression.

Fig. S22. Correlation between chromatin accessibility and gene expression.

Fig. S23. Chromatin accessibility analysis of differentially expressed genes.

Fig. S24. Analysis of GC content and codon bias in the coding region of PmCSL.

Fig. S25. Phosphorylation site analysis of PmCSL amino acid sequence.

Fig. S26. Hydrophilic and hydrophobic analysis of PmCSL protein.

Fig. S27. The PmCSL phylogenetic trees retrieved by BLASTP with E-value <1e^–10.

Fig. S28. The relative expression of PmCSL at different times following cold hardening.

Fig. S29. PCR detection of transgenic plants.

Table S1. Primer sequences used in this study.

Table S2. Comprehensive evaluation of physiological and biochemical indices of P. mume during cold hardening.

Table S3. Statistical results of the RNA-seq dataset.

Table S4. Annotation of genes related to low temperature response and signal transduction pathway in P. mume.

Table S5. Statistical results of the ATAC-seq dataset.

Table S6. The distribution and associated genes of transposase hypersensitive sites.

Table S7. Functional annotation of high chromatin openness genes.

Table S8. The information of differential proximal associated genes in the CF down-regulated group.

Table S9. The DPGs between the CG and CF groups.

Author contributions

PL and TZ conceived and drafted the manuscript; PL performed the experiments; PL, LL, WL, SA, and LQ analysed the data; JW and TC provided plant resources. TZ and QZ contributed to the conception of the study and finalized the manuscript. All authors read and approved the final manuscript.

Conflict of interest

The authors have no conflict of interest to declare.

Funding

The research was supported by the National Key R and D Program of China (2019YFD1001500), the National Natural Science Foundation of China (32071816), and the Special Fund for Beijing Common Construction Project.

Data availability

The data that support the findings of this study, including the raw RNA-seq and ATAC-seq data are openly available in the NCBI-SRA under BioProject ID PRJNA851928.

References