Comparative transcriptome analysis of two varieties of common bean (Phaseolus vulgaris L.) to identify candidate drought resistance genes

Abstract

In this work, we studied the physiological, biochemical and molecular mechanisms underlying the response to polyethylene glycol (PEG) in a drought-resistant variety ‘HZYTJ’ and the sensitive variety ‘WZWJD’ of common beans. The two bean varieties were subjected to a 10% PEG simulated drought stress treatment, and the resultant physiological and biochemical characteristics of the leaves and roots were compared. The superoxide dismutase (SOD), peroxidase (POD), catalase (CAT) and other indicators were higher in ‘HZYTJ’ than in ‘WZWJD’; additionally, the ‘HZYTJ’ exhibited slower changes in response to stress. The proline content in ‘HZYTJ’ changed slowly compared to that in ‘WZWJD’. The transcriptome analysis of the two varieties showed that 1,361 and 1,447 genes were differentially expressed in ‘HZYTJ’ and ‘WZWJD’, respectively. DEGs were also analyzed via the KEGG pathway database, and 10 significant metabolic pathways were found to be involved in ‘HZYTJ’, while nine were found in ‘WZWJD’. Using the GO function database, we found that 111 significant GO terms were involved in ‘HZYTJ’, whereas 125 were involved in ‘WZWJD’. Through analysis of transcription factors, 487 DEGs were identified, of which bHLH transcription factors were the most abundant, followed by ERF, MYB-related, B3 and other transcription factor families. The selected and screened genes with DEGs were validated by qRT–PCR. These results provide a basis for further research into the drought resistance of the common bean.

Keywords:

• Phaseolus vulgaris L.
• drought stress
• transcriptome
• drought resistance

Introduction

The common bean (Phaseolus vulgaris L.) belongs to Leguminosae, Papilionoideae, and Phaseolus [1]. The environmental conditions necessary for the successful cultivation of the common bean are complex; specifically, its cultivation is easily impacted by a variety of biotic and abiotic factors, and the impact of drought is particularly serious. Drought can directly affect the emergence rate of common beans and the growth of seedlings, thereby affecting the yield of common beans. Within the context of anthropogenic climate change, drought spells have become more severe in some common-bean-producing areas [2]. The intrinsic responses of the common bean to drought have shown that the magnitude of drought resistance varies depending on the intensity of the drought at hand. Thus, drought resistance has garnered the attention of many researchers across the globe, and a series of studies have been conducted in common bean, to name a few [3–7].

The common bean has acquired a complex set of physiological mechanisms linked to stress resistance. The acquisition of such capabilities can be attributed to the complex environment in which the common bean grows. During periods of drought, plants exhibit physiological and biochemical adaptive responses, such as a decline in the photosynthetic rate, the accumulation of ABA and osmotic adjustment substances, and the operation of antioxidant systems; overall, these changes help the plant to effectively cope with the drought stress [5, 8]. Many studies have pointed to the uniqueness of the adaptive responses of the common bean (i.e. and its varieties) under drought stress [6–12].

Transcriptome sequencing (RNA-seq) is a high-throughput sequencing technology that detects the gene transcription of an organism under specific environmental conditions. Currently, RNA-seq technology is widely used not only in plants but also in the mining of plant drought resistance genes. Transcriptome sequencing was used to compare and analyze the gene expression levels in many plant species under drought stress to identify drought resistance genes [13–16].

In this study, we performed RNA-seq on the two common bean varieties in three different states under PEG-simulated drought stress. Using high-throughput sequencing combined with bioinformatic analyses, we studied the differences among the two varieties at the transcriptional level. By determining and analyzing the physiological responses of the leaves and roots of two common bean varieties with different drough tolerance, this study will provide a theoretical basis for improving the quality and yield of common beans.

Materials and methods

Common bean varieties

The varieties of the common bean used in this study were the ‘Hong zi ya ta jia’ (‘HZYTJ’) and the ‘Wang zhong wang jia dou’ (‘WZWJD’) [17]. The ‘HZYTJ’ is a variety with relatively strong drought resistance. This variety is characterized by high yield potential and strong disease resistance. The ‘WZWJD’ variety has relatively weak drought resistance. The characteristic is that it is an early maturing variety with relatively low yield. The study was conducted in the artificial climate room of the intelligent multispan greenhouse of Shanxi Agricultural University, Taigu County, Shanxi Province (altitude 790 m, longitude 112°37’26", latitude 37°29’17").

Treatment

The seeds (i.e. with full grains, 100 seeds) of the two common bean varieties were selected and soaked in warm water. Once the true leaves of the common bean varieties were fully expanded, six seedlings for each treatment were removed and planted in 35-L plastic dark water tanks with 1/2 Hoagland solution [18, 19]. Hydroponic seedlings were treated with 10% PEG-6000 for 6 h to simulate drought stress, with three biological replicates (Supplemental Figure S1). Samples (leaves and roots) were then collected at three time points, specifically, when untreated (0 h), when wilting after treatment (6 h after treatment), and after transfer to 1/2 Hoagland solution rehydration treatment (12 h after transfer) to determine if the stress would be reflected in various physiological indicators (superoxide dismutase (SOD), peroxidase (POD) and catalase (CAT) activity and malondialdehyde (MDA) content). Three biological replicates were run for each physiological indicator.

