Comparative transcriptome analysis of cold-tolerant and -sensitive asparagus bean under chilling stress and recovery

Plant material

We selected the asparagus bean cultivars Dubai bean and Ningjiang 3 that are respectively cold-tolerant and -sensitive to chilling temperatures when adults (Huang et al., 2018). Seeds were soaked at 25 °C for 4 h and then sown in watered plugs containing vermiculite and perlite (V:V = 1:1). The plugs were placed in a growth chamber at 25 °C under 14 h light/10 h dark cycles with light intensity of 30,000 lx and watered every two days after seedling emergence. After 2 weeks, growth chamber temperature was reduced to 5 ° C, and light intensity was set to 6,000 lx. Seedling leaves were harvested 0, 3, 12, and 24 h later. Growth conditions were returned to 25 °C and 30,000 lx and seedling leaves were harvested at 3, 12, and 24 h later. Three biological replicates of each cultivar were sampled at each of seven time points.

Cultivar phenotyping

Half of the leaf samples were frozen in liquid nitrogen and stored at −80 °C. Fresh leaves were weighed; proline and superoxide dismutase (SOD) contents in the remaining samples were evaluated using respective assay kits (Cat. Nos. A107-1-1 and A101-1-1; Nanjing Jiancheng Bioengineering Institute, Nanjing, China) as described by the manufacturer.

Extraction of RNA and library construction

Total RNA was extracted using the RNAiso Plus kit (Takara Bio Inc., Kusatsu, Japan) as described by the manufacturer. The integrity and quality of the RNA was confirmed by 1.5% agarose gel electrophoresis and using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA).

Since 21 RNA samples were extracted from each cultivar, we constructed a total of 42 sequencing libraries using NEBNext^® Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs Inc., Ipswich, MA, USA) as described by the manufacturer. Libraries were sequenced on an Illumina Hiseq X platform, and 125–150 bp paired-end reads were generated.

Data analysis

Raw sequencing reads were filtered by trimming adapters and removing low quality reads. Clean reads were mapped to the cowpea genome using HISAT2 software (Kim, Langmead & Salzberg, 2015; Lonardi et al., 2019), and reads mapped to each gene were counted using HTSeq v. 0.6.1. The fragments per kilobase of transcript per million mapped reads (FPKM) of each gene was calculated according to their length and the number of reads mapped to them (Anders, Theodor & Huber, 2015). Correlations between samples were determined using Pearson correlation coefficients (PCC). Principal component analysis (PCA) was performed using Factoextra and FactoMineR package in R software.

Differentially expressed genes (DEGs) between Dubai bean and Ningjiang 3 under cold stress and recovery were analyzed using the DESeq package in R (ver. 1.18.0). Resulting P values were adjusted using the Benjamini and Hochberg approach to control false discovery rates (FDRs). Genes with adjusted P < 0.05, and log2 fold change ≥ 2 were considered differentially expressed.

Gene Ontology (GO) enrichment of DEG sets was analyzed using the BiNGO plugin of Cytoscape (Maere, Heymens & Kuiper, 2005). Transcription factors (TF) were identified in DEGs using iTAK (Zheng et al., 2016). Pathways of gene sets were analyzed using MapMan v. 3.6.0RC1 tool. Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was analyzed using TBtools (Chen et al., 2020).

Quantitative RT-PCR analysis

We validated the RNA-seq results using 15 cold-responsive genes, including 10 AP2/ERF transcription factors, four circadian rhythm-related genes, and ICE1. Table S1 shows the primers that were designed using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/) and synthesized (Invitrogen; Thermo Fisher Scientific Inc.). Quantitative RT-PCR was performed using three technical replicates and SYBR Premix Ex Taq TM II (Takara Bio Inc.) on a Bio-Rad CFX96 Real-Time PCR system (Bio-Rad Laboratories Inc., Hercules, CA, USA). The PCR amplification conditions were 95 °C for 3 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. The relative expression of the selected genes was normalized to that of the cowpea actin gene (GenBank ID: LOC114185033) and calculated using the 2^−ΔΔCT method.

