Identification of common and specific cold resistance pathways from cold tolerant and non-cold tolerant mango varieties

Plant materials and cold stress

Two mango varieties, Jinhuang (JH) and Tainong (TN), were used to analyze the resistance mechanisms against cold stress. A total of 15 seedlings of JH (61.20 cm ± 8.35) and 15 seedlings of TN (57.00 cm ± 6.75) were obtained from Hainan Baizhou Agriculture Co., LTD, Danzhou, China. Each seedling was being added 50 mL of H₂O and was pre-treatment in a greenhouse at 24 °C for 7 days. Subsequently, the 15 seedlings of JH were devided into three groups on average, namely JH4-1 group, JH4-2 group and JH24 group, while the 15 seedlings of TN were devided into three groups on average as well, namely TN4-1 group, TN4-2 group and TN24 group.

Among these groups, JH4-1, JH4-2, TN4-1, and TN4-2 were placed in a low temperature culture room at 4 °C for 36 h (JH4-1 and TN4-1 groups) or 72 h (JH4-2 and TN4-2 groups) for cold stress treatment. The two control groups of JH24 and TN24 were maintained in the greenhouse at 24 °C for 72 h. The leaves of treated samples were collected at 36 h (JH4-1 and TN4-1 groups) or 72 h (JH4-2 and TN4-2 groups). Control samples were also collected from the JH24 group and TN24 group. Five leaf samples were collected from each group. All the samples were immediately frozen in liquid nitrogen and stored at −80 °C for further use.

Superoxide dismutase and peroxidase

Before carrying out transcriptome sequencing, the physiological and biochemical indicators of two mango varieties should be determined. In the study, two time points of superoxide dismutase (SOD) and peroxidase (POD) activities of JH mango and TN mango varieties were used to judge this standard under adverse conditions. Each leaf of 30 mango plants was selected after cold stress, and every leaves from the same group were mixed for evaluating the physiological indices in both mango varieties. The activities of SOD (Catalog No. G0101F, Grace, Suzhou, China) and POD (Catalog No. G0107F, Grace, Suzhou, China) were measured using the kits provided by Suzhou Grace Biotechnology Co., Ltd. The statistical significance of SOD or POD value in different groups was determined by Student’s t test. Values of P < 0.05 were considered significant.

RNA preparation and qualification

Total RNAs from mango leaves were extracted by using RNAprep Pure Plant Plus Kit (Polysaccharides & Polyphenolics-rich) according to the manufacturer’s instruction (Catalog No. DP441, Tiangen, Beijing, China), and further treated with DNase I (Catalog No. RT411, Tiangen, Beijing, China) to remove DNA contamination. The RNA concentration and purity were measured using Nanodrop2000 (Thermo Fisher Scientific, MA, USA). The RNA integrity (RIN) was assessed using the Caliper LabChip GX system (Perkin Elmer, Waltham, MA, USA).

Library preparation for transcriptome sequencing

A total amount of ~3.5 μg RNA per sample was used as input material for the cDNA library construction. Sequencing libraries were generated using Hieff NGS Ultima Dual-mode mRNA Library Prep Kit (Catalog No. 13533ES96, Yeasen, Shanghai, China) as following three steps: firstly, the enrichment, purification, and fragmentation of mRNA were performed; secondly, the first strand cDNA was synthesized for preparing the second strand cDNA, which includes end repair and dA-tailing; thirdly, adaptor primers adpter3 and adpter5 (Table 1) were linked to the two ends of ds cDNA and the products were purified. Finally, the cDNA libraries were further amplified and purified for the library quality inspection by using Qsep400 standard DNA clip kit in Qsep-400 instrument.

The prepared libraries were sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) using NovaSeq 6000 S4 Reagent Kit (Illumina, San Diego, CA, USA), and 150-bp paired-end reads were generated (Modi et al., 2021).

Data quality control of the transcriptomic data

Based on sequencing-by-synthesis (SBS) technology, the significant amounts of raw data with high-quality were generated. Raw data was saved in fastq format. Raw reads in the fastq format were firstly cleaned to remove the adaptor sequences and low quality reads (ploy-Ns ratio greater than 10%, or bp value of Q ≤ 10 greater than 50%) for generating the clean data. Q30 and GC content of the clean data were calculated to evaluate the overall quality of the clean reads. All the downstream analyses were based on the clean data with high quality (Liu et al., 2021a).

