Transcriptomic and physiological analysis reveals crucial biological pathways associated with low-temperature stress in Tunisian soft-seed pomegranate (Punica granatum L.)

ABSTRACT

Low temperature severely affects the growth of pomegranate in the early spring during the production process under protected cultivation. To understand the molecular responses to cold stress in Tunisian soft-seed pomegranate, this study investigated the transcriptome profiles and physiological changes of pomegranate leaves exposed to cold stress (6 °C) and freezing stress (0°C). Some potential cold response/resistance genes involved in plant hormone signal transduction, photosynthetic systems and carbon fixation in the C4 pathway, and sucrose and galactose metabolism were identified. In addition, an analysis of physiological indicators indicated that both stresses caused cell membrane damage; the accumulation of soluble sugar, soluble protein and proline; and the occurrence of photoinhibition owing to the damage in photosynthetic apparatus and the decrease in light energy conversion efficiency and electron transfer rate as shown by the decrease in net photosynthetic rate [Pn], potential maximum photochemical efficiency of PSII [Fv/Fm], actual photochemical efficiency of PSII [YII] and photochemical quenching coefficient [qP], and the effect was much moresevere in pomegranate under freezing stress. This study results offer useful information to understand the molecular mechanism of pomegranate response to cold stress and also lay a foundation for the selection of major candidate genes to conduct molecular breeding for cold tolerance in pomegranate.

KEYWORDS:

• Pomegranate
• cold stress
• transcriptome
• differentially expressed genes
• biological pathways

1. Introduction

Temperature is one of the key factors that affects the growth and development of crops. Low temperature not only affects the cellular metabolism of plants, resulting in slow growth and development, but also severely restricts their geographical distribution (Sharma and Pandey 2015). Plants have evolved complex cold tolerance mechanisms to avoid the negative effects of low temperature by managing low-temperature stress, including osmotic regulators, protective enzyme systems and biofilm systems (Thomashow 1999). Simultaneously, the subjection of plants to low-temperature stress can induce a range of signaling molecules, initiate or inhibit the expression of different genes to change the characteristics of cell membranes, regulate stomatal opening and closing, and regulate the synthesis and degradation of sugar and fatty acids, which induces the foundation of cold resistance mechanisms to improve the adaptability of plants to unfavorable temperatures (Chinnusamy et al. 2007).

Pomegranate (Punica granatum L.) is a fruit tree with extremely high economic value and ecological, ornamental and social benefits, which is native to Iran and widely planted in tropical and subtropical regions (Pourghayoumi et al. 2017). The primary pomegranate producing areas are India, China, Iran, the United States and Turkey (Silva et al. 2013). In recent years, Tunisian soft-seed pomegranate has become popular among consumers, and the largest planting areas are in China owing to its good taste, soft and edible kernel and high nutritional values. However, since Tunisian soft-seeded pomegranate is intolerant to cold, freezing injury can occur when the temperature is below −10°C and lasts for more than half a day in the winter (Feng et al. 2021). Therefore, in the last few years, the protected cultivation, including solar greenhouses and plastic-house cultivation, of Tunisian soft-seeded pomegranate has developed rapidly, which not only effectively avoids freezing injury, but also results in early market supply and significantly improves economic benefits. However, during the process of production, Tunisian soft-seed pomegranate is susceptible to low temperatures in early spring and extreme weather, such as cold spells in late spring, which results in the occurrence of cold injury of pomegranate buds and new shoots when the temperature is below 3°C and retardation in plant growth and development when the temperatures is below 10°C. Thus, cold weather causes at least 20–30% yield losses every year (Shen et al. 2021b). Therefore, it is of great significance to understand the mechanism by which pomegranate plants respond to low temperature and freezing stress and improve their cold resistance through gene transfer technology.

Plants need to undergo a complex signal transduction process from sensing low-temperature signals to making physiological and biochemical changes to defend or adapt to low-temperature stress, which is regulated by multiple genes and metabolic processes (Guo et al. 2018a). In recent years, RNA-Seq has successfully revealed the cold tolerance or adaptation mechanism of many plant species, such as Arabidopsis thaliana (Gehan et al. 2015), rice (Oryza sativa) (Dametto et al. 2015), kiwifruit (Actinidia deliciosa) (Sun et al. 2021), and pineapple (Ananas comosus) (Chen et al. 2016), because of the high accuracy and sensitivity of gene discovery (Sharma et al. 2018). Many cold stress response genes identified are widely involved in the metabolism of biological macromolecules, such as sugars, lipids, proteins and nucleic acids, and secondary production metabolism, energy metabolism and other metabolic pathways (Singh et al. 2020). Simultaneously, RNA-Seq has also identified many genes that respond to low temperature that are related to transcription factors (TFs) (Miura et al. 2007; Nakashima et al. 2009; Song et al. 2019), plant hormone signal transduction (Shi et al. 2012; Mehrotra et al. 2020), sucrose and starch (Li et al. 2021), and photosynthesis (Chen et al. 2018), which change significantly under low-temperature stress. However, to the best of our knowledge, studies on the regulatory mechanisms of pomegranate response to low-temperature stress have thus far primarily been focused on physiology and biochemistry (Moradi et al. 2017). There are still large gaps in understanding the response of pomegranate plants to low-temperature stress at the regulation of transcription. The purpose of this study was to explore the genomic characteristics of pomegranate adaption to cold. For this reason, this study sequenced and annotated the Tunisian soft-seeded pomegranate transcriptome under normal (control) and low-temperature conditions (6 and 0°C) using RNA-Seq and public databases. This study identified many key genes associated with responses to cold and analyzed vital crucial biological pathways that involve the key genes in pomegranate and also verified the physiological changes that some key genes may regulate. This study results can offer useful information to understand the molecular mechanism of pomegranate response to cold stress. It also lays a foundation for the protected cultivation of pomegranate and guides the future breeding of cold-resistant cultivars.

2. Materials and methods