Cultivar phenotyping under chilling stress and recovery

One symptom of chilling injury is a water deficit in shoots due to imbalanced water transport and transpiration (Kang & Saltveit, 2001). Under cold stress at 5 °C, the leaves of both Dubai bean and Ningjiang 3 gradually drooped and wilted (Fig. 1). After 3 h of recovery at 25 °C, the leaves of both varieties continued to droop. At 12 and 24 h of recovery, the leaves of Dubai bean remained drooped and crinkled, whereas most Ningjiang 3 leaves not only stopped wilting, but also showed signs of resumed growth.

Changes in leaf fresh weight and proline and SOD contents

The fresh weight of Dubai bean and Ningjiang 3 leaves progressively decreased under cold stress at 5 °C (Fig. 2). After 3 h of recovery at 25 °C, the leaf fresh weight of Dubai bean and Ningjiang 3 was unchanged and slightly decreased, respectively. At recovery for 12 and 24 h, the fresh weight of Dubai bean leaves continued to decline, whereas that of Ningjiang 3 increased and was only 10% lower than in cold stress for 0 h (Table S2).

Antioxidant enzymes, such as SOD, catalase (CAT), and ascorbate peroxidase (APX), help plants clear excessive reactive oxygen species (ROS) produced at low temperatures (Yu et al., 2006). In addition, compatible solutes, such as proline and carbohydrates, play important roles in cell osmotic adjustment and maintaining membrane integrity under cold stress (Liu et al., 2012). We examined proline and SOD contents after cold stress at 5 °C and after recovery for 3 h at 25 °C (Fig. 2). We found a gradual increase in the proline and SOD content of leaves from both plants, but these values declined or remained unchanged after recovery for 12 and 24 h (Tables S3 and S4). However, proline content at C3, R12 and R24 time points and SOD content at C3, C12, C24 and R3 time points were significantly higher in Ningjiang 3 than in Dubai bean (P < 0.01).

Taken together, we initially inferred that cold tolerance is better for Ningjiang 3 than Dubai bean at the seedling stage based on the phenotype, leaf fresh weight and proline and SOD contents.

Global transcriptome results for Dubai bean and Ningjiang 3

We used Illumina RNA-seq to investigate transcriptome dynamics at seven time points during cold stress and recovery using total RNA isolated from leaves of Dubai bean and Ningjiang 3. Three independent biological replicates of leaves at each time point were analyzed (21 samples each in total).

Over 2 billion high-quality reads (average ∼49 million reads per sample) were generated from each plant at different time points (Table S5) and mapped to the cowpea genome (Lonardi et al., 2019) at a mean rate of 93.72% (range: 89.41%–96.06%). The R package RNASeqPower (Hart et al., 2013) was used to calculate the statistical power of this experiment. A power of 0.65 was found for genes that presented fold-change greater than 2 with a sample size of three (the number of biological replicates used in this study).

The reads mapped to each gene were counted using HTSeq v. 0.6.1, and the FPKM of each gene was calculated (Anders, Theodor & Huber, 2015). Genes with FPKM >1 were considered to be expressed. At seven time points of cold stress and recovery, 16,211–17,961 and 16,174–17,297 genes were respectively expressed in Dubai bean and Ningjiang 3 (Fig. 3A). The expression of 10%–12% of genes was very high (FPKM ≥ 50) at different time points (Fig. 3B). The expression of genes with high (FPKM 10 to ≥ 50), moderate (FPKM 2 to ≤ 10), and low (FPKM 0.1 to ≤ 2), respectively was 27%–34%, 25%–29%, and 28%–34%. For genes showing high/very high expression, more genes were expressed in Dubai bean during cold stress and in Ningjiang 3 during recovery.

Comparison of global transcriptome between two cultivars