Mapping analysis of the transcriptomic sequencings

Reference genome of mango was downloaded for the RNA-Seq mapping analysis. Version information of the reference genome is Mangifera_indica.v4.0.genome.fa (Wang et al., 2020). HISAT2 (Kim et al., 2019), a Burrows-Wheeler Transform and Ferragina-Manzini (FM) index based search, was used to map clean reads in the reference genome to obtain the localization information of reads on the mango genome. StringTie was applied to assemble the mapped reads for subsequent analysis (Pertea et al., 2015). It utilizes a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantify transcripts representing multiple spliced variants for each gene locus.

Identification of differentially expressed genes

Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) was applied to measure the expression level of a gene or a transcript by StringTie using maximum flow algorithm. In the study, the differential expression genes (DEGs) analysis was performed on comparing the JH4-1 group with the JH24 control group and comparing the JH4-2 group with the JH24 group, respectively. Meanwhile, the DEGs analysis was also performed on comparing the TN4-1 group with the TN24 control group and comparing the JH4-1 group with the JH24 group, respectively. The four comparing groups were named JH24 vs. JH4-1, JH24 vs. JH4-2, TN24 vs. TN4-1, TN24 vs. TN4-2.

In the study, transcripts that increased or decreased with a fold change (FC) ≥ 2 and false discovery rate (FDR) < 0.01 were considered to be differentially expressed. DESeq2 accepted input of the clean reads and the DEGs between the experimental group and control group were screened (Liu et al., 2021b). The hierarchical clustering map was used to show the distribution of DEGs, and cluster analysis was used to judge the expression pattern of each gene.

Verification of transcription levels by real-time quantitative reverse transcription PCR

To confirm gene expression differences of four comparing groups (JH24 vs. JH4-1, JH24 vs. JH4-2, TN24 vs. TN4-1, TN24 vs. TN4-2) by transcriptomic sequencing, four DEGs of probable xyloglucan endotransglucosylase/hydrolase protein 23 (XTH23), probable terpene synthase 12 (TS12), BON1-associated protein 2-like (BON1-2), RADIALIS-like 4 protein (RADIALIS-4) were selected for real-time quantitative reverse transcription PCR (qRT-PCR) analysis using the primers in Table 1. Total RNA was extracted from plant leaves by using a FastPure Universal Plant Total RNA Isolation Kit (Vazyme, Nanjing, China) according to the manufacturer’s instructions. Two μg RNA (360~500 ng/μL) per sample was transcribed into first strand cDNA by using a HiScript III 1st Strand cDNA Synthesis Kit (+gDNA wiper) (Vazyme, Nanjing, China). Five leaf samples of each group were mixed into a mixture, and three technical replicates were performed for each mixed sample. The amplification program was conducted on a StepOne real-time PCR system (Applied Biosystems, Waltham, MA, USA) and the cycle condition was as follows: 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s, and 60 °C for 1 min. The internal control gene, β-actin, was used to normalize each gene. Transcription levels of selected DEGs were evaluated based on the 2^−ΔΔCt method. The statistical significance of qRT-PCR was determined by Student’s t test. Values of P < 0.05 were considered significant.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses of DEGs

GO enrichment analysis of four comparing groups (JH24 vs. JH4-1, JH24 vs. JH4-2, TN24 vs. TN4-1, TN24 vs. TN4-2) was conducted by the topGO to analyze the functional classification of DEGs. GO terms were extracted from the best hits obtained from BLASTx against the non-redundant database using Blast2GO (Conesa et al., 2005). Then the obtained DEGs were sorted by GO categories using in-house Perl scripts. GO terms with a P value below 0.05 were considered significantly enriched.