2.1. Plant materials and experimental treatment

Tunisian soft-seed pomegranate seedlings with strong growth ability and no disease or mechanical damage were selected as experimental materials. They were obtained in Xingyang, Henan Province, China. This experiment was conducted in the experimental station of the Agricultural College of Shihezi University (Shihezi, China). The plants were cut uniformly and rooted in the sand bed of a greenhouse. They were then transferred to plant growth bags for unified cultivation and management. After 60 days, seedlings with the same growth trend were transferred to an artificial climate chamber (light cycle 16 h/8 h, day and night temperature 25°C/20°C, light 100 μmol·m⁻²·s⁻¹, relative humidity approximately 70%) for 30 days. The experiment was divided into three groups: the normal temperature control group (25°C, CK), cold stress treatment group (6°C) and freezing stress treatment group (0°C), with three replicates in each group. The seedlings of the control group were always cultured adaptively in the incubator. Except for the temperature, the other conditions of the two stress groups were the same as those of the control group. The stress group was reduced to the target temperature by continuously cooling for 3°C/h. After treatment for 48 h, samples were taken, and the leaves were cut, immediately frozen in liquid nitrogen and then stored at −80°C.

2.2. RNA preparation and library preparation for transcriptome sequencing

Sample RNA was extracted using a TianGen columnar plant RNA extraction kit (TianGen, Beijing, China). Total RNA was isolated using the TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA), and then the concentration, quality and integrity were determined by NanoDrop spectrophotometry (Thermo Fisher Scientific, Waltham, MA, USA). A TruSeq RNA sample preparation kit (Illumina, Ltd., San Diego, CA, USA) was used to generate the sequencing library. In addition, mRNA was purified using an AMPure XP system (Beckman Coulter, Brea, CA, USA). The product was purified using the AMPure XP system and quantified by Agilent high-sensitivity DNA analysis on a Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). The library was then sequenced by Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China) on an Illumina HiSeq 2500 platform.

2.3. Gene function annotation and expression analysis

FASTP was used to obtain clean readings by filtering low-quality reads. HISAT2 and default parameters were used to map all the clean data to the pomegranate reference genome. HTSeq statistics were used to compare the reads count value of each gene as the original expression of the genes. The reads count actively correlated with the authentic level of expression of the genes and positively correlated with the true level of expression of the gene, as well as its length and sequencing depth. The amount of expression was standardized by fragments per kilobase of transcript per million mapped reads (FPKM). The genes with FPKM > 1 were considered to be expressed. A larger FPKM value indicates a higher level of expression. DESeq (Wang et al. 2010) software was used to identify differentially expressed genes (DEGs) among the treatment group. The screening condition was that the threshold of false discovery rate (FDR) ≤ 0.05 and |log2 of fold change|>1. The DEGs were enriched and analyzed through Gene Ontology (GO) (Ashburner et al. 2000) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2004) database.

2.4. Validation of RNA-Seq data

The expression of 12 DEGs in Tunisian soft-seed pomegranate treated at control and low temperature conditions (6°C and 0°C) was analyzed by real-time quantitative reverse transcription PCR (qRT-PCR) to verify the RNA-Seq data (Trapnell et al. 2013). The extracted RNA samples were reverse transcribed using a HiFiScript cDNA synthesis kit (CoWin BioSciences, Cambridge, MA, USA), and the real-time fluorescence quantitative detection was conducted. Table 1 lists the primers used in qRT-PCR.

Table 1. The genes and primers used for qRT-PCR analysis.

Gene	Forward primer sequence	Reverse primer sequence
bglX (LOC116187211)	ATGCTTGAGAAGATTGATGC	GATCGTATAAGTTAGCACTG
bglB (LOC116214899)	GTGATTCTGTCAGAGAACGG	CGAAGTTATCAAGCAATGAC
TPS (LOC116210539)	ATCGGATGAGGACATGTTCG	TGAAGGACCTTGACAACTTC
PetH (LOC116214148)	TGCCAACTCTGATAGCCTAC	AAGATCTTGTCGCTGTACTCC
BAK1 (LOC116205131)	ATGATGAGGATGTCATGTTG	AGTAGAGCGACTTGGATCAG
BAK1 (LOC116205129)	TGGAGCAACTGATCCAAGTG	TCTGCCACTCCTTCCACCTC
GID1 (LOC116188586)	GAAGTACTTTGTCACGATCC	TTCAGGTAGAGCAGCTTCAC
ABF (LOC116205012)	AGGCATATACTGCAGAGCTG	GTCACAGGAACAACAATCTC
PYL (LOC116203110)	ACTTCTGGATGACAATGAGC	ACTTGAGGTTGCACCTGATC
ERF1 (LOC116197399)	GGTCGCTGGTGGCTGAGGAC	GTCCTCCAAGACCACCACAT
IAA (LOC116209597)	AGGAAGCTCATGGATCCAGC	ACTCGAACTCGTGCATGTAG
SAUR (LOC116208602)	TTGTTGCACTGGGATATCTC	AGGATCCTCTGCAGCTCCTG
Actin (LOC116211173)	TTCCCTGGTATTGCTGACAG	CGTACTCACCCTTGGAGATC

Note: actin: actin gene; bglX: beta-glucosidase isoform X1; bglB: beta-glucosidase 44-like isoform X1 gene; TPS: alpha-trehalose-phosphate synthase 9 gene; PetH: ferredoxin–NADP reductase; BAK1: brassinosteroid inensitive 1-associated receptor kinase 1-like isoform X1 gene and brassinosteroid insensitive 1-associated receptor kinase 1; GID1: gibberellin receptor GID1B-like gene; ABF: bZIP transcription factor 12-like gene; PYL: abscisic acid receptor PYL8 gene; ERF1: ethylene-response factor C3-like gene; IAA: auxin-responsive protein IAA17-like gene; SAUR: auxin-responsive protein SAUR71 gene.

2.5. Determination of photosynthetic parameters and physiological indices

