Full text
Vicia sativa, an annual legume species, is generally cultivated worldwide due to its stronger survivability in extreme conditions. As an economically important legume, V. sativa is commonly used as fodder for livestock and as green manure to improve soil fertility (Dong et al., 2019; Thomas et al., 2022; Zhang, Li, et al., 2020); for instance, the use of V. sativa can reduce soil erosion processes during the first stages of vineyard plantation through soil quality improvements and by reducing soil and water losses (Rodrigo-Comino et al., 2020). In addition, V. sativa has been utilized in health-promoting foods for human consumption as it contains high levels of proteins and other important nutraceutical components, including vitamins, β-carotenes, flavonoids, catechins, and condensed tannins (Fu et al., 2020; Pastor-Cavada et al., 2011). As a widespread stressor, cold stress seriously affects the growth and reproduction of crops in high-altitude and temperate regions (Kazemi-Shahandashti & Maali-Amiri, 2018; Min, Liu, et al., 2020). Although V. sativa is well adapted to various environmental stressors, cold stress is a critical factor that affects its utilization (Shi et al., 2020; Zhou et al., 2022). Therefore, determining the molecular mechanism of V. sativa's response to cold stress is important for enabling breeders to cultivate cold-resistant varieties.
Obvious physiological and cellular perturbations occur in many plant species under cold stress, and plants usually continue to survive in harsh environments by maintaining cell and membrane stability and protein structure (Chen et al., 2018). The three main cold-responsive genes in plants are Inducer of CBF Expression (ICE), C-repeat Binding Factors (CBFs), and Cold-Regulated Genes (CORs); these three aforenamed key players model an imperative signaling pathway, the ICE-CBF-COR cascade, which is a cold response pathway that alleviates cold stress in plants (Hwarari et al., 2022). Previous studies have reported that the ICE-CBF-COR signaling pathway is considered to be important during the cold acclimation process (Liu et al., 2018; Ritonga & Chen, 2020; Wi et al., 2022). As a pioneer of cold acclimation, ICE acts upstream to induce and regulate the expression of the CBFs, which are members of the APETALA 2/Ethylene Response Factor (AP2/ERF) transcription factors (TFs) family, and promote the widespread repair of cellular metabolic pathways and the reconstruction of tissue structure during cold acclimation (Kang et al., 2016; Leuendorf et al., 2020). Furthermore, CBFs can bind to cis-elements (CCGAC) in the promoter of COR genes and activate their expression (Song et al., 2021).
In addition to the ICE-CBF-COR genes, the response of plants to cold stress is also regulated by other factors, such as hormones, circadian rhythm, and light (Ritonga & Chen, 2020). For example, the main response of Lolium perenne to cold acclimation (3°C for 7 days) and freezing stress (−7°C for 4 days) was elevated abscisic acid (ABA), jasmonic acid (JA), and salicylic acid (SA) in different tissues of three genotypes, including sensitive, highly tolerant, and semi-tolerant genotypes (Prerostova et al., 2021). Similar results appeared in Oryza sativa (Wang et al., 2018), Arabidopsis thaliana (Wu et al., 2019), Camellia sinensis (Zhou et al., 2020), and Artemisia annua (Liu et al., 2017), while the opposite trend was also shown in other species or tissues, such as Zea mays and Vitis vinifera (Gao-Takai et al., 2019; Pal et al., 2020). Previous studies have suggested that ABA, JA, and SA can enhance cold resistance in multiple plants through the regulation of CBF transcription, sucrose metabolism, and soluble sugar content (Hu et al., 2017; Lee & Seo, 2015; Zhao et al., 2021). Additionally, light is a critical environmental signal that regulates plant growth, development, and stress response. In addition to red light, recent research has uncovered a new mechanistic framework in A. thaliana; freezing tolerance in Arabidopsis was enhanced via a blue-light signal mediated by CRYPTOCHROME2 (Li, Zhong, Qu, et al., 2021).
Compared with the microarray method, high-throughput sequencing has significantly improved the efficiency of transcript collection, which has been used to reveal the underlying regulatory networks of multiple species in response to cold stress, such as Medicago sativa (Zhou et al., 2018), Zanthoxylum bungeanum (Tian et al., 2021), Picea abies (Vergara et al., 2022), and Fragaria nilgerrensis (Liu et al., 2021). In a recent study, high-throughput transcriptome sequencing revealed differences in the molecular mechanisms of the cold stress response between V. sativa's leaves and roots, and the overexpression of six cold-resistant candidate genes enhanced cold tolerance in yeast (Min, Liu, et al., 2020). Moreover, a protocol of ethyl methane sulfonate (EMS)-induced mutations in V. sativa was established, and two mutants without anthocyanin accumulation showed increased cold tolerance (Shi et al., 2020). These results lay the foundation for understanding the underlying mechanisms and novel genes involved in V. sativa's response to cold stress. In this study, a total of three cold-resistant and three cold-sensitive V. sativa accessions were selected to investigate the difference in transcription and metabolism changes between cold-resistant and cold-sensitive accessions during cold stress. The present study will help to reveal the molecular regulation mechanism of V. sativa's response to cold stress and provide new candidate genes for the cultivation of cold-resistant V. sativa.
MATERIALS AND METHODS Plant materials and stress treatmentsA total of six V. sativa accessions were used in the present study, including three cold-resistant (Lan1, Lan2, and Lan3) and three cold-sensitive (368, 521, and 538) accessions. Lan1, Lan2, and Lan3 were cultivated on the Qinghai Tibet Plateau in our laboratory; the average altitude in this area is more than 4000 m, and the annual average temperature is less than 0°C. Among the 541 V. sativa accessions provided by the National Plant Germplasm System of the United States, we found that nine accessions grew slowly on the Qinghai Tibet Plateau. Additionally, a hydroponic experiment showed that 368 (PI 383798), 521 (PI 173158), and 538 (PI 173160) had the highest wilting degrees after cold treatment for 24 h (Figure S1a), and the weight loss ratios of these three accessions were the highest (Figure S1b). Upon collecting the experimental samples in this study, we performed a hydroponic experiment using three cold-resistant and three cold-sensitive V. sativa accessions (Figure S2). Before germination, the seeds were vernalized at 4°C; the seedlings (3 days old) were transplanted into a 1/2 MS nutrient solution (pH 5.8) and grown in an artificial climate incubator at 22°C (16 h light/8 h dark). After growing in the incubator for 2 weeks, the seedlings were removed from the incubator and grown at 4°C (16 h light/8 h dark). A total of six cold treatments were conducted in this study, including five cold treatments (for durations of 3, 6, 12, 24, and 48 h) and one control (cold treatment for 0 h).
Physiological experimentsFor each cold treatment duration, leaf materials were quickly harvested from 60 V. sativa seedlings, and were divided into three biological replicates for the physiological experiments. Some fresh leaf materials of each sample were used to determine the chlorophyll content and electrolyte leakage, as described previously (Gitelson et al., 2003; Su et al., 2015); the remaining leaf materials were quickly flash-frozen in liquid nitrogen and stored in a freezer at −80°C until other physiological indicators were determined, such as the contents of malondialdehyde (MDA), hydrogen peroxide (H2O2), soluble sugar, and proline (Pro). These four physiological indicators were measured according to the method described in a previous report (Zhou et al., 2018). SPSS software (19.0, Chicago, Illinois, USA) was used to perform one-way analysis of variance (ANOVA) on all data, and the differences between treatments were assessed using Duncan's multiple range test at a significance level of 0.05.
RNA extraction and cDNA library constructionThree seedlings (with each replicate containing one seedling) of six V. sativa accessions were collected after cold treatments (0, 6, and 48 h), and the shoot tissues were pooled into a frozen tube pipe for RNA extraction. According to the manufacturer's instructions, the TRIzol method was used to isolate the total RNA of each sample, and RNA concentration was determined using a NanoDrop ND8000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Only samples with RNA concentrations greater than 600 ng/μL were submitted for transcriptome sequencing. An mRNA-Seq Sample Preparation Kit (NEBNext® Ultra™ RNA Library Prep Kit for Illumina®, NEB, USA) was used to construct 54 cDNA libraries in this study, according to the manufacturer's instructions. All cDNA libraries were sequenced by the Beijing Biomarker Institution (
In order to obtain clean reads, sequencing read filtering was performed to remove some low-quality reads from the raw reads, such as reads with high unknown base (N) content and that were adaptor-polluted. Afterward, Trinity software was used to assemble the transcriptome sequences (Grabherr et al., 2011), and all the parameters used default values. Moreover, the functional annotation of the transcriptome was performed based on five public databases, including the NCBI nonredundant protein sequences (NR), Clusters of Orthologous Groups of proteins (COG), Gene Ontology (GO), Eukaryotic Orthologous Group (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases.
qRT-PCR verificationA total of ten unigenes were randomly selected to verify the RNA-seq results in the present study using qRT-PCR analysis. qRT-PCR was performed using a CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, CA, USA). The cDNA sequences of all genes were synthesized using 2 × SG Fast qPCR Master Mix (Sangon Biotech, Shanghai, China), and CFX Manager software (Bio-Rad, CA, USA) was used to collect qRT-PCR data. The qRT-PCR parameters were as follows: 95°C for 3 min, 39 cycles of 95°C for 10 s, and 55°C for 30 s. In this study, Primer6 software (Premier Biosoft International, Palo Alto, CA, USA) was used to design the specific primers (Table S1). According to the 2−ΔΔCt method, the relative expression levels of all genes were normalized to the Vsactin gene (Min, Lin, et al., 2020).
Differentially expressed genes analysisIn the present study, RSEM software was used to identify the transcript expression levels, which were calculated using the fragments per kilobase of transcript per million mapped transcripts (FPKM) method. Moreover, the NOISeq package was used to identify the differentially expressed genes (DEGs) for each treatment. DEGs were restricted, with an absolute fold change value of ≥2 and an adjusted p-value of ≤0.05 as the thresholds; this was determined via pairwise comparisons of cold-resistant samples with cold-sensitive samples for each treatment duration. PCA was performed to observe differences in the expression levels of unigenes among different samples, and the pheatmap R package was used to map the expression patterns of DEGs during cold treatment. Based on the GO and KEGG databases, the biological functions of the DEGs were annotated through GO and KEGG pathway enrichment analyses. Transcription factors (TFs) in the DEGs were identified according to the regulations described in the PlnTFDB database. The R package WGCNA V4.0.3 was used to further search the cold resistance genes of V. sativa via weighted gene co-expression network analysis (WGCNA) (Langfelder & Horvath, 2008).
Widely targeted metabolome analysisA total of six replicate samples each of the shoot tissues of “Lan3” and “521” were collected after cold treatment for 0, 6, and 48 h, and widely targeted metabolome analysis was performed using a liquid chromatography-mass spectrometry (LC-MS) system provided by the Beijing Biomarker Institution (Beijing, China). The LC-MS system comprised an Acquity I-Class PLUS ultra-high-performance liquid chromatography tandem Xevo G2-XS QT with a high-resolution mass spectrometer. According to the standard procedure, MassLynx v4.2 software was used to collect the raw data, and Progenesis QI software was used for peak extraction and peak alignment. Afterward, all metabolites were identified using the METLIN database and a self-built database, which was provided by the Beijing Biomarker Institution. In order to observe differences in metabolic composition among the different samples, PCA was performed using MetaboAnalyst 5.0 (
Based on the expression patterns of DEGs during cold treatment, we found that most genes were induced during cold treatment; thus, five cold resistance candidate genes were selected to test their cold resistance in transgenic yeast, as their expression levels differed between cold-resistant and cold-sensitive materials, including ERFs (c20836.graph_c0 and c54604.graphc0), DREB (c48739.graphc0), TIFY 5A (c41936.graphc0), and TSUD_290770 (c51037.graph_c0) proteins. Primer6 software was used to design the specific primers (Table S1), and the cDNA sequences of these genes were amplified using PrimeSTAR® Max DNA Polymerase (Takara, Beijing, China). The full-length cDNA sequences were cloned into pYES2 (Invitrogen, Carlsbad, USA), which is a yeast expression vector. According to the lithium acetate method described in a previous study (Gietz & Woods, 2002), we transformed an empty pYES2 and five expression vector plasmids in the yeast strain INVSc1. Moreover, freezing stress evaluation was performed according to a previously described method (Min, Liu, et al., 2020).
RESULTS Physiology assayIn order to investigate the differences in the physiological changes between cold-resistant and cold-sensitive V. sativa accessions during cold stress, a total of six physiological indexes were determined, including ion leakage and the contents of chlorophyll, H2O2, MDA, Pro, and soluble sugar. As shown in Figure 1, ion leakage and the contents of chlorophyll, H2O2, MDA, and Pro under cold treatment were higher than those under normal temperature, and significant increases were observed during the 6 h cold treatment; however, the soluble sugar content increased gradually during the cold treatment. Moreover, ion leakage and the contents of H2O2 and MDA in cold-sensitive accessions were significantly higher than in cold-resistant accessions during the cold treatment (Figure 1a–c), and the opposite trend appeared for the other three indicators (Figure 1d–f). The leaves of the three cold-resistant accessions were not significantly different from those under normal conditions, while the leaves of the three cold-sensitive accessions showed severe wilting and curling after 6 h of cold treatment (Figure S2); this was consistent with the change trends of some physiological indexes, namely ion leakage and the contents of H2O2 and MDA.
Enlarge this image.FIGURE 1. Analyses of dynamic physiological responses under continuous cold stress. (a) H2O2 content. (b) Ion leakage. (c) MDA content. (d) Chlorophyll content. (e) Proline content. (f) Soluble sugar content. Each value represents the mean of three replicates ± SD (standard deviation), shown by a vertical error bar. “*”, “**” and “***” above the bars indicate significant differences at the 0.05, 0.01, and 0.001 levels according to Duncan's multiple range test, respectively.
For high-throughput RNA-seq, a total of 54 cDNA libraries were constructed in the present study. A total of 409.78 Gb of clean data were obtained, and the percentage of Q30 in each library was higher than 93.79%. All RNA-seq data in our study were uploaded to the NCBI SRA database (SRR20083627–SRR20083644). Finally, a total of 85,707 unigenes were assembled, with average and N50 lengths of 905.84 and 1835 bp, respectively. Then, the clean reads were mapped onto the assembled unigenes using Bowtie2 software. Overall, the average mapped reads and mapped ratio were 21.70 million and 85.51% of all libraries, respectively (Table S2). Based on the five public databases (NR, GO, KOG, COG, and KEGG), the molecular functions of all unigenes were annotated using BLASTx (E-value ≤10−5). The results of 40,663 (47.44%), 25,584 (29.85%), 23,919 (27.91%), 13,771 (16.07%), and 15,084 (17.60%) unigenes were successfully annotated in the NR, GO, KOG, COG, and KEGG databases, respectively (Figure S3). A total of 5669 (6.61%) unigenes were annotated in all five databases, and 43,962 (51.29%) unigenes were annotated in at least one database. In addition, we randomly selected ten unigenes to verify the reliability of the transcriptome data using a qRT-PCR experiment (Table S1), and the qRT-PCR results of these unigenes were generally consistent with our transcriptome data (Figure S4).
Identification and expression pattern analysis of DEGsThe PCA plot showed high reproducibility between the biological replicates, and all samples showed significant differences between the three cold treatment durations (Figure 2a). In order to further understand the transcriptomic characteristics of V. sativa in response to cold stress, all the unigenes of the cold-resistant and cold-sensitive accessions under each treatment duration were compared. A total of 7745 unigenes were differently expressed between the cold-resistant and cold-sensitive accessions during cold stress, and 4175, 4737, and 4657 DEGs were obtained during cold treatments that lasted for 0, 6, and 48 h, respectively. According to the expression pattern of these DEGs, which were divided into two categories, including up-regulated and down-regulated DEGs (Figure 2b,c), and a total of 1812 unigenes were differently expressed between the cold-resistant and cold-sensitive accessions for these three durations (Figure 2d). Additionally, a total of 6690, 6788, 6033, 7260, 5302, and 6542 DEGs were obtained in Lan1, Lan2, Lan3, 368, 521, and 538 during cold treatment (Figure S5a), respectively, and there were 4177 DEGs among the three cold-resistant accessions (Figure S5b) and 3694 DEGs among the three cold-sensitive accessions (Figures S5c and S6).
Enlarge this image.FIGURE 2. Summary of the differentially expressed genes. (a) Score plot of the PCA analysis. (b) Heat map of all DEGs. (c) A summary of the numbers of up- and down-regulated DEGs after 0, 6, and 48 h of cold stress. (d) Venn diagram represents the number of overlapping DEGs between cold-resistant (CR) and cold-sensitive (CS) V. sativa germplasms.
In the present study, a total of 55 GO categories were assigned to the 7745 DEGs, including 16 cellular component (CC), 12 molecular function (MF), and 27 biological process (BP) categories (Figure S6). Under the CC category, “membrane” was the largest group (1896, 41.5%), followed by “cell” (1815, 39.7%), “cell part” (1792, 39.2%), and “membrane part” (1578, 34.5%). Moreover, the DEGs associated with “catalytic activity” (2503, 54.8%) and “binding” (2140, 46.8%) represented the two most abundant categories in the MF group. Furthermore, “metabolic process” (2516, 55.0%), “cellular process” (2278, 49.8%), and “biological regulation” (723, 15.8%) were more enriched than other terms in the BP group. Additionally, GO category enrichment analysis was performed to determine whether the accumulated transcripts were functionally involved in cold response processes. As a result, “integral component of membrane,” “catalytic activity,” and “metabolic process” were the most significantly enriched categories in BP, MF, and CC, respectively (Figure 3), and “metabolic process” and the “catalytic activity” categories contained numerous genes and had a small p-value.
Enlarge this image.FIGURE 3. GO enrichment analysis of the top 20 most strongly represented categories. The names of the GO categories are listed along the y-axis. The degree of GO enrichment is represented by the −Log10(p-value) (column chart) and the number of transcripts (line chart) enriched in each category.
Interestingly, similar results also appeared in the 4177 DEGs among the three cold-resistant accessions, and 3694 DEGs among the three cold-sensitive accessions (Figure S7a,b). In the “metabolic process” category, there were 664 DEGs shared in the three types of DEG (7745, 4177, and 3694 DEGs), and the UDP-glucosyltransferase family (Fei et al., 2021), LRR receptor-like kinase family (De la Rosa et al., 2020), and cytochrome P450 family 71 protein (Ciacka et al., 2022) contained the largest number of genes. In the “catalytic activity” category, there were 643 DEGs shared among the three types of DEG, and the UDP-glucosyltransferase family (Fei et al., 2021), GDSL-like lipase/acylhydrolase (Dong et al., 2017), and LRR receptor-like kinase family (De la Rosa et al., 2020) contained the largest number of genes. It is important to note that the UDP-glucosyltransferase family and LRR receptor-like kinase family genes in the “metabolic process” and “catalytic activity” categories were identical. As shown in Figure S8, most of the UDP-glucosyltransferase family, LRR receptor-like kinase family, and cytochrome P450 family 71 protein genes were induced by cold treatment, while most of the GDSL-like lipase/acylhydrolase genes showed the opposite trend.
Identification of TFsAs important upstream regulators, TFs can directly regulate the expression level of stress-responsive genes and play critical roles in plant responses to abiotic stresses. In the present study, a total of 241 DEGs were detected as TFs (Table S3), which were classified into 27 families (Figure 4a). Among these TF families, the MYB family contained the largest number of members (56 DEGs), followed by the bHLH (43 DEGs), AP2/ERF (Hu et al., 2017), zf-CCCH (Hou et al., 2022), and zf-C2H2 (Gao-Takai et al., 2019) families. Moreover, we investigated the expression levels of MYB and bHLH TFs during cold treatment, and the expression levels of most of the TFs were induced during cold treatment (Figure 4b,c). Additionally, we found that 24 MYB and 14 bHLH TFs underwent steady upregulation during cold treatment.
Enlarge this image.FIGURE 4. Summary of the transcription factors in DEGs. (a) Distribution of transcription factors in DEGs between the cold-resistant and cold-sensitive accessions during cold stress. (b) Heat map of the MYB transcription factors identified in the DEGs. (c) Heat map of the bHLH transcription factors identified in the DEGs. The orange, red, and light blue in the heat map represent the expression levels of transcription factors.
A total of 2734 metabolites were detected using untargeted LC-MS analysis, and the PCA plot shows that there were significant differences between the cold-resistant (Lan3) and cold-sensitive (521) V. sativa accessions (Figure 5a). All metabolome analysis data were uploaded to the MetaboLights database (MTBLS7363). A total of 1729 differential metabolites were obtained between cold-resistant and cold-sensitive V. sativa accessions, and their relative concentrations showed obvious differences (Figure 5b). Moreover, 969, 1042, and 1087 differential metabolites were found after cold treatment for 0, 6, and 48 h, respectively, and 448 differential metabolites were common to these three durations (Figure 5c,d). Among these 1729 differential metabolites, a total of 253 metabolites were identified (Table S4), and these metabolites were assigned to 78 KEGG pathways (Table S5).
Enlarge this image.FIGURE 5. Summary of the differential metabolites between “Lan3” and “521.” (a) Score plot of the PCA analysis. (b) Heat map of all differential metabolites. (c) A summary of the numbers of up- and down-regulated differential metabolites after 0, 6, and 48 h of cold stress. (d) Venn diagram represents the number of overlapping differential metabolites between “Lan3” and “521.”
In order to annotate the biological functions of DEGs and differential metabolites between “Lan3” and “521,” conjoint KEGG pathway enrichment analysis was performed using the BMKcloud platform (
Enlarge this image.FIGURE 6. The conjoint KEGG pathway enrichment analysis of DEGs and differentially metabolites between cold-resistant and cold-sensitive V. sativa germplasms after 0 h (a), 6 h (b), and 48 h (c) of cold treatment.
A total of 18 differential metabolites and 61 DEGs were obtained in the “phenylpropanoid biosynthesis” pathway (Figure 7, Table S6). Among these 18 metabolites, the contents of naringenin, caffeic acid, p-coumaric acid, and 4-hydroxybenzaldehyde in “Lan3” were significantly higher than in “521” under three cold treatment durations. Significantly, p-coumaric acid and naringenin were shown to be key metabolites in the flavonoid biosynthesis pathway. In addition, the Pearson correlation coefficient between 18 differential metabolites and 61 DEGs related to phenylpropane biosynthesis was calculated. As shown in Figure S11, a total of one, seven, six, and nine DEGs were significantly correlated with naringenin, caffeic acid, p-coumaric acid, and 4-hydroxybenzaldehyde, respectively. Thus, we speculate that these genes may regulate the phenylpropane metabolic pathway of V. sativa.
Enlarge this image.FIGURE 7. The contents of 18 differential metabolites related to the biosynthesis of phenylpropanoids in “Lan3” and “521.” The colors indicate the abundance of metabolites calculated as Log2(Value).
Weighted gene co-expression network analysis was performed to understand the variation in the regulatory network of phenylpropane biosynthesis between cold-resistant and cold-sensitive V. sativa accessions. As a result, a total of 14 distinct modules were identified based on their expression patterns. According to the color of each module in the figure, these modules were named turquoise, blue, brown, yellow, green, red, black, pink, magenta, purple, green-yellow, tan, salmon, and cyan, respectively (Figure S12). Of these 14 modules, the turquoise module had the most genes (1840), followed by the blue (648), brown (625), yellow (547), green (533), red (499), black (328), pink (175), magenta (168), purple (157), green-yellow (125), tan (84), salmon (Sun et al., 2021), and cyan (Liu et al., 2018) modules. Additionally, the correlations between these modules and the 18 differential metabolites related to phenylpropane biosynthesis were analyzed, and the correlation coefficients ranged from −0.92 to 0.99 (Figure 8). Among these 14 modules, the tan module had the highest correlation coefficient and showed significant negative associations (p < 0.1) with naringenin and p-coumaric acid, and the yellow and pink modules showed significant positive correlations (p < 0.1) with naringenin and p-coumaric acid, respectively.
Enlarge this image.FIGURE 8. WGCNA of phenylpropane biosynthesis genes. Each column corresponds to a specific metabolite, and each row corresponds to a module. The values within each cell indicate the correlation coefficient between the module and the metabolites, and the p-values are shown in parentheses. A high degree of correlation is indicated by dark red.
A total of five DEGs were selected to further investigate whether they are involved in the response to cold stress, and their expression levels were significantly induced by cold stress in cold-resistant and cold-sensitive V. sativa accessions, especially after cold treatment for 6 h (Figure 9a). Therefore, these DEGs were overexpressed in the S. cerevisiae yeast strain INVSc1, and an empty plasmid (pYES2) was used as a negative control. Then their growth characteristics were observed. As shown in Figure 9b, the growth rates of empty and recombinant plasmids were not significantly different under normal conditions. However, the recombinant plasmid cells showed higher growth rates compared to the empty plasmid cells under cold treatment conditions (−20°C for 36 h). Among these five recombinant plasmid cells, pYES2-c54604.graph_c0, pYES2-c20836.graph_c0, pYES2-c41936.graph_c0, and pYES2-c51037.graph_c0 were still able to grow in a 10−5 dilution concentration, but the pYES2 cell and the one remaining transgenic yeast (pYES2-c48739.graph_c0) cell barely grew in a 10−4 dilution. These results suggest that the cold resistance candidate genes obtained in the present study have the potential to improve cold resistance in V. sativa.
Enlarge this image.FIGURE 9. Overexpression of four DEGs in transgenic yeast enhances cold tolerance. (a) Heat map of the five transcription factors in six V. sativa accessions under cold treatment. Orange and blue indicate a relative increase and decrease in expression (log2FPKM), respectively. (b) Phenotypic growth assays of INVSc1 cells transformed with the pYES2 empty vector, pYES2-c54604.graph_c0, pYES2-c20836.graph_c0, pYES2- c41936.graph_c0, pYES2-c51037.graph_c0, and pYES2-c48739.graph_c0 were spotted on SC-Ura medium in 2 mL aliquots of ten-fold serially diluted (An et al., 2020) cultures under cold stress.
As the preferred method for understanding gene expression at different developmental stages in nonmodel plants, next-generation sequencing (NGS) has previously been used to reveal regulatory networks involved in abiotic stress in V. sativa (Min, Lin, et al., 2020; Min, Liu, et al., 2020). In the present study, a total of 54 cDNA libraries were constructed and 85,707 unigenes were assembled, with average and N50 lengths of 905.84 and 1835 bp, respectively; these values are higher than those reported in previous studies (Dong et al., 2017; Rui et al., 2017), but lower than those previously reported in V. sativa (Min, Lin, et al., 2020; Min, Liu, et al., 2020). Among these 85,707 unigenes, a total of 43,962 (51.29%) were annotated in five public databases, which is lower than those identified in previous studies (De la Rosa et al., 2020; Zhu et al., 2019). This result may be due to the differences in the samples used for sequencing and some tissue-specific genes. Moreover, the qRT-PCR results of 10 randomly selected unigenes were largely consistent with our transcriptome data (Figure S4), which suggests that the sequencing data in the present study were reliable and can be used for subsequent analysis.
Physiological changes inIt is widely known that ion leakage and MDA content can be used to reflect membrane integrity and the extent of oxidation injury and cellular damage during abiotic stress, respectively (Zhou et al., 2018). In both the cold-resistant or cold-sensitive V. sativa accessions, ion leakage and MDA content increased under cold stress, and their peak values were reached after 6 h of cold treatment (Figure 1b,c); this is consistent with the results of a previous report on the responses of Kandelia obovate and Zanthoxylum bungeanum to cold stress (Fei et al., 2021; Tian et al., 2021). However, the ion leakage and MDA content in the cold-sensitive accessions were higher than those in the cold-resistant accessions, which suggests that cold stress is more harmful to the cold-sensitive accessions. In the antioxidant enzyme system, the degradation of H2O2 by catalase (CAT), peroxidase (POD), and ascorbate peroxidase (APX), with changes in the activity of the relevant enzymes consistent with the H2O2 content (Zhao et al., 2023). And glutathione reductase (GR) is another enzyme belonging to the antioxidant system as well as the system of physiological GSH:GSSG ratio sustainment; this enzyme uses NADPH as an electron donor to reduce GSSG to GSH, which is a component of a nonenzymatic antioxidant system (Ciacka et al., 2022). In the present study, we found that the expression level of VsAPX1, VsGR1, VsCAT2, and VsCAT3 gradually increased during cold treatment, and VsCAT1 and four VsPOD genes were highly expressed in the early stage of cold treatment (Figure S13a). The upregulation of these genes can maintain the stability of plant cell membranes under abiotic stress by scavenging reactive oxygen species.
The plant's defense mechanisms are activated when it is subjected to abiotic stress; these mechanisms can include the synthesis of proline and soluble sugar, which can regulate the cell's osmotic potential, maintain turgor pressure, and partially equip plants to withstand dehydrative stresses. In the present study, we found that the proline content increased significantly during cold treatment (Figure 1e), and similar results were reported in M. sativa, Forsythia suspensa, and Z. bungeanum (Li, Shi, & Cushman, 2021; Tian et al., 2021; Zhou et al., 2018). Additionally, a total of 18 genes related to arginine and proline metabolism were obtained in all DEGs, and D1-Pyrroline-5-carboxylate synthetase (P5CS, c53263.graph_c0) and aldehyde dehydrogenase (ALDH, c50916.graph_c0) involved in proline synthesis were induced by cold treatment in V. sativa, while the expression level of proline dehydrogenase (c50343.graph_c0) was decreased during cold treatment (Figure S14), which is consistent with the change trend of proline content during cold treatment (Figure 1e). Furthermore, we investigated the expression level of several osmotic regulator genes in V. sativa during cold stress, such as galactinol synthase (GolS1), trehalase (TRE), and trehalose-5-phosphate phosphatase (TPP), which are involved in trehalose and raffinose synthesis (He et al., 2022; Qiu et al., 2022). The expression level of VsGolS1, VsTRE, and two VsTPP genes (c57302.graph_c1 and c60158.graph_c0) significantly increased after 48 h of cold treatment, and two other VsTPP genes (c33025.graph_c1 and c52674.graph_c0) were significantly induced after 6 h of cold treatment (Figure S13b). The upregulation of these genes promoted the accumulation of soluble sugar (Figure 1f) and ultimately maintained the osmotic pressure of cells during cold stress.
Plant hormones involved in the response ofPhytohormones are small-molecular-weight compounds that can regulate a wide range of metabolic processes at extremely low concentrations, and they act as a central regulator of the plant cold stress response (Eremina et al., 2016). In the present study, the KEGG enrichment analysis showed that the plant hormone signal transduction pathway was significantly enriched in the response of V. sativa to cold stress (Figures S9 and S10), which is consistent with previous results reported for M. sativa and Phyllostachys violascens (Hou et al., 2022; Zhou et al., 2018). Moreover, a total of 13 differential metabolites related to the biosynthesis of plant hormones and 72 DEGs related to plant hormone signal transduction were obtained in the present study (Figure S15, Table S7), and most DEGs were significantly correlated with these metabolites (Figure S16). Among the 13 differential metabolites, the contents of 11 metabolites in “Lan3” were significantly higher than those in “521” 48 h after cold treatment; this suggests that the synthesis of plant hormones in cold-sensitive accessions may be more severely affected by cold stress than in cold-resistant accessions.
A recent study reported that ABA and SA were significantly accumulated in Zea mays seedlings under low-temperature stress, and auxin and JA were decreased (Xiang et al., 2021). However, the opposite trend was also shown in other species, such as Zea mays and Vitis vinifera (Gao-Takai et al., 2019; Pal et al., 2020). A total of 16 DEGs related to ABA signal transduction were obtained in the present study, and the expression level of most genes was induced during cold treatment except for three genes (Figure S17a). The pyrabactin resistance (PYR)/PYR-like (PYL)/regulatory component of ABA receptor (RCAR) proteins that have been characterized as ABA receptors function as the core components in the ABA signaling pathway, and previous studies reported that PYL genes are positive regulators related to stress response in ABA signaling transduction, such as PYL4 and PYL9 (Nai et al., 2022; Ren et al., 2022). In our study, we found that three PYR/PYL genes (c25218.graph_c0, c34730.graph_c0, and c49758.graph_c0) were significantly induced by cold treatment (Figure S18a), which suggested that the ABA signaling pathway plays an important role in the response of V. sativa to cold stress.
In the past decade, a number of cold-responsive (COR) genes have been identified and characterized in both dicotyledonous and monocotyledonous plants, and COR gene expression is mediated by both ABA-dependent and ABA-independent pathways (Sun et al., 2009). Besides PYR/PYL genes, we investigated the expression level of other genes in the ABA-dependent pathway during cold stress, such as protein phosphatase 2C (PP2C), sucrose nonfermenting-1-related protein kinase 2 (SnRK2), and ABA responsive element binding protein (AREB)/ABRE binding factors (ABF) genes. As a result, the expression level of 12 PP2C genes significantly decreased during cold treatment, especially after 48 h of cold treatment (Figure S18b) and the expression level of five SnRK2 and four AREB/ABF genes significantly increased during cold treatment (Figure S18c,d). Moreover, we found that 11 dehydration-responsive element binding (DREB)/c-repeat binding transcription factors (CBF) in ABA-independent pathway were significantly induced during cold treatment (Figure S18e). These results suggested that some transcriptional regulatory pathway in V. sativa during cold treatment are ABA dependent, and others are ABA independent.
Additionally, an Arabidopsis transformant overexpressing the cytokinin (CK) biosynthetic gene isopentenyl transferase (DEX:IPT) showed enhanced cold tolerance, which was associated with elevated CK and SA levels in the shoots and auxin in the apices; moreover, a transformant overexpressing the CK degradation gene HvCKX2 (DEX:CKX) had weaker cold tolerance, accompanied by lowered levels of CKs and auxins (Prerostova et al., 2021). In the present study, trans-Zeatin content in both cold-resistant and cold-sensitive accessions was gradually increased during cold stress, but the increased ratio in the cold-sensitive accessions was significantly higher than that in the cold-resistant accessions (Figure S15). Moreover, a total of six DEGs related to CK signal transduction was obtained in the present study, and the expression level of most genes were induced during cold treatment except for one gene (Figure S17b). Cytokinin response regulators (RRs) are important components of the two component signal systems, which are involved in the regulation of plant growth and development and in the response to abiotic stress (Yang et al., 2016). Response regulators majorly consist of type-A and type-B RRs, and it has been shown that type-A RRs are extensively involved in abiotic stress responses, such as drought and cold stress responses (Zeng et al., 2021). In the present study, a total of five VsRRs were obtained, including RR2 (c57730.graph_c0, c50419.graph_c0, and c50419.graph_c1), RR8 (c55995.graph_c1), and B−type RRs (c41958.graph_c1). As shown in Figure S17b, except for c57730.graph_c0, the expression of all other genes is induced by cold treatment, but the cold resistance function of homologous genes of these genes has not been reported yet. These results show that the responses of plant hormones to cold stress vary among species.
Role of the phenylpropanoid pathway in the response ofThe phenylpropionate pathway is probably the most studied pathway in plant secondary metabolism, and it is activated under multiple abiotic stresses, such as drought, salt, heat, cold, heavy metals, and UV radiation (Sharma et al., 2019). In the present study, KEGG enrichment analysis showed that the biosynthesis of phenylpropanoids was the most significantly enriched pathway in all the DEGs and differential metabolites between cold-resistant and cold-sensitive V. sativa accessions (Figure 6); this is consistent with the results previously reported for V. sativa (Min, Liu, et al., 2020), and similar results have also been reported in other species, such as Actinidia chinensis, Fagopyrum tataricum, and Phyllostachys violascens (Hou et al., 2022; Jeon et al., 2018; Sun et al., 2021). In plants growing in challenging environments, the accumulation of phenylpropane compounds can protect them from stress damages as they possess antioxidant attributes and other protective properties (Wang et al., 2021). Among these 384 identified metabolites, a total of 18 differential metabolites related to the biosynthesis of phenylpropanoids were obtained; moreover, the contents of two key metabolites (p-coumaric acid and naringenin) in the flavonoid biosynthesis pathway in “Lan3” were significantly higher than those in “521” under the three cold treatment durations (Figure 10). This may be the main factor responsible for the difference in cold resistance between “Lan3” and “521.”
Enlarge this image.FIGURE 10. Biosynthesis of phenylpropanoids pathway in V. sativa. Heatmap of transcript levels for genes and contents for metabolites that are involved in the biosynthesis of phenylpropanoids.
Phenylalanine ammonia lyase (PAL), cinnamate-4-hydroxylase (C4H), and 4-coumarate-CoA ligase (4CL) are the major enzymes involved in the phenylpropanoid pathway, and PAL is an entry-point enzyme (Wang et al., 2021). Briefly, L-phenylalanine is catalyzed by PAL to yield trans-cinnamic acid and ammonia, which are subsequently transformed into 4-coumaroyl-CoA; then, the 4-coumaroyl-CoA is converted into numerous phenylpropanoid compounds, such as flavones, anthocyanins, and naringenin, via a series of catalytic reactions. Previous studies have reported that cold stress induces the expression of several structural genes in the phenylpropanoid pathway, including chalcone synthase (CHS) and 4CL; this leads to the accumulation of flavonoids and lignin to promote the adaptation of Arabidopsis, Eriobotrya japonica, and Malus domestica to low-temperature environments (An et al., 2020; Catala et al., 2011; Zhang, Yin, et al., 2020). In our study, we found that the expression levels of seven DEGs in the phenylpropanoid pathway were significantly increased after cold treatment, and the transcript levels of these seven genes in cold-resistant accessions were significantly higher than those in cold-sensitive accessions at three cold treatment durations (Figure 10).
Flavonoids are synthesized through the phenylpropanoid pathway, transforming phenylalanine into 4-coumaroyl-CoA, which finally enters the flavonoid biosynthesis pathway. In the present study, the KEGG enrichment analysis results indicate that the flavonoid biosynthesis pathway is significantly enriched during the cold treatment (Figure 6), and a total of 19 DEGs related to flavonoid biosynthesis were obtained (Figure S19). Flavanone 3-hydroxylase (F3H) is one of the key flavonoid biosynthesis genes, and the heterologous expression of DoF3H in Escherichia coli conferred it higher tolerance to salt and cold stresses (Si et al., 2023). We found two F3H genes in these 19 DEGs, and the expression level of these two F3H genes was induced by the cold treatment, especially after 48 h of cold treatment (Figure S19). These results suggest that the above genes in the phenylpropanoid pathway may enhance cold resistance in V. sativa.
TFs related to cold resistance inTranscription factors play critical roles in plant responses to abiotic stresses and directly regulate the expression levels of stress-responsive genes (Zhou et al., 2018). In the present study, a total of 241 TFs belonging to 27 families were detected. MYB and bHLH represented the most abundant categories, comprising 56 and 43 transcripts (Figure 4a), respectively, and the expression levels of these DEGs were significantly induced by cold stress in the cold-resistant and cold-sensitive V. sativa accessions (Figure 4b,c). In addition, we found that the expression level of c38862.graph_c0 was steady upregulation during cold treatment in V. sativa (Figure 4b), and the homologous gene of this gene in rice and Arabidopsis is MYB4. Previous studies have shown that overexpression of the rice Osmyb4 gene increases chilling and freezing tolerance of Arabidopsis plants (Vannini et al., 2004), and the Arabidopsis R2R3-MYB transcription factor AtMYB4 plays crucial roles in the regulation of phenylpropanoid metabolism and are involved in the response to abiotic stress (Li, Shi, Li, et al., 2021). However, the homologous gene MYB3 of c60205.graph_c0, another gene stably upregulated during cold treatment (Figure 4b), is a newly identified repressor in Arabidopsis phenylpropanoid biosynthesis (Zhou et al., 2017). A previous study reported that MtMYB3 binds directly to MYB cis-elements in the MtCBF4 promoter, leading to the inhibition of MtCBF4 expression and decreasing the freezing tolerance of M. truncatula plants (Zhang et al., 2016). Therefore, we speculate that these MYB transcription factors may also regulate the response of V. sativa to cold stress through the phenylpropane biosynthesis pathway.
The AP2/ERF TF family is one of the most important cold stress-related TF families in plants, and along with other TF families, such as WRKY, bHLH, bZIP, MYB, and NAC, it enhances cold stress tolerance (Ritonga et al., 2021). For example, a recent study reported that the AP2/ERF family protein MdABI4 can interact with a key regulatory protein, MdICE1, which is involved in cold stress response to positively regulate cold tolerance in apples (An et al., 2022). A total of five DEGs were selected to further investigate the involvement of the identified DEGs in cold stress; moreover, two ERF TFs (c20836.graph_c0 and c54604.graph_c0) were overexpressed in yeast, and it was indicated that their biological functions is to confer cold tolerance (Figure 9b). These findings provide a genetic basis and abundant gene resources for further study of the molecular mechanisms of the cold stress response, laying a foundation for cold-resistant V. sativa breeding research.
CONCLUSIONVicia sativa is cultivated worldwide due to its ability to grow in extreme abiotic stress conditions, and it is an important protein source for livestock and humans. However, cold stress is one of the most important limitations to its utilization at high altitudes. In the current study, during cold treatment, ion leakage and the contents of H2O2 and MDA in the cold-sensitive accessions were significantly higher than in the cold-resistant accessions, which suggests that cold stress is more harmful to cold-sensitive accessions. In order to determine the molecular regulatory mechanism of V. sativa during cold stress, the shoot tissues of three cold-resistant and three cold-sensitive accessions were collected after cold (4°C) treatment for 0, 6, and 48 h, followed by integrative analyses of their transcriptomes and metabolomes. KEGG enrichment showed that the biosynthesis of phenylpropanoids was the most enriched pathway in the cold stress response of V. sativa. A total of 18 differential metabolites were obtained in the biosynthesis of phenylpropanoids pathway, and the contents of p-coumaric acid and naringenin in “Lan3” were significantly higher than those in “521” under the three cold treatment durations. Moreover, the expression patterns of seven DEGs in the phenylpropanoid pathway were consistent with those of the above two metabolites. In conclusion, the results obtained in this study showed that Ca2+, ROS, and hormonal signaling couple with TFs play important roles in regulating the phenylpropane and flavonoids biosynthetic pathways to respond to cold stress (Figure 11). These results reveal that the phenylpropanoid pathway may be involved in the cold stress response of V. sativa.
Enlarge this image.FIGURE 11. Models describing the signaling pathways involved in the acquisition of cold tolerance in V. sativa.
Qiang Zhou: Investigation, Data curation, Writing—original draft. Yue Cui: Data curation, Formal analysis. Shuwei Dong: Data curation. Dong Luo: Formal analysis. Longfa Fang: Formal analysis. Zunji Shi: Formal analysis. Wenxian Liu: Formal analysis. Zengyu Wang: Resources. Zhibiao Nan: Resources. Zhipeng Liu: Conceptualization, Writing—review and editing, Funding acquisition, Supervision.
ACKNOWLEDGMENTSThis research was supported by the National Natural Science Foundation of China (32071862 and 32201444), the Gansu Provincial Science and Technology Major Projects (22ZD6NA007), the Qinghai Province Handsome Scientist Responsibility System Project (2023-NK-147-1), and the China Postdoctoral Science Foundation (2020M683609).
CONFLICT OF INTEREST STATEMENTThe authors have stated explicitly that there are no conflicts of interest in connection with this article.
DATA AVAILABILITY STATEMENTAll sequencing reads are available in NCBI SRA: SRR20083627–SRR20083644 (
© 2023. This work is published under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
; Cui, Yue 2 ; Dong, Shuwei 2 ; Luo, Dong 2 ; Fang, Longfa 2 ; Shi, Zunji 2 ; Liu, Wenxian 2 ; Wang, Zengyu 3 ; Zhibiao Nan 2 ; Liu, Zhipeng 2
1 State Key Laboratory of Herbage Improvement and Grassland Agro-ecosystems, Key Laboratory of Grassland Livestock Industry Innovation, Ministry of Agriculture and Rural Affairs, Engineering Research Center of Grassland Industry, Ministry of Education, College of Pastoral Agriculture Science and Technology, Lanzhou University, Lanzhou, China; College of Ecology, Lanzhou University, Lanzhou, China
2 State Key Laboratory of Herbage Improvement and Grassland Agro-ecosystems, Key Laboratory of Grassland Livestock Industry Innovation, Ministry of Agriculture and Rural Affairs, Engineering Research Center of Grassland Industry, Ministry of Education, College of Pastoral Agriculture Science and Technology, Lanzhou University, Lanzhou, China
3 Grassland Agri-Husbandry Research Center, College of Grassland Science, Qingdao Agricultural University, Qingdao, China