We investigated global differences in the Dubai bean and Ningjiang 3 transcriptomes during cold stress and recovery. We calculated Pearson correlation coefficients (PCC) of all samples from seven time points based on the mean FPKM of three biological repeats for all expressed genes (Fig. 4A). The PCC values between the two cultivars under cold stress at 5 °C for 0, 3, 12, and 24 h, were 0.93, 0.96, 0.97, and 0.88, respectively, and 0.94, 0.79, and 0.81 at 25 °C recovery for 3, 12, and 24 h, respectively. The correlation between samples of the two cultivars under 5 °C and recovery at 25 °C for 3 h was close (PCC >0.88) but not at 12 h of recovery at 25 °C (PCC = 0.79).

We conducted PCA using the factoextra and FactoMineR packages in R (Fig. 4B). Samples with closer distances suggested more similar transcriptional programs. Samples under cold stress at 5 °C were distributed on the left side of the plot, while samples under recovery at 25 °C were distributed on the right, indicating that the low-temperature response and recovery at room temperature are different processes. The cultivars under cold stress 5 °C and after 3 h of recovery at 25 °C were closely grouped, but more distanced at 12 h of recovery at 25 °C (Fig. 4A). The results indicated that the two cultivars are more likely to have similar transcriptomes and functions/activities under cold stress at 5 °C and recovery for 3 h at 25 °C. In contrast, transcriptional activity at 12 h of recovery was somewhat distinct and might explain the difference in cold tolerance between the two asparagus bean cultivars.

Comparison of DEGs between Dubai bean and Ningjiang 3

Figure 5A shows the upregulated expression of 18,257 genes (including 745 TF-encoding genes) and the downregulated expression of 15,142 genes (including 572 TF-encoding genes) at different stages of cold stress and recovery in Dubai bean compared with Ningjiang 3. Most (8,954 and 7,455) genes were differentially expressed after recovery for 12 (R12), and 24 (R24) hours, whereas the minimum (1,947) was differentially expressed after cold stress for 3 h (C3).

Transcription factor prediction of DEGs between Dubai bean and Ningjiang 3

Transcription factors (TFs) play a vital role in regulating biotic and abiotic stress response networks (Liu et al., 2019). We investigated differentially expressed TFs to better understand the molecular mechanisms of the low-temperature responses in the two varieties (Fig. 5B). The number of DEGs in some families differed between the two cultivars (Table S6). For instance, more TFs (including bHLH, NAC, WRKY, bZIP, and AP2/ERF) were upregulated in Dubai bean, whereas most downregulated TFs in Dubai bean were in the C2H2, Aux/IAA, TCP, HB-KNOX, and DDT families. Transcription factors in the zf-HD and FAR1 families were upregulated in Dubai bean and Ningjiang 3, respectively.

GO enrichment of DEGs between Dubai bean and Ningjiang 3

We investigated the potential functions of DEGs by analyzing GO enrichment (Fig. 6). A few GO terms were enriched under cold stress, such as responses to auxin and microtubule-based movement. However, more terms, such as protein phosphorylation, photosynthesis, oxidation–reduction process, and cellular glucan metabolic process were enriched during recovery, suggesting their functional participation in cold tolerance of asparagus bean.

Comparison of DEG pathways between Dubai bean compared with Ningjiang 3

We classified the functions of the DEGs using Mercator4 v. 3.0 to further determine functional groups (Lohse et al., 2014). The log2 fold changes of DEGs between Dubai bean and Ningjiang 3 at each stage were loaded into the Overview pathway using the MapMan tool, and the functions of DEGs were classified (Fig. 7A). The top five groups with most DEGs were RNA biosynthesis, protein modification, solute transport, cell cycle organization, and cell wall organization. The ratios of upregulated genes in Dubai bean were essentially stable in the first three groups, fluctuated in cell wall organization, and gradually declined in cell cycle organization.

The ratios of upregulated genes associated with photosynthesis in Dubai bean under cold stress cold stress were >65% at 0, 3, 12 h and 7% at 24 h, 83% within 3 h of recovery, then respectively decreased to 20% and 5% at 12 and 24 h of recovery. The trend was the same for coenzyme metabolism. Some other pathways including chromatin organization, RNA processing, DNA damage response, protein biosynthesis, and protein homeostasis, gradually declined, indicating that these functions were adversely affected in Dubai bean at low temperatures.