The leaves of seedlings were used to determine the photosynthetic parameters and physiological indices. The relative electrolyte conductivity (EC) was measured with a conductivity meter. The contents of proline were determined by sulfosalicylic acid colorimetry (Zhang 1990). The contents of soluble protein and soluble sugar were determined using the Coomassie brilliant blue method (Chen and Wang 2015) and anthrone colorimetry (Chen and Wang 2015), respectively. Chlorophyll fluorescence was measured by determining the maximum photochemical efficiency (Fv/Fm), photochemical quenching coefficient (qP) and non-photochemical quenching coefficient (NPQ) of PSII in leaves of pomegranate seedlings using a pulse modulated chlorophyll fluorescence imaging system Max-Imaging-PAM. Each analysis consisted of three biological and technical replications.

2.6. Statistical analyzes

Excel 2010 software was used for data statistics, R language (mango, UK) software was used for data analysis and mapping, one-way analysis of variance (ANOVA) method was used to test the average value of each treatment with Duncan's multiple range method, the least significant difference (LSD) method was used to analyze the difference significance between different treatments (p<0.05). The numbers used in the tables and figures are means ± SE (N = 3).

Results

3.1. Statistics of transcriptome sequencing data

Transcriptome sequencing was performed on an Illumina HiSeq 2500 platform. Transcriptome sequencing of nine samples from three different temperature treatments was conducted. Each sample had three biological repeats, and 57.78Gb of clean data were obtained. Details of the reading sequence are shown in Table 2. The percentages of Q20 bases and Q30 bases of all samples > 97% and 93%, respectively. The clean reads of each sample were compared with the reference genome of pomegranate, and the alignment efficiency > 96%, indicating that the transcriptome data obtained was of good quality. The sequencing results conformed to the requirements of further analysis of pomegranate at the transcriptome level.

Table 2. Statistics of sequencing data.

Sample	Raw reads	Clean reads	Clean bases	Raw bases	Q20%	Q30%	Mapped ratio
CK-1	48,917,258	44,758,382	6.71G	7,337,588,700	97.68	93.82	96.69%
CK-2	53,218,738	48,581,648	7.29G	7,982,810,700	97.48	93.46	96.37%
CK-3	47,076,612	43,003,958	6.45G	7,061,491,800	97.36	93.17	96.25%
6°C-1	46,779,854	42,913,962	6.44G	7,016,978,100	97.51	93.57	96.73%
6°C-2	45,637,330	41,366,752	6.21G	6,845,599,500	97.43	93.38	96.22%
6°C-3	48,402,262	44,132,750	6.62G	7,260,339,300	97.69	93.87	96.33%
0°C-1	46,472,240	42,328,138	6.35G	6,970,836,000	97.56	93.62	96.31%
0°C-2	43,650,204	40,247,860	6.04G	6,547,530,600	97.46	93.41	96.58%
0°C-3	41,767,864	37,790,386	5.67G	6,265,179,600	97.35	93.25	96.01%
Average	46,880,262	42,791,537	6.42G	7,032,039,367	97.50	93.51	96.39%

3.2. Differential expression gene (DEG) analysis

Two comparison groups, including 6°C vs CK and 0°C vs CK, produced a large number of DEGs. As shown in Figure 1A, 7772 DEGs (3465 upregulated and 4307 downregulated) in 6°C vs CK and 3246 DEGs (1682 upregulated and 1564 downregulated) in 0°C vs CK were found, indicating that cold stress (6°C) induced numerous DEGs. In addition, as shown in Figure 1B, there were 1561 genes commonly expressed in both comparison groups, and these genes are involved in the response of plants to cold and freezing stress. We also identified 6211 and 1685 genes expressed only under cold stress and freezing stress treatments, respectively. It is apparent that 6211 genes are unique during the process of cold stress and only function under this treatment. They do not function during the freezing stress period. As shown in Figure 1C, cluster analysis was carried out according to the expression level of the same gene in different samples and the expression patterns of different genes in the same sample.

Figure 1. Histogram and Venn diagram of DEGs during cold stress and freezing stress. (A) Histogram shows the count of differentially expressed upregulated and downregulated genes under cold stress and freezing stress. (B) The numbers in the blue and orange circles represent the number of unique differential genes in the two treatment groups, and the overlap part represents the common differential genes between the two control groups. DEGs, differentially expressed genes. (C) Cluster analysis of differentially expressed genes under different treatments. Horizontal genes, each column is a sample, red indicates high expression genes, and blue indicates low expression genes.

3.3. Functional classification by GO

To determine the function of the DEGs identified, GO enrichment was performed using an FDR adjusted p-value of <0.05 as the cut-off. The DEGs were divided into functional groups that consisted of three categories, including Cell Component (CC), Molecular Function (MF), and Biological Processes (BP). The results clearly state that the DEGs were enriched into 219 and 192 GO terms in the comparison of 6°C vs CK and 0°C vs CK, respectively. The top 30 enriched GO terms are shown in Figure S1. Most of the DEGs from the 6°C vs CK were primarily enriched for photosynthesis, structural constituent of ribosome and cellular amide metabolic processes. Many DEGs from 0°C vs CK were involved in the membrane composition, tetrapyrrole binding and cell wall organization or biogenesis.

3.4. KEGG pathway analysis for the DEGs

This study explored the biological function of these DEGs in more detail by conducting a KEGG pathway enrichment analysis of DEGs using clusterProfiler. Compared with the control group, all the DEGs were allocated to 126 and 110 pathways under cold stress and freezing stress, respectively. Among these pathways, 19 and 13 were markedly enriched (p < 0.05), respectively. The major pathways involved in 6°C vs CK were ribosome, photosynthesis – antenna proteins and ribosome biogenesis in eukaryotes. In 0°C vs CK, the results demonstrated that the DEGs were primarily enriched in plant-pathogen interaction, phenylpropanoid biosynthesis, flavonoid biosynthesis and starch and sucrose metabolism. The three most significantly enriched pathways were photosynthesis, starch and sucrose metabolism and protein processing in endoplasmic reticulum for both stress treatments (Figure 2). This suggests that these pathways could play a positive role in the response to cold stress and freezing stress.