Meanwhile, the KEGG database (http://www.genome.jp/kegg/genome) was used to analyze the molecular pathway of these four comparing groups, and DEGs with a P value below 0.05 in the KEGG pathways considered to be significantly enriched using ClusterProfiler software (Wu et al., 2021).

Phylogenetic tree construction of ERFs

ERF is a class of protein belonging to the AP2/ERF transcription factor superfamily involved in the regulation of biotic and abiotic stresses in plants. In order to explore the up-regulated ERFs (log₂FC > 2) in cold-tolerant mango varieties, all ERF-related genes with differential expression were screened from JH24 vs. JH4-2 group. Phylogenetic tree was constructed based on the amino acid sequences to determine the phylogenetic relationship of all ERFs in JH mango under cold stress. The protein informations of 43 ERFs were listed in Table S1. All amino acid sequences were aligned using ClustalX, then passed to MEGA 11 for tree building using maximum likelihood method with 1,000 bootstrap replicates. Phylogenetic analysis was conducted in MEGA11 (Tamura, Stecher & Kumar, 2021).

The activities of SOD and POD of mango plants under cold stress

The phenotypic observation showed that JH mango plants of both JH4-1 group and JH4-2 group displayed asymptomatic damage, but leaf tips displayed browning after transferred back at room temperature (RT) for 3–5 days (Fig. 1A). Meanwhile, no obvious symptoms were observed on TN mango plants of TN4-1 group or TN4-2 group, but leaf bases displayed wither and rot at 3–5 days after being transferred back at RT (Fig. 1B).

At physiological levels, JH mango plants performed better than TN mango plants responses to the low temperature treatment. For instance, comparing TN4-1 group with JH4-1 group, and comparing TN4-2 group with JH4-2 group, the activities of SOD were found to be higher in JH mango plants than that of TN mango plants under cold stress. However, there was no significant difference in SOD activities between the two mango varieties at 24 °C for 72 h (Fig. 1C). Similar results were observed in the POD activities of the two mango varieties. In detail, the activities of POD were found to be significant differences between the JH mango plants and the TN mango plants at 4 °C treatment for 36 h or 4 °C treatment for 72 h, while no significant difference of the two groups at 24 °C for 72 h (Fig. 1D). These findings showed that JH mango is a more cold-tolerant variety than TN and could mitigate the adverse effects of low temperature by increasing the activities of antioxidant defense enzymes.

Transcriptomic data quality and mapping analysis

The high-quality of 30 total RNA samples from mango leaves (Conc = 150.60~514.70 ng/μL; OD_260/280 = 2.01~2.15; RIN = 8.10~8.70) were obtained and conformed to the cDNA library preparation requirement. A large number of raw reads were generated from an Illumina NovaSeq 6000 platform and about 20,000,000 clean reads of each sample were further generated. In total, the GC content was 42.03~43.85% and Q30 percentage was over 93.10% for each sample (Table 2), indicating that high quality transcriptomic sequencing clean data was obtained. The clean data of each sample is available in the GenBase in National Genomics Data Center, Beijing Institute of Genomics, Chinese Academy of Sciences/China National Center for Bioinformation, under accession number CRA016366 that is publicly accessible at https://ngdc.cncb.ac.cn/gsa/browse/CRA016366 (CNCB-NGDC Members & Partners, 2024).

Screening and cluster analysis of the transcriptomic DEGs

After low temperature treatment, a large number of DEGs were identified in both varieties. For the JH mango variety, 8,337 DEGs (4,051 up-regulated and 4,286 down-regulated) were identified from the JH24 vs. JH4-1 using DESeq2 software based on FC and FDR values, while 10,683 DEGs (5,552 up-regulated and 5,131 down-regulated) were further identified in JH24 vs. JH4-2, indicating more DEGs were screened when treated at 4 °C for longer duration (Fig. 2 and Table 3). Further analysis revealed that both the highest up-regulated genes were XTH23, with log₂FC of 12.93 and 13.93, respectively. The highest down-regulated genes were both of TS12, with log₂FC of −11.82 and −10.47, respectively. The qRT-PCR showed the consistent results (Fig. 3).

For the TN mango variety, 7,996 DEGs (3,763 up-regulated and 4,233 down-regulated) were identified from TN24 vs. TN4-1, while 10,723 DEGs (5,674 up-regulated and 5,049 down-regulated) were further identified in TN24 vs. TN4-2, also indicating more DEGs were screened when treated at 4 °C for longer duration (Fig. 2 and Table 3). Further analysis revealed that the highest up-regulated gene was XTH23 (log₂FC of 13.04) in TN24 vs. TN4-1, while it was BON1-2 (log₂FC of 13.91) in TN24 vs. TN4-2. The highest down-regulated genes were both of RADIALIS-4, with log₂FC of −11.03 and −10.97, respectively. This outcome was also confirmed by the qRT-PCR (Fig. 3).

GO and KEGG enrichment analyses of DEGs

Functional analysis of the DEGs showed that 31 GO terms (q value < 0.05) were annotated in JH24 vs. JH4-1, and nine of them associated with photosynthesis or chloroplast. Furthermore, 62.5% (10/16) of the GO terms were annotated into photosynthesis or chloroplast pathways in JH24 vs. JH4-2 when JH treated at 4 °C for longer duration. In TN24 vs. TN4-1, nine of them (9/46, 19.6%), the same GO terms as JH24 vs. JH4-1, were annotated. In addition, eight (8/16, 50.0%) photosynthesis or chloroplast associated GO terms were identified when TN treated at 4 °C for longer duration (Figs. 4 and 5). These results indicate that the most annotated photosynthesis-related or chloroplast-related GO terms respond to the cold stress in mango plants.

In general, more pathways were enriched and annotated when two varities were treated at low temperature for longer duration. KEGG enrichment annotation analysis was performed on DEGs from different comparison groups of JH24 vs. JH4-1, JH24 vs. JH4-2, TN24 vs. TN4-1, and TN24 vs. TN4-2. The results showed that 132, 134, 131 and 133 pathways were annotated, respectively (Fig. 6). Furthermore, three significant pathways were found in JH24 vs. JH4-1 by enrichment analysis, namely photosynthesis-antenna proteins pathway, photosynthesis pathway, and valine, leucine and isoleucine degradation pathway. In JH24 vs. JH4-2, in addition to annotating to these three pathways, four KEGGs of thiamine metabolism pathway, inositol phosphate metabolism pathway, glycine, serine and threonine metabolism pathway, and glyoxylate and dicarboxylate metabolism pathway were annotated. Similarly, three enrichment pathways were found in TN24 vs. TN4-1, including photosynthesis pathway, photosynthesis-antenna proteins pathway and proteasome pathway. More enrichment pathways were found in TN24 vs. TN4-2, including oxidative phosphorylation pathway, ubiquitin mediated proteolysis pathway and inositol phosphate metabolism pathway when treated at 4 °C for longer duration (Fig. 7).

Common and specific cold resistance pathways

Two common pathways, photosynthesis-antenna proteins pathway and photosynthesis pathway, were found by comparing the KEGG enrichments from different comparison groups of JH24 vs. JH4-1, JH24 vs. JH4-2, TN24 vs. TN4-1, and TN24 vs. TN4-2. Further analysis of the annotated pathways showed that they have the same up-down genes mode pattern, suggesting that these two pathways are the primary protection response pathways to cold resistance in mango (Figs. S1 and S2). However, in addition to these two common pathways, it was found that the cold-resistant JH mango was annotated to a meaningful pathway, valine, leucine and isoleucine degradation pathway when treated with low temperature for 36 h. Furthermore, more specific pathways were identified and obtained in JH when treated with low temperature for 72 h, such as thiamine metabolism pathway, inositol phosphate metabolism pathway, glycine, serine and threonine metabolism pathway, and glyoxylate and dicarboxylate metabolism pathway (Table 4), indicating that JH mango variety has the specific cold resistance regulation pathways to cold stress.

Phylogenetic relationship of ERFs in JH mango

Compared with the control group, a total of 43 ERFs genes were significantly annotated into DEGs in JH mango after cold-treated for 72 h. The log₂FC values of these DEGs were ≥2.0, and three of them ERF109-1, ERF017-1 and ERF017-2 were ≥10.0. In addition, the log₂FC values of 24 ERFs were 5~10, and the other 16 ERFs were ≤5. The phylogenetic tree showed that the homologous genes were classify into a sub-cluster, i.e., ERF105-1~6, ERF109-1~6, ERF004-1~5, indicating homologous proteins have similar role or complementary function (Fig. 8). Interestingly, the up-regulated levels of ERF105-1~6 were relatively high and consistent, with log₂FC values ranging from 5.52~7.60. Furthermore, the expression levels of ERF017-1~3 were highly up-regulated, and the log₂FC value of ERF017-1, ERF017-2 were 12.14 and 10.75, respectively. These significantly up-regulated ERFs genes may play important regulatory roles in plant to cold resistance.