Sample preparation

The samples were homogenized in liquid nitrogen using a precooled mortar and pestle at 4 °C. The homogenized samples in PBS (pH 6.0) were then centrifuged it at 4,000 g for 15 min. The homogenate was filtered through two layers of filter paper. The supernatant was immediately used for enzyme activity analysis.

Measurement of physiological indicators

The antioxidant SOD, POD and CAT enzyme activities were determined according to the methods described [20–22]. The main reagents for enzyme activities determination were purchased from the Solarbio Company (Beijing, China). We quantified MDA content by the thiobarbituric acid assay with slight modifications [23]. In brief, we ground about 500 mg of sample in liquid nitrogen and homogenized it in 0.1% (w/v) trichloroacetic acid (TCA) solution (0.5 mL). For the reaction mixture, 1 mL supernatant was added to 20% TCA (2 mL) containing 0.6% thiobarbituric acid (TBA) in a clean glass tube. We cooled it immediately thereafter, then centrifuged it at 3,000 g for 10 min at 4 °C. We measured the absorbance of the clarified supernatant at 600, 532 and 450 nm. The soluble protein content in leaves and roots was determined by Coomassie brilliant blue [24]. The samples (leaves and roots) (0.1 g) were ground into a homogenate using 5 mL of 50 mmol L⁻¹ phosphate buffer (pH 7.8). The supernatant (0.1 mL) was transferred to a test tube and 4 mL Coomassie brilliant blue (0.01%) and 1 mL distilled water were added. After standing for 3 min, absorbance was read at 595 nm in a UV spectrophotometer. The soluble sugar content was determined by anthrone colorimetry [25]. The reaction mixture containing 1 mL supernatant and 5 mL anthrone (100 mg anthrone + 100 mL 80% H₂SO₄) was heated at 100◦C for 10 min, and then absorbance was read at 620 nm. The proline content was determined according to the method described by Bates et al. [26]. In brief, fresh leaves or roots (0.5 g) were ground with 5 mL of 3% sulfosalicylic acid. The extract was filtered and acid ninhydrin solution (2 mL) was added to the filtrate. The absorbance of the solution was read at 520 nm. The content was determined using a proline standard curve.

RNA extraction

Samples from the seedling leaves representing the drought stress of ‘HZYTJ’, namely, 0 h (CK1), 6 h (T1) and 12 h (T3); and those representing the drought stress of ‘WZWJD’, namely, 0 h (CK2), 6 h (T2) and 12 h (T4) were immediately frozen in liquid nitrogen and stored at −80 °C. Plant RNA extraction kit (Tiangen Company, Beijing, China) was used to extract total RNA from the leaves of the two common bean varieties.

RNA library construction and sequencing methods

The concentration, purity and integrity of the extracted RNA were then determined by a spectrophotometer and an Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA) before freeze storage. The mRNA-seq library was constructed using the Illumina’s TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA). A total amount of 1 μg RNA per sample was used to prepare the library with insert sizes of 350 bp and sequenced on a NovaSeq 6000 (Illumina, San Diego, CA, USA). The samples were sequenced on NovaSeq 6000 at Bioeditas Technology Corporation (Shaanxi, China). The original image files were converted into raw data, and then the raw reads obtained through sequencing were filtered to remove the following reads: reads with connector contamination and reads with an N content greater than 5%. Overall, these steps were taken to obtain high-quality reads that would be subsequently spliced. The transcript sequence obtained through splicing was classified and characterized according to different genome annotation information. HISAT2 tools were used to map the reference genome (ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/499/845/GCF_000499845.1_PhaVulg1_0/GCF_000499845.1_PhaVulg1_0_genomic.gff.gz).

Differentially expressed gene (DEG) analysis