3.5. Analysis of the DEGs involved the cold and freezing stress responses

3.5.1 DEGs related to transcription factors

In the comparison groups of 6°C vs CK and 0°C vs CK, 1938 TFs in 55 TF families and 762 TFs in 52 TF families were found to be involved in the responses to cold and freezing stress, respectively (Figure 3). Different members of major TF families were upregulated or downregulated under cold and freezing stress, such as bHLH, NAC, ERF, WRKY, C2H2, MYB, MYB_related and bZIP. Based on the number of genes in each TF family identified, the bHLH and ERF families were the largest classes in response to cold and freezing stress, with 168 TFs (70 up- and 98 downregulated, 8.6%) and 85 TFs (57 up- and 28 downregulated, 11.1%), respectively, followed by the NAC family (81 up- and 66 downregulated, 7.5%) and the ERF family (75 up- and 61 downregulated, 7%) under cold stress (Figure 3A), the bHLH family (26 up- and 32 downregulated, 7.6%) and the NAC family (32 up- and 26 downregulated, 7.6%) under freezing stress (Figure 3B).

Figure 2. DEG pathway allocation base on the basis of KEGG under cold stress (A) and freezing stress (B), respectively. DEG, differentially expressed gene; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 3. Analysis of TFs that are associated with the response to cold stress and freezing stress. The number of upregulated and downregulated genes in each bar is shown in orange and blue, respectively. TFs, transcription factors.

3.5.2. DEGs related to the osmoregulation system

The KEGG pathway enrichment analysis indicated that the ‘starch and sucrose metabolism’ (Ko00500) and ‘galactose metabolism’ (Ko00052) pathways enriched many DGEs that regulated and relieved the osmotic stress caused by low temperature. A total of 72 and 38 DEGs were identified in two metabolic pathways, respectively. In Starch and sucrose metabolism pathway, two α-trehalose-phosphate synthase genes (LOC116203594 and LOC116210539), one trehalose-phosphate phosphatase A gene (LOC116206238), one β-glucosidase 40-like gene (LOC116212511) and one hexokinase-3-like gene (LOC116211459) were upregulated in response to both stresses (Figure 4, Table S1). In the galactose metabolism pathway, four DEGs were upregulated under both cold and freezing conditions, including one sucrose galactosyltransferase1 gene (LOC116192468), one hexokinase-3-like gene (LOC116211459) and two galactinol synthase 2-like genes (LOC116215737 and LOC116212865) (Figure 4, Table S1).

Figure 4. Heat map of the DEGs involved in the osmoregulation system. The heat map was generated from the log2FC. Changes in the level of expression are represented by a change in color. Blue indicates a lower level of expression, whereas red indicates a higher level of expression. DEGs, differentially expressed genes.

3.5.3. Key DEGs related to the hormone signal transduction pathway

The perception and transmission of cold signals are important for plants during low temperature stress. There were 98 and 52 DEGs involved in the plant hormone signal transduction pathway (ko04075) under cold and freezing stress, respectively. The primary DEG pathways included auxin (IAA), brassinosteroid (BR), abscisic acid (ABA), ethylene (ETH) and gibberellin (GA) pathways, indicating that multiple processes in plant hormone signal transduction pathways in the leaves of pomegranate are activated under cold and freezing stress (Figure 5). This study found that the genes in relation to IAA signal transduction pathway were the most numerous with 31 (7 up- and 24 downregulated) and 33 (12 up- and 11 downregulated) DEGs identified under cold and freezing stress, respectively. Among them, five DEGs were upregulated under both cold and freezing conditions, including IAA (LOC116209597, LOC116188243, and LOC116188658) and SAUR (LOC116208602 and LOC116207863). In addition, this study identified 23 and 4 DEGs related to ABA biosynthesis under cold and freezing stress, respectively. Among the 23 DEGs under cold stress, four PP2C genes (LOC116196060, LOC116198959, LOC116204689 and LOC116203216) were significantly upregulated by 7.93-, 5.83-, 4.6- and 4.36-fold, respectively. Under freezing stress, all four DEGs, including ABF (LOC116202490), PYL(LOC116202208 and LOC116198629) and PP2C (LOC116203216) were upregulated. In the BR signal transduction pathway, 19 and 10 DEGs under cold and freeing stress also were found, respectively. The DEGs were primarily involved in BR associated receptor kinase (BAK1), BR-signaling kinase (BSK) and cyclin D3 (CYCD3), and most them were down-regulated. There are only two BSK (LOC116193567, LOC116192908) under cold stress and one BAK1 (LOC116196446) and one BSK (LOC116213130) under freezing stress were up-regulated (Figure 5, Table S2).

Figure 5. Schematic diagram of DEGs involved in hormone-related pathways. The heat map was generated from the log2FC. Changes in the level of expression are represented by a change in color. Blue indicates a lower level of expression level, whereas red indicates a higher level of expression.

3.5.4. DEGs are vitally interrelated in photosynthesis, photosynthesis-antenna proteins and carbon fixation in photosynthetic organisms

Photosynthesis is sensitive to low-temperature stress. The results of the KEGG pathway analysis showed that many DEGs were involved in photosynthesis. This study identified 14 and 4 DEGs related to light-harvesting complex I and II chlorophyll a/b binding protein (LHCA and LHCB) under cold and freezing conditions, respectively. The expression of all 14 DEGs were downregulated under cold stress, particularly LHCA2 (LOC116188492), LHCA4 (LOC116189511), LHCA5 (LOC116215606) and LHCB4 (LOC116209479) for which expression were significantly downregulated by 5.26-, 3.98-, 4.30-and 4.09-fold, respectively. However, the levels of expression of all four DEGs, including LHCA2 (LOC116212466), LHCB1 (LOC116195644), LHCB2 (LOC116203120) and LHCB5 (LOC116204457), were significantly upregulated and higher under freezing stress than cold stress (Figure 6, Table S3).