We investigated the metabolic pathways responsible for the differences in cold tolerance between Dubai bean and Ningjiang 3 at 12 h of recovery. We loaded the log2 fold changes in DEGs onto the X4.1 Metabolism overview R1.0 pathway using MapMan (Fig. 7B), because the correlation between the two varieties was the lowest (PCC = 0.79) at this stage. The differential activity of metabolic pathways at 12 h of recovery was consistent with the PCA results. Most of the genes involved in photosynthesis, including those associated with photophosphorylation and the Calvin cycle, were downregulated in Dubai bean, indicating impaired photosynthetic function. In addition, most genes involved in cell wall synthesis and modification, amino acid synthesis, and nucleotide metabolism were also downregulated in Dubai bean. These results indicated that the disrupted cell functions of Dubai bean had not yet returned to normal despite 12 h of recovery under normal temperatures.

Identification of DEGs in Dubai bean and Ningjiang 3

We analyzed DEGs in each cultivar at different stages of cold stress and recovery by plotting Venn diagrams using Venny2.1 (Fig. 8). Compared with 0 h, 3773, 6286, and 7362 genes in Dubai bean were differentially expressed after cold exposure for 3, 12, 24 h, respectively. We identified 2,026 DEGs that were common to all three stages. Furthermore, 2,374, 2,894, and 5,527 DEGs were identified in Ningjiang 3 exposed to 5 °C for 3, 12, 24 h, respectively, and 1,636 genes were common to all three stages. We identified 7,005 and 5,994 DEGs that were common to Dubai bean and Ningjiang 3, respectively, under all three recovery stages.

KEGG analysis of DEG enrichment in Dubai bean and Ningjiang 3

We blasted common DEGs into the KEGG database to identify the primary biological pathways in which these common DEGs were involved (Fig. 9). Enriched pathways unique to Dubai bean during cold exposure included chromosome and associated proteins, enzymes with EC numbers, and glycosylphosphatidylinositol (GPI)-anchored proteins. Those in Ningjiang 3 included plant hormone signal transduction, pentose and glucuronate interconversions, and porphyrin metabolism. Enriched KEGG pathways common to both cultivars included the MAPK signaling pathway, plant-pathogen interaction, and transcription factors. Ribosome biogenesis, photosynthesis, chaperones and folding catalysts, and glyoxylate and dicarboxylate metabolism were enriched KEGG pathways in both cultivars during recovery. These annotations provide a basis for understanding the common and different biological pathways involved in the two asparagus bean cultivars under cold stress and recovery.

Enriched transcription factors in Dubai bean and Ningjiang 3

The results of the KEGG pathway enrichment analysis revealed TFs that were enriched in DEGs that were common of both cultivars under cold stress, indicating their involvement in the low-temperature responses of asparagus bean. We identified 36 differentially expressed TFs in this pathway (Table 1); 28 and 25 were in Dubai bean and Ningjiang 3, respectively. These TFs belonged to the MYB (12), AP2/ERF (11), WRKY (6), HSF (3), HD-ZIP (2), E2F/DP (1), and NF-YB (1) families. Homologous genes of TFs, comprising ERF (LOC114194502), MYB (LOC114169018 and LOC114168556), and three WRKY (LOC114194094, LOC114184078 and LOC114188798) are involved in the cold response of Arabidopsis, according to annotations in The Arabidopsis Information Resource (TAIR) database.

Validation by qRT-PCR

We validated the accuracy and reliability of RNA-seq data using 15 cold-responsive genes comprising 10 AP2/ERF genes encoding TFs, four associated with the circadian clock, and ICE1 (Xu et al., 2011a; Xu et al., 2011b; Dong, Farré & Thomashow, 2011; Chinnusamy et al., 2003). The expression of these genes determined by qRT-PCR was similar to those determined using RNA-seq (Fig. S1). The mean correlation coefficient between the RNA-seq and real-time PCR results was 0.83 (0.66–0.97), indicating that the RNA-seq data accurately reflected the abundance of transcripts (Tables S7 and S8).