To identify the differentially expressed genes (DEGs) between the compared samples, FDR <0.001 and |log2Ratio|≥2 found by DESeq was set as the threshold for significant differential expression, at a p value <0.05 [27, 28]. Differential expression analysis of the two cultivars was performed using the DEGseq R package. These genes were then analysed using two major databases, namely, Gene Ontology (GO) enrichment significance analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathway database (http://www.genome.jp/kegg/) [29, 30]. Finally, the differential genes were enriched and analyzed. The GOseq R packages and KOBAS software were used to test the statistical enrichment of DEGs [31].

Gene validation and expression analysis by qRT–PCR

Quantitative real-time reverse transcription PCR (qRT–PCR) was performed to confirm the results of the DEG analysis. A Hiscript II Q RT SuperMix (Vazyme Company, Nanjing, China) was used for qPCR reverse transcription for total RNA reverse transcription to obtain cDNA as a template. The primers designed by Primer Premier 5 used were as shown in Supplemental Table S1. The qRT-PCR was performed on the Roche Light Cycler® 480 system (Roche Diagnostics, CA, USA). The nine genes selected were tested for functional verification, primers were designed using Primers software, the Actin gene was used as the internal reference gene, and the 2 × SYBR® Green dye method was used for qRT-PCR analysis [32]. The qRT–PCR reactions were repeated three times. The 2^-ΔΔCT method was finally used to calculate the relative expression of genes.

Data analysis

We performed statistical analyses by one-way analysis of variance (ANOVA) in SPSS (version 25.0) with Duncan’s test. The different letters indicate significant differences at p< 0.05 (Student’s t-test).

Results

Effect of PEG-simulated drought stress and rehydration on the activity of common bean antioxidant enzymes

As shown in Figure 1(A), the SOD activity (U·mg⁻¹ FW (fresh weight) in the leaves and roots of the two common bean varieties during the 10% PEG simulated drought stress treatment and rehydration first declined and then increased as the treatment progressed. Additionally, the change in the SOD activity of the common bean leaf was greater than that of the root. Specifically, the SOD activity in the leaves of ‘HZYTJ’ decreased by 43.72% 6 h into the treatment, to 2.69 U·mg⁻¹, and increased 12 h into the treatment to 4.84 U·mg⁻¹. The SOD activity after 6 h of treatment was significantly lower than that at 0 h and 12 h (p < 0.05). The SOD activity in the leaves of ‘WZWJD’ was lower than that in the leaves of ‘HZYTJ’. After 6 h, the SOD activity in the leaves of ‘WZWJD’ decreased by 49.47%, to 2.39 U·mg⁻¹, and after 12 h, it increased to 4.66 U·mg⁻¹. The SOD activity at 6 h was significantly lower than that at 0 h and 12 h (p < 0.05).

Figure 1. Changes in SOD (A), POD (B) and CAT (C) activity (U·mg⁻¹ fresh weight) in leaves and roots of common bean varieties ‘WZWJD’ and ‘HZYTJ’ during 10% PEG-simulated drought stress and rewatering. Samples (leaves and roots) were collected at three time points, specifically, when untreated (0 h), when wilting after treatment (6 h after treatment), and after transfer to 1/2 Hoagland solution rehydration treatment (12 h after transfer). Error bars are standard errors of the mean from three biological replicates. The different letters indicate significant differences at p < 0.05 (student’s t-test).

The SOD activity in the roots of ‘HZYTJ’ fluctuated but, there was no significant difference across the three treatment times. The SOD activity in the roots of ‘WZWJD’ decreased by 31.12% at 6 h, to 2.70 U·mg⁻¹, and rose to 3.36 U·mg⁻¹ at 12 h. The SOD activity at 6 h was significantly lower than that at 0 h and 12 h (p < 0.05).

As shown in Figure 1(B), as the treatment progressed, the POD activity in common bean leaves and roots first decreased and then increased; specifically, the POD activity in common bean leaves was significantly lower than that in roots. The change in the POD activity in the leaves was slower than that observed in the roots. In the leaves, the POD activity of ‘HZYTJ’ at 6 h was significantly lower than the POD activity at 0 h and 12 h. At 6 h, the POD activity of ‘HZYTJ’ was reduced by 19.90%, to 4,400 U·mg⁻¹, and at 12 h, the POD activity increased to 5,427 U·mg⁻¹. The POD activity in the ‘WZWJD’ leaves was significantly different across the three treatment times. The POD activity at 12 h was the highest, at 2,920 U·mg⁻¹, which was significantly higher than the POD activity at 0 h and the POD at 6 h (p < 0.05). The lowest POD activity in the ‘WZWJD’ leaves at 6 h was 1,970 U·mg⁻¹.

In the roots, there were significant differences in the POD activity of the two varieties at the three time points. The POD activity at 0 h was significantly higher than that at 12 h and 6 h, and the POD activity at 12 h was significantly higher than that at 6 h. The POD activity of ‘HZYTJ’ was 14,800 U·mg⁻¹ at 0 h, and the 12 h POD activity increased by 24.21% compared to the 6 h POD activity. The POD activity of the ‘WZWJD’ root at 0 h was 13,753 U·mg⁻¹, and at 12 h it increased by 31.57% compared with the POD activity at 6 h (p < 0.05).

As shown in Figure 1(C), the CAT activity in the common bean leaves was higher than that in roots. In the leaves, the CAT activities of the two common bean varieties across the three treatment times were significantly different from each other. The CAT activity of ‘HZYTJ’ at 6 h was 4.24 Umg·g⁻¹, which was 46.53% lower than the CAT activity at 0 h, which increased to 9.05 Umg·g⁻¹ at 12 h. The CAT activity of ‘WZWJD’ at 6 h was 2.80 Umg·g⁻¹, which was 36.51% lower than the CAT activity at 0 h, which increased to 6.57 Umg·g⁻¹ at 12 h (p < 0.05).

In the roots, there were significant differences in CAT activity across the three time points of ‘HZYTJ’. The CAT activity at 0 h (i.e. 3.96 Umg·g⁻¹) was significantly higher than the CAT activity at 6 h and 12 h (p < 0.05); the minimum CAT activity at 6 h was 1.78 Umg·g⁻¹. There was no significant difference in the CAT activity of ‘WZWJD’ across the three time points; the highest CAT activity at 12 h was 2.17 Umg·g⁻¹, and the lowest at 6 h was 1.67 Umg·g⁻¹ (p > 0.05).

PEG simulated drought stress and the effect of rehydration on osmotic regulators of common bean

As shown in Figure 2(A), the soluble protein content in common bean leaves was significantly higher than that of the roots (p > 0.05). Additionally, in the leaves, the soluble protein content of ‘HZYTJ’ was significantly different across the three time points. As the treatment progressed, the soluble protein content gradually increased; the content at 12 h was 0.30 mg·g⁻¹, and the minimum value at 0 h was 0.24 mg·g⁻¹. There was no significant difference in the soluble protein content in the leaves of ‘WZWJD’ at 0 h and 12 h, and at 6 h the content was the lowest at 0.16 mg·g⁻¹.

Figure 2. Changes in soluble protein (A) and sugar (B) in leaves and roots of common bean varieties ‘WZWJD’ and ‘HZYTJ’ during 10% PEG-simulated drought stress and rewatering. Samples (leaves and roots) were collected at three time points, specifically, when untreated (0 h), when wilting after treatment (6 h after treatment), and after transfer to 1/2 Hoagland solution rehydration treatment (12 h after transfer). Error bars are standard errors of the mean from three biological replicates. The different letters indicate significant differences at p < 0.05 (student’s t-test).

In the roots, there were significant differences in the soluble protein content of ‘HZYTJ’; the highest content at 12 h was 0.13 mg·g⁻¹, and the lowest content at 0 h was 0.08 mg·g⁻¹. There was no significant difference in the soluble protein content in the roots of ‘WZWJD’ at 0 h and 6 h (p > 0.05), and it reached its highest at 12 h, which was 0.16 mg·g⁻¹.

As shown in Figure 2(B), the soluble sugar content in common bean leaves was significantly lower than that of the roots. In the leaves, the soluble sugar content of the two common bean varieties was significantly different. The soluble sugar content in ‘HZYTJ’ was 0.055 mg·g⁻¹ at 0 h, and the minimum was 0.045 mg·g⁻¹ at 6 h. The content of soluble sugar in the ‘WZWJD’ leaves was different from that of ‘HZYTJ’. The soluble sugar content in ‘WZWJD’ leaves at 12 h (0.049 mg·g⁻¹) was significantly higher than that at 6 h (0.033 mg·g⁻¹) and 0 h (0.031 mg·g⁻¹).

In ‘HZYTJ’ roots, there was a significant difference in the soluble sugar content. The content at 6 h (0.167 mg·g⁻¹) was significantly higher than the content at 0 h and 12 h (p < 0.05), and there was non-significant difference between the content of soluble sugar at 0 h and 12 h (p > 0.05). There were significant differences in the soluble sugar content in ‘WZWJD’ across the three time points. The highest soluble sugar content was 0.119 mg·g⁻¹ at 6 h, and the lowest was 0.078 mg·g⁻¹ at 12 h.

As shown in Figure 3(A), the proline content in common bean leaves was significantly higher than that in roots. In the leaves of ‘HZYTJ’, the proline content at 6 h (0.0052 mg·g⁻¹) was significantly higher than the proline content at 0 h and 12 h (p < 0.05), and the proline content at 0 h was significantly lower than that at 12 h (p < 0.05). The proline content in ‘HZYTJ’ leaves at 0 h was 0.0019 mg·g⁻¹. The proline content of ‘WZWJD’ at 6 h (0.0045 mg·g⁻¹) was significantly higher than that at 0 h and 12 h. There was no significant difference in the proline content between 0 h and 12 h in the leaves.

Figure 3. Changes in the proline contents (mg·g⁻¹ fresh weight) (A) and MDA contents (Nmol·mg⁻¹ fresh weight) (B) in leaves and roots of common bean varieties ‘WZWJD’ and ‘HZYTJ’ during 10% PEG-simulated drought stress and rewatering. Samples (leaves and roots) were collected at three time points, specifically, when untreated (0 h), when wilting after treatment (6 h after treatment), and after transfer to 1/2 Hoagland solution rehydration treatment (12 h after transfer). Error bars are standard errors of the mean from three biological replicates. The different letters indicate significant differences at p < 0.05 (student’s t-test).

In the roots of ‘HZYTJ’, the proline content at 0 h was significantly higher that at 6 h and 12 h and was 0.0017 mg·g⁻¹. There was no significant difference between the proline content at 6 h and 12 h. The proline content in the roots of ‘WZWJD’ at 12 h (0.00168 mg·g⁻¹) was significantly higher than that at 0 h and 6 h, and the proline content at 0 h (0.0006 mg·g⁻¹) was significantly lower than that at 6 h (0.00114 mg·g⁻¹).

The effect of PEG-simulated drought stress and rehydration on the malondialdehyde content of common bean

As shown in Figure 3(B), under 10% PEG-simulated drought stress, the MDA content of roots was significantly higher than that of leaves. In the leaves, the MDA content of ‘HZYTJ’ was significantly different at 0 h and 6 h. The MDA content at 6 h was 3.87 nmol·mg⁻¹ fresh weight, and the minimum content at 0 h was 2.58 nmol·mg⁻¹. The MDA content of ‘WZWJD’ was significantly different among the three time points. The MDA content at 6 h was the highest at 4.79 nmol·mg⁻¹, and the MDA content at 12 h was the lowest at 1.91 nmol·mg⁻¹ (p < 0.05).

In the roots, the MDA content of ‘HZYTJ’ at 6 h was significantly higher than the MDA content at 0 h and 12 h, which was 11.10 nmol·mg⁻¹, and the MDA content at 0 h (13.35 nmol·mg⁻¹) was significantly lower than that at 6 h and 12 h (p < 0.05). The MDA content of ‘WZWJD’ at 6 h was significantly higher than that at 0 h and 12 h, which was 11.91 nmol·mg⁻¹, and the MDA content at 12 h was significantly lower than that at 0 h and 6 h, which was 7.84 nmol·mg⁻¹ (p < 0.05).

Transcriptome analysis of common bean under drought stress

Data quality assessment and comparison analysis

The RNA library was constructed using three sampling time points from each of the two varieties (‘HZYTJ’ and ‘WZWJD‘). For transcriptome sequencing, the raw reads produced by each library were filtered. We found that the proportion of clean reads of the six samples was higher than 89%, which will allow for a follow-up assembly. The Q30 values of the samples were all greater than 94%. The above-mentioned results show that the sequencing accuracy was high, which means the sequences can be used for subsequent analyses.

Sample correlation analysis and cluster analysis

Correlation and cluster pattern analyses were performed on 18 samples after three repetitions of six samples. As shown in Figure 4, samples CK13 and T43 were abnormal in their clustering. Therefore, to improve the accuracy of the results, this study removed samples CK13 and T43 prior to the samples being processed for subsequent analysis.

Figure 4. Correlation of RNA-seq data among different samples of common bean varieties ‘WZWJD’ and ‘HZYTJ’. Heatmap (A) of correlation coefficient values across samples based on RNA-seq FPKM. Cluster tree (B) based on the distances of expressed genes. CK11, CK12 and CK13 are the three biological replicates of ‘HZYTJ’ sampled at 0 h; T11, T12 and T13 at 6 h; T31, T32 and T33 at 12 h. CK21, CK22 and CK23 are the three biological replicates of ‘WZWJD’ sampled at 0 h; T21, T22 and T23 at 6 h; T41, T42 and T43 at 12 h. Samples (leaves and roots) were collected at three time points, specifically, when untreated (0 h), when wilting after treatment (6 h after treatment), and after transfer to 1/2 Hoagland solution rehydration treatment (12 h after transfer).

Differentially expressed gene (DEG) analysis

In this experiment, 24,241 genes in 16 samples were clearly annotated or speculatively analysed. After screening and analysis, 1,361 differentially expressed genes were obtained from ‘HZYTJ’ after PEG stress, of which 726 were upregulated genes and 635 were downregulated genes. In ‘WZWJD’, 1,447 differentially expressed genes were obtained, and among them, 506 and 941 genes were upregulated and downregulated, respectively. Overall, 545 genes were shared by the two varieties.

Functionally related genes usually have similar expression patterns under the same conditions, as shown in Figure 5. Through hierarchical cluster analysis of the differentially expressed genes of two bean varieties at different times, a heatmap of the differential expression of the two varieties was generated. According to the expression level of the differentially expressed genes in each sample, we took the logarithm to the base of 2 and then used K-means clustering to determine the clustering of the genes. The differential genes were classified into six different clusters according to their expression patterns (Figure 6). Based on the results, class 2, class 3 and class 6 were analysed further. After the DEG analyses of the two common bean varieties, 599 DEGs between ‘HZYTJ’ with PEG and without PEG and 544 DEGs between ‘WZWJD’ with PEG and without PEG, respectively.

Figure 5. Heatmap analysis of DEGs. DEGs analysis of the ‘HZYTJ’ transcriptome (A) and ‘WZWJD’ transcriptome (B). The DEGs gene expression level is shown by the logarithmic normalized RNA-seq FPKM value.

Figure 6. Six clusters of coexpressed genes and the K-means clustering pattern analysis of the DEGs. The gene expression level is shown by the logarithmic normalized RNA-seq FPKM value. The grey lines in the figure show the expression patterns of genes in each cluster, and the dark purple line represents the average expression of all genes in the cluster in the sample.

DEGs functional annotation

The selected DEGs were compared with the databases (GO and KEGG), and the functions of unknown genes were inferred according to genes whose functions are known (i.e. genes in the database). The GO (Gene Ontology) database divides the function of genes into three categories, namely, biological process, molecular function, and cellular component. GO functional analysis was performed on the two common bean varieties in this study, and the results are shown in Supplemental Figure S2. The annotated genes of ‘HZYTJ’ (T13-C1) and ‘WZWJD’ (T24-C2) were subdivided into 51 and 49 GO functional classes, respectively. We found that the two common bean varieties had most of their DEGs involved in biological processes.

According to the calculation of the number of DEGs annotated into each GO entry, the hypergeometric test method was used to find the GO entries that were significantly enriched with DEGs. In ‘HZYTJ’, 111 significant GO classification terms were involved; among them, there were 67 metabolic pathways in the biological processes category, accounting for 60.36% of the significant GO classification terms, and a total of 4,891 DEGs; the cell components category had 34 metabolic pathways, accounting for 30.63% of the significant GO classification terms, and a total of 3,834 DEGs; and the molecular functions category had 10 metabolic pathways, accounting for 9.01% of the significant GO classification terms, and a total of 1,431 DEGs. In the biological processes category, there were many different genes linked to responses to abiotic stimuli, oxalate metabolism, signal transduction, response to chemical substances and oxygenated compounds. In the cell components category, there were many different genes, such as those linked to membrane components, membrane inherent components, extracellular regions, and plasma membrane components. In the molecular functions category, there were genes linked to oxidoreductase activity, DNA binding transcription factor activity, transcription regulation activity, transmembrane transport protein activity and transportation activity.

In ‘WZWJD’, a total of 125 significant GO classification terms were involved, of which there were 61 metabolic pathways in biological processes, accounting for 48.8% of the significant GO classification terms, and a total of 2,836 DEGs; there were 52 metabolic pathways in cell components, accounting for 41.6% of the significant GO classification terms, and a total of 8,031 DEGs; and molecular functions had 12 metabolic pathways, accounting for 9.6% of the significant GO classification terms, and a total of 2,025 DEGs. In the biological processes category, there were many different genes linked to responses to non-biological stimuli, small molecule catabolism, oxalic acid metabolism, organic acid metabolism, and carboxylic acid metabolism. In the cell components category, there were many DEGs, such as those linked to biofilm, membrane components, membrane inherent components, membrane components, and cytoplasmic parts. In the molecular functions category, there were many different genes associated with calcium ion binding, catalytic activity, oxidoreductase activity, chlorophyll binding and tetrapyrrole conjugates.

Using the hypergeometric test method, the path with q ≤ 0.05 was defined as the path of significant enrichment in differentially expressed genes, and significant enrichment analysis of metabolic pathways was performed on the two samples under PEG stress for significantly differentially expressed genes (Tables 1 and 2). Ten significant metabolic pathways were found in ‘HZYTJ’, namely, plant hormone signal transduction, MAPK signalling pathway in plants, peroxisomes, fatty acid metabolism, valine, leucine acid and isoleucine degradation, keratin, amber and wax biosynthesis, glyoxylic acid and dicarboxylic acid metabolism, fatty acid degradation, photosynthesis-antenna protein, and α-linolenic acid metabolism. There were 215 genes in these pathways. Among them, there were 70 differentially expressed genes in the plant hormone signal transduction pathway.

Table 1. Analysis of the main metabolic pathways of ‘HZYTJ’ under drought stress. The expression difference multiple log2 (fold change) ≥ 2, false discovery rate (FDR) <0.01, and p value <0.05 were used as screening criteria, and the DEGs genes were selected.

Name	p Value	Corrected p value	DEGs
Plant hormone signal transduction	8.64E-12	1.17E-09	70
Photosynthesis-antenna proteins	9.89E-09	6.73E-07	12
Valine, leucine and isoleucine degradation	0.000346	0.009419	14
Cutin, suberin and wax biosynthesis	0.000533	0.010358	14
MAPK signaling pathway-plant	0.000509	0.010358	32
Glyoxylate and dicarboxylate metabolism	0.001404	0.023872	14
Fatty acid degradation	0.001764	0.025469	13
Peroxisome	0.001873	0.025469	19
alpha-Linolenic acid metabolism	0.003194	0.038546	12
Fatty acid metabolism	0.003401	0.038546	15

Table 2. Analysis of the main metabolic pathways of ‘WZWJD’ under drought stress. The expression difference multiple log2 (fold change) ≥ 2, false discovery rate (FDR) <0.01, and p value <0.05 were used as screening criteria, and the DEGs genes were selected.

Name	p Value	Corrected p value	DEGs
Photosynthesis-antenna proteins	2.65E-12	3.60E-10	15
Glyoxylate and dicarboxylate metabolism	1.67E-08	7.58E-07	23
Photosynthesis	4.08E-07	1.39E-05	15
Porphyrin and chlorophyll metabolism	1.09E-05	0.000298	17
Carotenoid biosynthesis	0.000142	0.002755	13
Carbon fixation in photosynthetic organisms	0.00126	0.019044	12
Plant-pathogen interaction	0.001255	0.019044	43
Carbon metabolism	0.001703	0.02203	39
Isoquinoline alkaloid biosynthesis	0.001782	0.02203	8

There were nine significant metabolic pathways in ‘WZWJD’, namely, plant pathogen interaction, carbon metabolism, glyoxylic acid and dicarboxylic acid metabolism, porphyrin and chlorophyll metabolism, photosynthesis, photosynthesis-antenna biosynthesis of protein and carotenoids, carbon fixation in photosynthetic organisms and biosynthesis of isoquinoline alkaloids. There were 185 genes in these pathways, most of which were linked to plant pathogens (i.e. 43) and carbon metabolism (i.e. 39). Both varieties exhibited metabolic pathways linked to glyoxylic acid and dicarboxylic acid and photosynthesis-antenna protein pathways.

Transcription factor analysis

In the analysis of transcription factors, 58 transcription factor families were identified. Among them, 487 differentially expressed genes were involved. Furthermore, bHLH transcription factors were the greatest (i.e. 35), followed by ERF, MYB-related, B3 and C2H2. Among other transcription factor families, there were more than 20 differentially expressed genes. The top 18 transcription factor families had more than 10 differentially expressed genes. There were nine transcription factor families with only one DEG. The detailed results are shown in Supplemental Table S2.

Gene validation and expression analysis by qRT–PCR

Nine genes exhibiting significant differential expression were selected and screened by qRT–PCR. The nine genes are XM_007143360.1 (PHAVU_007G071000g), EIXI-2 (XM_007143566.1) (PHAVU_007G087800g), CALM (XM_007144763.1) (PHAVU_007G187200g), frmA, ADH5, adhC (XM_007141519.1) (PHAVU_008G207900g), CML (XM_007137421.1) (PHAVU_009G130600g), SAUR (XM_007131669.1) (PHAVU_011G037100g), SAUR (XM_007131677.1) (PHAVU_011G037900g), PXG (XM_007153204.1) (PHAVU_003G0207000g) and GH3 (XM_007155808.1) (PHAVU_003G238700g). The remaining eight genes (except for XM_007143360.1) can be annotated in KEGG, and nine genes can be annotated in the GO function database. Six of these nine genes encode transcription factors. The specific results are shown in Table 3.

Table 3. The KEGG enrichment of DEGs for qRT–PCR validation. The expression difference multiple log2 (foldchange) ≥ 2, false discovery rate (FDR) <0.01, and p value <0.05 were used as screening criteria.

Gene	Gene name	Up/down	KEGG description	TF family
PHAVU_007G071000g	XM_007143360.1	up	–	G2-like
PHAVU_007G087800g	XM_007143566.1	down	EIX1_2; EIX receptor 1/2	C3H
PHAVU_007G187200g	XM_007144763.1	down	CALM; calmodulin	HB-other
PHAVU_008G207900g	XM_007141519.1	up	frmA, ADH5, adhC; S-(hydroxymethyl)glutathione dehydrogenase / alcohol dehydrogenase
PHAVU_009G130600g	XM_007137421.1	up	CML; calcium-binding protein CML	FAR1
PHAVU_011G037100g	XM_007131669.1	down	SAUR; SAUR family protein	bHLH
PHAVU_011G037900g	XM_007131677.1	down	SAUR; SAUR family protein	bHLH
PHAVU_003G0207000g	XM_007153204.1	up	PXG; peroxygenase [EC:1.11.2.3]
PHAVU_003G238700g	XM_007155808.1	up	GH3; auxin responsive GH3 gene family

After comparing the results of qRT–PCR and RNA-seq, we found that the expression levels of these nine genes were different, but the overall expression trends of the genes were largely the same, which is indicative of the biological meaningfulness of the transcriptome sequencing results across the two common bean varieties (Figure 7).

Figure 7. Relative expression levels of DEGs by RNA-seq (A) and qRT–PCR (B).

Discussion

Relationship between antioxidant enzyme activity and drought resistance in common soybean

As a global issue that is increasing in severity and frequency within the context of climate change, drought is worthy of in-depth study, especially in the context of agriculture. The drought resistance of common bean has a long evolutionary history in most areas where it occurs [33, 34]. Plant growth and development are slowed down under drought stress [35]; additionally, the intrinsic antioxidant defense system will be activated, increasing antioxidant enzyme activity, thereby reducing damage to the plant [36–38]. In this study, we showed that the SOD and POD activities of common bean decreased under drought stress and then increased after rewatering; these findings are inconsistent with those of Wang et al. [39]. This study showed that the SOD activity in the roots was lower than that in the leaves and that the SOD activity in the roots changed more slowly with increasing treatment time than that of leaves, but the POD activity of roots was much greater than that of leaves. In this study, CAT activity first decreased and then increased, this result is in agreement with a previous report [40]. In this study, the soluble protein, soluble sugar and proline contents in common bean leaves and roots under PEG stress increased, consistent with a previous study [41]. In this study, the soluble protein and proline contents in the leaves were significantly higher than that of roots, but the soluble sugar content in the leaves was significantly lower than in the roots. Under drought stress, the MDA content of common bean plants increased, similar to a previous observation [42], and the MDA content in leaves was significantly lower than that of roots.

Screening of common bean for drought tolerance by RNA-Seq

With the continuous improvement of analytical techniques within the context of molecular biology, transcriptome sequencing technology will continue being widely used in plant and animal genomics research. In this study, RNA-Seq was performed to analyse the transcriptome of the two common bean varieties under drought stress. Overall, this approach enabled the construction of a transcriptome database of the common bean varieties, ‘HZYTJ’ and ‘WZWJD’, under drought stress. We found 1,361 DEGs in ‘HZYTJ’, of which most were upregulated, and there were 1,447 DEGs in ‘WZWJD’, with most of them being downregulated. The differentially expressed genes were mainly involved in biological metabolism. In the GO function database, the metabolic processes with most DEGs were linked to abiotic stimulus responses and oxalate metabolic processes. This shows that under drought stress, different common bean varieties could have the same metabolic response. There are many different genes involved in the process of signal transduction, reaction to chemical substances and reaction to oxygenated compounds in ‘HZYTJ’, and there are also small molecule catabolism processes in ‘WZWJD’. There are numerous genes involved in organic acid metabolism and carboxylic acid metabolism by GO and KEGG enrichment analysis, which indicates that under drought stress, the more drought-resistant bean variety and the more drought-sensitive bean variety employ different metabolic processes. By analysing the KEGG metabolism pathway, we found 10 significant metabolic pathways in ‘HZYTJ’, involving a total of 215 genes; there were nine significant metabolic pathways in ‘WZWJD’, involving a total of 185 genes. The metabolic pathways shared between the two varieties are photosynthesis-antenna protein, and glyoxylic acid and dicarboxylic acid metabolism pathways. This shows that under drought stress, the metabolism of these two pathways across the two varieties would change. Through an analysis of the factors, we found that bHLH transcription factors were the most abundant, followed by transcription factor families such as ERF, MYB-related, B3 and C2H2. During qRT–PCR, nine genes were randomly selected for verification, and the results of this verification were consistent with those of the transcriptome sequencing results, indicating that the transcriptome sequencing was reliable and therefore yielded biologically meaningful results.

The EIX1-2 (XM_007143566.1) gene is a receptor protein and belongs to the C3H family of transcription factors. It participates in the defence against plant pathogen under drought stress in Arabidopsis [43] and tomato (Solanum lycopersicum) [44–46]. The CALM (XM_007144763.1) gene encodes calmodulin, and it belongs to the HB-other family of transcription factors. This gene participates in the metabolic process of the MAPK signalling pathway in plants. CALM is involved in the antioxidant effect of Cassia leaves [47]. The expression of CALM is upregulated in the immune, detoxification or apoptotic response caused by Cr-induced plant organ damage [48]. A study found that CALM is involved in the Ca²⁺ signal transduction pathway by studying the resistance response mechanism of Malus halliana seedlings under salt-alkali stress [49]. The frmA, ADH5 and adhC (XM_007141519.1) genes are S-(hydroxymethyl) glutathione dehydrogenase/alcohol dehydrogenase and are involved in carbon metabolism. Searching through the corresponding KO number on the KEGG pathway, we found that many studies have been reported [50, 51]. The calcium-binding protein CML (XM_007137421.1) is in a family of transcription factors, FAR1. This protein is involved in plant–pathogen interactions. The CML gene has been studied in many plants. Dopamine can activate the Ca²⁺ signal transduction pathway by increasing the expression of CNGG and CAM/CML family genes, thereby improving the drought tolerance of apples [52]. Dubrovina et al. [53] cloned 54 CML mRNA transcripts from wild grapes and found that the grapes have significant stress tolerance. Through real-time quantitative RT–PCR, they found that the CML gene family actively responded to abiotic stress. Another study screened the CML gene in the calcium signal pathway by studying the transcriptome expressed during the low temperature-induced dormancy of garlic stems [54]. A study on the defence response of Chinese cabbage to soft rot through transcriptome analysis found that CML is a downstream immune defence gene triggered by the pathogen mode [55].

SAUR (XM_007131669.1) and SAUR (XM_007131677.1) are both SAUR family proteins that belong to the bHLH class of transcription factors. These proteins participate in the metabolic process of plant hormone signal transduction. For example, the SAUR gene was upregulated in response to barley waterlogging [56]. SAUR was also involved in the response of yam tuber expansion [57]. PXG (XM_007153204.1) is a peroxidase that participates in the biosynthesis pathway of cutin, amber and wax. Analysis of the evolution and function of PXG in fungi and other nonplant clades shows that PXG in fungi and plants have similar but different effects [58]. The XM_007143361.1 gene was not annotated in KEGG. It is annotated as an oxalate-CoA ligase in the UniProt database. The GH3 (XM_007155808.1) gene belongs to the auxin-responsive GH3 gene family and is involved in signal transduction. To date, many GH3 families have been identified across different plant species, such as the apple (Malus domestica) [59], maize (Zea mays) [60], citrus (Citrus sinensis) [61], rice (Oryza sativa) [62] and Kiwifruit (Actinidia chinensis) [63].

With the continued improvement of transcriptome technology, the cost of sequencing has gradually decreased; therefore, whole genome sequencing of the common bean has been completed. These advancements have allowed for more in-depth research into the molecular mechanisms underlying the drought resistance of the common bean. This study lays a foundation for the further study of the drought resistance of the common bean.

Conclusions

This study analyzed the physiological responses of two common soybean varieties under PEG treatment. The transcriptome analysis of the two varieties showed that 1,361 and 1,447 genes were differentially expressed in ‘HZYTJ’ (strong drought resistance) and ‘WZWJD’ (weak drought resistance), respectively. The selected and screened genes were validated by qRT–PCR. This study will provide candidate genes for genetic engineering towards drought resistance breeding of common bean.

Authors’ contributions

S.Z. and D.Z. conceived and designed the study with the help of M.L.; S.Z., D.Z., T.L., T.W., G.X. performed the experiments and analyzed the data; S. Z and D.Z. drafted and wrote the manuscript; S. Z, D.Z. and M.L. revised the manuscript. All authors have read and agreed to the published version of the manuscript.

Acknowledgements

We thank Bioeditas Technology Corporation (Shaanxi, China; www.bioeditas.cn) for performing RNA sequencing and providing technical assistance during the data analysis.

Disclosure statement

No potential conflict of interest was reported by the authors. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Data availability statement

The raw RNA sequencing data to NCBI at (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA946172) under accession numbers (SRX19727897-SRX19727914). All data that support the findings reported in this study are available from the corresponding author Meilan Li upon reasonable request ([email protected] or [email protected]).

Additional information

Funding

This work was supported by Shanxi Province’s key research and development projects, dryland vegetable bean quality efficiency machine integrated green production technology research and demonstration (201703D211002-7-4).