Figure 6. Effects of different temperature treatments on the expression of DEGs involved in the photosynthetic pathway. A heat map of DEGs involved in the photosynthetic system, including photosynthesis (A), photosynthesis – antenna proteins (B) and carbon fixation in photosynthetic organisms (C). DEGs, differentially expressed genes.

A total of 39 and 17 DEGs associated with photosystem I, photosystem II complex, cytochrome b6-f complex, ATP synthase, and oxygen-evolving complex were also identified under cold and freezing conditions, respectively. Among the 39 DEGs, the levels of expression of psbD (CGW82_pgp069), petF (LOC116201807 and LOC116210390) were significantly upregulated, and 36 DEGs were downregulated under cold stress. All 17 genes were upregulated under freezing stress, particularly psbD (CGW82_pgp069), ATPF0B (GW82_pgp078) and petF (LOC116210390) for which expression was significantly upregulated by 5.36-, 5.25- and 4.58-fold, respectively (Figure 6, Table S3).

In addition, 37 and 3 DGEs that mapped to carbon fixation in the photosynthetic organisms pathway were found under cold and freezing stress, respectively. Under cold stress, 14 of the 26 downregulated DEGs were related to the Calvin cycle, primarily including rbcS (LOC116204098), SBPase (LOC116195124), FBP (LOC116193140 and LOC116203977), PRK (LOC116188159 and LOC116204435), PGK (LOC116205969) and rpiA (LOC116214078 and LOC116208224). Six of 11 upregulated DEGs were related to the C4 pathway, primarily including PPC (LOC116202419 and LOC116200004), GOT1 (LOC116187710 and LOC116192188), maeB (LOC116201315) and MDH2 (LOC116209134). However, under freezing stress, two rbcL genes (LOC116190969 and CGW82_pgp056) were upregulated by 5.21- and 4.39-fold, respectively, and only one PPC gene (LOC116213702) was downregulated (Figure 6, Table S3).

3.6. Verification of DEGs using qRT-PCR

The transcript abundance of 12 selected DEGs was verified by qRT-PCR to check the authenticity of RNA-Seq data. All 12 regulatory genes were obviously induced by cold conditions, including genes involved in the osmotic regulation system, photosynthetic system and various hormone metabolic pathways. The results show that there is a strong correlation with the RNA-Seq data, thus, proving the reliability of data (Figure 7).

Figure 7. A correlation of the levels of expression of 12 DEGs verified between RNA-Seq and qRT-PCR. The Y-axis on the left shows the relative levels of gene expression (2^−ΔΔCT) analyzed by qRT-PCR (orange pillars), while the Y-axis on the right shows the corresponding RNA-Seq expression data (bule pillars). The X-axis shows the different temperature treatments. DEGs, differentially expressed genes; qRT-PCR, real-time quantitative reverse transcription PCR.

3.7. Net photosynthetic rate and chlorophyll fluorescence parameters

Compared with the control, the net photosynthetic rate (Pn), potential maximum photochemical efficiency of PSII (Fv/Fm), photochemical quenching coefficient (qP) and actual maximum photochemical efficiency of PS II (Y II) decreased significantly, and the non-photochemical quenching coefficient (NPQ) increased in the leaves of pomegranate seedlings under both stresses; however, the values of Pn, Fv/Fm and qP were higher under cold stress than those under freezing stress (Figure 8).

Figure 8. Determination of net photosynthetic rate and chlorophyll fluorescence parameters. Subsequent growth and development of soft-seed pomegranate seedlings after stress treatments image (A), the image of Fv/Fm (B), Fv/Fm (C), Pn (D), YII(E), qP F), and NPQ (G) numerical values. Data represent the means and standard errors of three replicates. The chart is marked with different lowercase letters, indicating that the difference is significant at P < 0.05. Fv/FM, potential maximum photochemical efficiency of PS II; NPQ, non-photochemical quenching coefficient; qP, photochemical quenching coefficient; Pn, net photosynthetic rate; Y II, actual maximum photochemical efficiency of PS II.

As shown in Figure 8A, the Fv/Fm fluorescence imaging results under cold and freezing stress were consistent with the values of Fv/Fm in Figure 8B. Compared with the control, the area of blue decreased, and the area of yellow-green increased in leaves under both stresses, particularly under freezing stress.

3.8. Relative electrolyte conductivity and the content of osmoregulation substances under cold and freezing stress

Compared with the control, the relative electrolyte conductivity (EC), the contents of proline, soluble protein and soluble sugar significantly increased in the leaves of pomegranate seedlings under both stresses, and the values of these indices were lower under cold stress than under freezing stress (Figure 9).

Figure 9. Analysis of relative electrolyte conductivity (A), proline (B), soluble protein (C) and soluble sugar (D) contents under low temperatures.

4. Discussion

Among various abiotic stresses, cold stress seriously affects crop growth and development and yield (Knight and Knight 2012). Pomegranate is an economically important fruit tree and widely planted in tropical and subtropical regions, which is susceptible to low temperatures in early spring and extreme weather, such as cold spells in the late spring (Jayesh and Kumar 2004). Understanding the molecular mechanism of pomegranate in response to cold stress would aid in gene function studies and the breeding of cold-tolerant varieties. This study reported the transcriptome information of Tunisian soft-seed pomegranate under cold and freezing stress by RNA-Seq. A total of 10,988 DEGs were obtained with 7772 and 3246 DEGs responding to cold and freezing stresses, respectively. Among them, 1561 DEGs were commonly regulated by both stresses, and 6211 DEGs were specifically expressed under cold stress, indicating that a larger number of genes was induced under cold stress compared with that under freezing stress, and there were many common and unique molecular response mechanisms to cold stress in pomegranate (Figure 1). GO and KEGG enrichment analyzes also indicated that most DEGs were classified in response to cold stress and involved in multiple metabolic pathways, including TFs, photosynthesis-antenna proteins and carbon fixation in photosynthetic organisms, starch and sucrose metabolism, and plant hormone signal transduction.

4.1. Abundant TFs that were activated under cold stress in pomegranate

When plants are exposed to cold stress, TFs can regulate the target genes related to cold resistance by binding specific cis regulatory elements in their promoters, thereby enhancing cold tolerance (Thomashow 1999; Ke et al. 2020). In recent years, TF families, including AP2/ERF, MYB, bHLH, NAC and WRKY, have attracted much attention as important regulators involved in plant stress responses (Bhardwaj et al. 2015; Palmeros-Suárez et al. 2015). This study found that the bHLH, NAC and ERF families were the three primary TF groups that responded to cold and freezing stress. The bHLH and NAC families are considered to be play important roles in response to abiotic stresses, including low-temperature stress. ERF is a TF of the AP2/ERF subfamily and is also involved in cold regulation (Wu et al. 2015b; Du et al. 2016). Previous studies identified many bHLHs, NACs and ERFs involved in the defenses to low temperature stress (Guo et al. 2021). The heterologous expression of the two bHLH genes VaICE1 and VaICE2 from Amur grape (Vitis amurensis Rupr.) in transgenic A. thaliana enhances the cold tolerance of plants (Xu et al. 2014). A total of 66 bHLH genes could be involved in the regulation of responses of sweet cherry (Prunus avium L.) to cold stress (Shen et al. 2021a). The overexpression of TERF2 confers cold tolerance in rice (Tian et al. 2011). The overexpression of DREB1 from several plant species in transgenic A. thaliana also conferred tolerance to cold stress (Yadav 2010). Most bHLHs and ERFs under both stresses and NACs under cold stress were upregulated in this study, suggesting that bHLHs, ERFs and NACs could play an important positive regulatory role in the response of pomegranate to cold stress (Dong et al. 2019). Moreover, in this study, 35 and 18 MYBs, and 39 and 28 WRKYs, were upregulated under cold and freezing stress, respectively. The MYB and WRKY families are also involved in the tolerance of plants to low temperature stress (Katiyar et al. 2012; Yang et al. 2020). These results suggest that multiple genes and TFs are involved in the regulatory mechanism of pomegranate in response to cold and freezing stress through different manners.

4.2. Genes associated with soluble sugar metabolic pathways are crucial for plants against cold stress

When plants are subjected to abiotic stress, the cells will rapidly accumulate osmotic regulatory substances, such as soluble sugar, soluble protein and proline, to maintain or reduce the osmotic potential of cells (Peng et al. 2012; Abeynayake et al. 2015). In addition, it has been reported that these substances also play important roles in maintaining the stability of cell membranes by activating reactive oxygen species (ROS) scavenging systems (Wu et al. 2015a). In this study, both stresses caused damage to the membrane systems and induced the accumulation of soluble sugar, soluble protein and proline in pomegranate leaves, particularly under freezing stress (Figure 4), suggesting that the accumulation of these osmoregulatory substances helped to alleviate the damage from cell dehydration caused by cold stress and freezing stress, but the cell membrane of pomegranate suffered more severe damage under freezing stress than under cold stress.

To manage the damage from dehydration caused by low-temperature stress, plants can induce the accumulation of soluble sugars to regulate their osmotic potential and maintain cell turgor (Cheong et al. 2019). In this study, we found many DEGs that participated in the accumulation of soluble sugars. Among them, α/β-amylase (α/βAMY), β-glucosidase (bglU, bglX and bglB), fructofuranosidase (sacA), and hexokinase (HK) were involved in the procession of degrading starch and cellulose, and the expression of the genes that encoded these enzymes were upregulated to differing degrees under cold stress. The results are similar to those of Xu et al. (Wu et al. 2015a) in Camellia weiningensis and Wu et al. (Fernandez et al. 2010) in tea oil Camellia (C. oleifera). Moreover, the genes related to the synthesis of trehalose and raffinose were also upregulated under cold stress, including trehalose phosphate synthase (TPS), trehalose phosphatase (TP), UDP-glucose 4-epimerase (galE), UTP–glucose-1-phosphate uridylyltransferase (UGP2), galactosyltransferase (GOLS) and galactinol synthase (GS). Trehalose is an important multifunctional disaccharide in response to abiotic stress (Fernandez et al. 2010). When plants are subjected to abiotic stress, trehalose forms a special protective structure on the surface of cell membrane to protect the activity of protein molecules (Mostofa et al. 2015). It has been shown that the metabolic pathway in which GolS is involved plays an important role in the storage and transportation of plant carbon assimilates, biological and abiotic stress responses and other life processes (Cong et al. 2019). The overexpression of GolS enhanced the tolerance to Zinc stress by increasing the levels of raffinose and ROS-scavenging capacity in rice and A. thaliana (Wang et al. 2016). In addition, this study found that UGP2, which catalyzes the formation of uridine diphosphate glucose (UDPG), was upregulated under cold stress and downregulated under freezing stress. UDPG participates in the synthesis of sucrose, and it helps to maintain the stability of soluble sugar content under cold stress. Therefore, cold stress induced pomegranate to establish a better osmotic adjustment system to resist cold damage. The metabolic mechanisms of various sugars in pomegranate under cold stress warrant further study, particularly in combination with metabolomics.

4.3. Plant hormone-dependent signaling pathways play vital roles in response to cold stress

The signaling pathways of plant hormones play vital roles in recognizing and transmitting cold signals, which then trigger a series of cold-induced responses (Guo et al. 2018b). In this study, 98 and 52 DEGs were involved in ABA, IAA, BR and GA pathways under cold and freezing stress, respectively. Among them, IAA, ABA ang BR are the primary responsive hormone signal transduction pathways under both stresses, and large numbers of DEGs are involved in them.

Auxin has been proven to be one type of phytohormone involved in plant responses to abiotic stress (Jain and Khurana 2009). In non-canonical auxin signal transduction pathways, auxin can induce the formation of a TIR1/AFB (transport inhibitor response 1/auxin signaling F-box)-Aux/IAA (auxin/indole-3-acetic acid) complex and the subsequent ubiquitination and degradation of Aux/IAAs, enabling the released ARFs to regulate the transcription of auxin-responsive genes, such as Aux/IAA, Gretchen Hagen3 (GH3) and small auxin-up RNA (SAUR). GH3 and SAURs are primarily responsible for regulating auxin homeostasis and responding to changes in auxin. They play an important role in the auxin signaling pathway, light signaling pathway and plant stress responses (Yu et al. 2022). In this study, two TIR1 and two ARF were inhibited under cold stress. One ARF was upregulated under cold stress but downregulated under freezing stress. Eight down- and three upregulated IAA genes were found under cold and freezing stress, respectively. One IAA, one AUX1 and five SAURs were downregulated, and two IAA and CH3 genes were upregulated under both stress conditions. Auxin/IAA acts as the hub of integrating genetic and environmental information to achieve plant stress resistance (Tiwari et al. 2004). The ARF family is an important family of TFs that are involved in the transcriptional regulation of auxin-induced genes. Recent studies have shown that repression of the auxin signaling pathway and the downregulation in the expression of genes positively responsive to the auxin signaling pathway, such as Aux ⁄IAA and ARF gene family members, are an important aspect of the defense response to abiotic stresses (Kudo et al. 2012; Magwanga et al. 2019). Therefore, the TIR1/AFB-Aux/IAA-ARF signal transduction pathway could play an important role in response to cold and freezing stresses in pomegranate.

Studies have shown that the accumulation of ABA in plants under low-temperature stress can transmit environmental stimuli to induce ABA-dependent signaling pathways (de Zelicourt et al. 2016). ABA can activate serine/threonine-protein kinase (SNRK2) via the binding of ABA receptor complex (ABA-PYR/PYL/RCAR) and protein phosphatases 2C(PP2C), and then regulate the expression of downstream genes, such as ABA binding factors/ABA responsive element-binding protein (ABF/AREB) and ABA insensitive 5 protein (ABI5) among others, thus, directly or indirectly regulating plant cold tolerance (Figure 7b) (Roskoski 2010; Yan et al. 2020). This study identified 18 DEGs (11 up- and 7 downregulated) in the ABA signaling pathway under cold stress. Among 11 upregulated genes, including two ABA receptor genes (PYL), five PP2C, one SnRK2, and four ABF. However, there were only three upregulated genes, including one PYL, one PP2C and one SnRK2, that were identified under freezing stress. Kamiyama et al. (2021) and Nian et al. (2021) reported that the enhanced levels of expression of SNRK2 and PYL can improve the cold tolerance of plants. Many studies have found that PP2C was involved in the regulation of responses to cold stress (Zhang et al. 2021). Therefore, this study results also proved that the ABA signal transduction pathway played an important role in the response of pomegranate to cold stress. However, in this study, the levels of expression of two PYL and one SNRK2 gene were downregulated under cold stress, and one SNRK2 gene (LOC116213130) was downregulated under cold stress but upregulated under freezing stress, indicating the complex and diverse roles of PYL and SNRK2 in ABA signal transduction under cold stress. However, the specific regulatory mechanisms merit further exploration.

BR also has a positive effect on resistance to low-temperature stress. A previous study reported that the application of BR in cucumber can increase its cold tolerance (Jiang et al. 2013). In the signal pathway of BR, BSK can trigger downstream signal transduction cascade via phosphorylated by BRI1. By contrast, BIN2 phosphorylates and inactivates BR response factors, BZR1 (BR residtant 1) and BZR2, to inhibite transcription of BR-responsive genes (Figure 7b). In this study, two up-regulated BSK (LOC116193567, LOC116192908) may function and play key role in pomegranate cold response. Therefore, the cold response in pomegranate is a complicated regulatory network, involving a variety of plant hormone signaling pathway.

4.4 Photoinhibition occurred under cold and freezing stress, but many genes in photosynthesis and the carbon fixation pathway may be crucial for plants to acclimate cold stress

Low temperature could induce photoinhibition, decrease the absorption and ability to capture light energy, electron transfer efficiency of the photosystem, and fixed quantum efficiency of carbon dioxide (Yang et al. 2019). In this study, Pn, Fv/Fm and qP in the leaves of pomegranate seedlings significantly decreased under both stresses, particularly under freezing stress, indicating that both stresses caused damage to the photosynthetic apparatus and inhibited the electron transfer rate and light energy conversion efficiency. These effects led to the occurrence of photoinhibition, and the effect was much more severe in pomegranate under freezing stress. In addition, the increase in NPQ under both stresses indicated that the self-protective mechanism was activated to reduce excess light energy by increasing the dissipation of heat (Figure 8).

Light energy is captured and converted into chemical energy in the first step of photosynthesis. In this study, under cold stress, all 14 DEGs mapped to the photosynthesis-antenna proteins pathway related to LHCA and LHCB that were downregulated under cold stress. This result is similar with those of Li et al. (2019) in maize. Wu et al. (2020) also reported that the expression of LHCA and LHCB significantly decreased under cold stress. Thus, downregulation of these LHCA and LHCB genes in this study suggests that cold stress decreased the ability to absorb and capture light energy. It is worth noting that all four genes LHCA2 (LOC116212466), LHCB1(LOC116195644), LHCB2(LOC116203120) and LHCB5(LOC116204457) were upregulated under freezing stress, which could be an adaptive protective response, and the specific reasons merit further study.

Under conditions of adversity, the plants must maintain as high an electron transfer rate and light energy conversion efficiency as possible to provide more NADP and ATP for carbon assimilation. In this study, 36 of 39 DEGs were related to the structure of photosystem complex, cytochrome b6/f complex, ATP synthase, and oxygen-evolving complex, and their levels of expression were significantly downregulated under cold stress. The downregulated DEGs primarily included ATPFOB and ATPF1D, which encode the ATPase subunit; PsbB, PsbO and PsbP, which encode the PSII oxygen-evolving enhancer protein gene; PsbW, which encodes the PSII reaction center W protein; Psb28, which encodes the photosystem II reaction center PSB28 protein and PsbQ, which encodes the subunit of PSII complex in PSII; pasB, PsaD, PsaE, pasF, PsaK, PsaL, PsaN and PsaO,

which encode the subunit of PSI complex in PSI; PetC which encodes the cytochrome b6/f iron-sulfur subunit; and petG, which encodes cytochrome b6. The changes in regulation of these genes indicated that cold stress could damage the structure of the photosystem and decrease the electron transport rate. This result is similar to those of Xu and Yu (2022) in C. weiningensis and Xu et al. (2020) in transgenic banana (Musa spp.). Nevertheless, one psbD that encodes ferredoxin-NADP⁺ reductase, CGW82_pgp069) and two petF that encode ferredoxin, LOC116201807 and LOC116210390) were upregulated under cold stress, indicating that they may have positive transcription regulation in response to cold stress. However, all 17 DEGs were upregulated under freezing stress, but the values of Pn, Fv/Fm and qP decreased (Figure 8). Thus, this study hypothesized that this could also be an adaptive response, and the enhanced expression of photosynthesis-related genes is insufficient to maintain a certain level of photosynthetic capacity under freezing stress. Further study is merited to clarify their roles under different low-temperature conditions.

In this study, it was found that 14 of the 26 downregulated DEGs that mapped to carbon fixation in the photosynthetic organisms pathway were related to the Calvin cycle, which encodes the ribulose bisphosphate carboxylase small chain (rbcS), sedoheptulose-1,7-bisphosphatase (SBPase), fructose-1,6-bisphosphatase (FBP), phosphoribulokinase (PRK), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), phosphoglycerate kinase (PGK), and ribose-5-phosphate isomerase (RPI). Many studies have shown that the decrease in plant carbon assimilation efficiency and rate of regeneration of ribulose 1,5-biphosphate (RuBP) at low temperature could be owing to the downregulation of key enzyme activities and levels of expression of gene transcription in the Calvin cycle (Peng et al. 2012; Abeynayake et al. 2015). Therefore, this study results suggest that cold stress decreased the efficiency of carbon assimilation. However, this study also found that six of the 11 upregulated DEGs under cold stress were related to the C4 pathway. Among them, two PPC (LOC116202419 and LOC116200004), two GOT1 (LOC116187710 and LOC116192188), one maeB (LOC116201315) and one MDH2 (LOC116209134) were induced by cold stress. Studies have shown that the activation of C4 pathway in C3 plants under stress conditions is an adaptive adjustment to the environment and survival optimization strategy, and the functions of PPC and maeB appear to be more important for plants under stress than under optimal growth conditions (Doubnerová and Ryslavá 2011). Therefore, these genes were believed to play an important role to maintain and enhance the efficiency of carbon assimilation in pomegranate under cold stress. Further research is required to clarify their functions in different cold conditions.

5. Conclusions

In this study, it reported the transcriptomic and physiological responses of Tunisian soft-seed pomegranate to cold and freezing stress. DEGs related to TFs, photosynthesis, osmotic regulation system, hormone signal transduction were found to be involved in the response to cold stress. Some potential cold response/resistance genes involved in plant hormone signal transduction, photosynthetic systems and carbon fixation in the C4 pathway, and sucrose and galactose metabolism were identified. In addition, analyzes of physiological indicators also proved that both stresses caused membrane damage; the accumulation of soluble sugar, soluble protein and proline; and the occurrence of photoinhibition owing to the damage in photosynthetic apparatus and the decrease in light energy conversion efficiency and electron transfer rate. This effect on pomegranate was much more severe under freezing stress. This study results provide useful information to understand the molecular mechanism of pomegranate response to cold stress and also lay a foundation for the selection of major candidate genes of cold tolerance molecular breeding in pomegranate.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was funded by Key Projects of Science and Technology Research and Development of Xinjiang Production and Construction Corps (2021AB015).

Notes on contributors

Sihui Guan

Sihui Guan is a graduate student in horticulture at Shihezi University's College of Agriculture. Her research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. Her research interests are focused on the study of fruit Tree stress physiology and Molecular Biology.

Yaqian Chai

Yaqian Chai is a graduate student in horticulture at Shihezi University's College of Agriculture. Her research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. Her research interests are focused on the study of fruit Tree stress physiology and Molecular Biology.

Qing Hao

Qing Hao is a researcher in Xinjiang Academy of Agricultural Sciences. His research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. His research interests are focused on the study of Genetic mechanism of pomegranate traits and the breeding of new varieties.

Yadong Ma

Yadong Ma is a teather at Shihezi University's College of Agriculture. His research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. His research interests are focused on the study of horticultural Plant Cultivation and Physiology.

Wenliang Wan

Wenliang Wan is a doctoral candidate at Shihezi University's College of Agriculture. His research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. His research interests are focused on the study of physiology and biochemistry of plant stress.

Huiying Liu

Huiying Liu is a professor in horticulture at Shihezi University's College of Agriculture. Her research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. Her research interests are focused on the study of growth and development of facility crops and its regulation.

Ming Diao

Ming Diao is a professor in horticulture at Shihezi University's College of Agriculture. His research is conducted in the research group “soft seed pomegranate” at the Key Laboratory of Special Fruits and Vegetables Cultivation Physiology and Germplasm Resources Utilization of Xinjiang Production and Construction Corps. His research interests are focused on the study of horticultural plant physiology and ecology.