Identification and transcriptomic profiling of salinity stress response genes in colored wheat mutant

Plant materials

One of the progenies resulting from the cross between ‘Woori-mil’ (obtained from the National Agrobiodiversity Center, RDA, Korea; accession no. IT172221) and ‘D-7’ (an inbred line developed by Korea University; Fleming4/3/PIO2580//T831032/Hamlet) exhibited color segregation. This specific progeny, carefully chosen, had spikes containing colored seeds, and these were utilized in the subsequent generation. Finally, we developed the common wheat (Triticum aestivum L., 2n = 6x = 42, AABBDD) inbred line K4191 (hereafter termed PL1), distinguished by its deep purple grain color. K4191 was derived from the F4:8 generation resulting from the cross between ‘Woori-mil’ and ‘D-7’, both of which have common seed color (Hong et al., 2019). To induce genetic variation and diversity the population of colored wheat, colored wheat seeds (PL1) were irradiated with 200 Gy gamma rays at a dose rate of 25 Gy/h using a ⁶⁰Co gamma irradiator (150 TBq of capacity; Noridon, Ottawa, ON, Canada) at the Korea Atomic Energy Research Institute. Subsequently, the irradiated seeds were planted at the radiation breeding research farm. Briefly, 1500 M₀ seeds were exposed to irradiation, and the resulting seeds were sown to generate the M₁ generation. Among these seeds, 287 phenotypically distinctive lines were carefully selected with one spike per plant, and mutation breeding spanning from M₀ to M₄ was performed as thoroughly described in a previous study (Hong et al., 2019). The resulting mutants were continuously cultivated up to the M₆ generation and carefully selected based on excellent agricultural traits, including flowering time, plant height, yield, and grain color. In total, 50 mutant lines displaying stable phenotypes for at least two generations were chosen for further salt-tolerance screening.

Salt stress treatment

The PL1 (control line, K4191) and PL6 (mutant line) seeds were surface-sterilized with 70% ethanol for 1 min and then washed with sterile distilled water. Subsequently, the seeds were placed on moist filter papers in a Petri dish (SPL Life Sciences) until the first leaf of the seedlings appeared. Next, the uniformly germinated seeds were transferred to Incu Tissue culture vessels (SPL Life Sciences) filled with half-strength Hoagland’s culture solution (Sigma-Aldrich, St. Louis, MO, USA). The solutions were replaced daily. The seedlings were grown for 7 days in a well-controlled chamber at 22 °C and 60% humidity, with a photoperiod regime of 16/8 h day/night at 200–300 µmol m⁻²s⁻¹ light. After 7 days of transplanting, the seedlings were subjected to a salt stress treatment of a total volume of 200 mL of the solution containing 150 mM NaCl. Following treatment with 150 mM NaCl, the wheat leaves were collected at 0, 3, 24, and 48 h. Both control and salt-stressed seedlings were collected individually. The samples were immediately frozen in liquid nitrogen and stored at −80 °C until use in further experiments.

Measurement of leaf Na⁺ and K⁺ contents

The wheat leaves were collected separately and immediately frozen in liquid nitrogen. Subsequently, the samples were freeze-dried for 3 days in a Freeze Dry System (IlshinBioBase, Dongducheonsi, Gyeongi, Korea). The freeze-dried samples were then finely ground into a powder using a mortar and pestle. For further analysis, 50 mg of the freeze-dried samples was weighed using an analytical balance and boiled for 2 h at 200 °C in three mL of HNO₃ (70%, v/v) for digestion. After digestion, the extracted samples were diluted with 5% HNO₃ and filtered through a hydrophilic polytetrafluoroethylene syringe filter (0.45-μM pore size, 25-mm diameter). The shoot Na⁺ and K⁺ contents were measured using inductively coupled plasma atomic emission spectroscopy (ICP-AES, 720 series; Agilent, Santa Clara, CA, USA) and quantitatively analyzed using a VistaChip CCD detector (Agilent).

Measurement of chlorophyll content

To determine the chlorophyll content, wheat seedling samples were extracted with 100% methanol at 4 °C. The sample extracts were then subjected to centrifugation at 12,000 × g for 10 min, and the supernatant was used for chlorophyll content analysis. The total chlorophyll, chlorophyll a, and chlorophyll b concentrations were determined by measuring the absorbance at 644.8 and 661.6 nm using a UV–VIS spectrophotometer (Lichtenthaler, 1987). The chlorophyll concentration was calculated using the following equations: (1)${\displaystyle C_{\mathrm{a}}=11.24\times A_{\mathrm{661.6}}-2.04\times A_{644.8}}$ (2)${\displaystyle C_{\mathrm{b}}=20.13\times A_{\mathrm{644.8}}-4.19\times A_{661.6}}$ (3)${\displaystyle C_{\mathrm{total}}=7.05\times A_{661.6}+18.09\times A_{644.8}}$ where C_a, C_b, and C_total denote the concentrations of chlorophyll a, chlorophyll b, and total chlorophyll, respectively.

RNA sequencing and gene expression analyses

Total RNA was extracted from the wheat leaves of both PL1 and PL6 at each timepoint (0, 3, 24, and 48 h) using TRIzol (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. Two independent biological replicates were performed for each timepoint and line to ensure the reliability and reproducibility of the RNA-seq data. Additionally, the extracted RNA samples were treated with DNase I to remove any potential genomic DNA contamination. The RNA quality was assessed using an Agilent 2100 bioanalyzer (Agilent Technologies, Amstelveen, The Netherlands), and RNA quantification was performed using an ND-2000 Spectrophotometer (Thermo Inc.; Wilmington, DE, USA). For constructing the RNA-seq paired-end libraries, 10 µg of total RNA extracted from the samples was used with the TruSeq RNA Sample Preparation Kit (Catalog #RS-122-2001; Illumina, San Diego, CA, USA). The mRNA was isolated using a Poly(A) RNA Selection Kit (LEXOGEN, Inc.; Vienna, Austria) and reverse-transcribed into cDNA following the manufacturer’s instructions. The libraries were assessed using the Agilent 2100 bioanalyzer, and the mean fragment size was evaluated using a DNA High Sensitivity Kit (Agilent, Santa Clara, CA, USA). High-throughput sequencing was conducted using the HiSeq 2000 platform (Illumina). Before alignment, adaptor sequences were removed, and sequence quality was evaluated using the Bbduk tool (minimum length >20 and Q >20; https://jgi.doe.gov/data-and-tools/software-tools/bbtools/bb-tools-user-guide/bbduk-guide/). The reads were aligned to the wheat genome sequence provided by the International Wheat Genome Sequencing Consortium (IWGSC) wheat reference sequence (IWGSC Reference Sequence v1.0; https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.0/) using the HISAT2 alignment program with default parameters (Kim, Langmead & Salzberg, 2015). Reads mapped to the exons of each gene were enumerated using the HTSeq v0.6.1 high-throughput sequencing framework (Anders, Pyl & Huber, 2015). Subsequently, the differentially expressed genes (DEGs) under salt stress and control conditions were identified using the EdgeR package (Robinson, McCarthy & Smyth, 2010). Upregulated and downregulated genes with a p-value of <0.05, false discovery rate (FDR) of <0.05, and an absolute fold change value of >2 were used for downstream functional analysis. The log2-transformed transcript per million values were calculated using TPMCalculator (Vera Alvarez et al., 2019), and heatmaps of DEGs under control and stress conditions were generated. Local BlastX was used with peptide sequences of the Poaceae family retrieved from the National Center for Biotechnology Information (NCBI) database using an e-value threshold of 1 × 10⁻⁵ to annotate the DEGs. For gene expression analysis, total RNA was used to synthesize first-strand cDNA using the Power cDNA Synthesis Kit (iNtRON Biotechnology, Gyeonggi-do, Korea). Reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was performed in a total volume of 20 μL containing 1 μL of cDNA template, 0.2 μM primers, and 10 μL of TB Green Premix Ex Taq II (Takara, Kusatsu, Shiga, Japan). RT-qPCR was conducted using a CFX96TM Real-time PCR system (Bio-Rad, Hercules, CA, USA) with the following program: 95 ° C for 5 min, followed by 40 cycles at 95 °C for 10 s and 65 °C for 30 s. Actin (AB181991) was used as an internal control. The primers used in this experiment are listed in Table S1.

Functional analysis of DEGs

All expressed genes under both control and stress conditions were subjected to Gene Set Enrichment Analysis (GSEA) using the GSEA software (Subramanian et al., 2005). The gene matrix transposed file format (.GMT) of wheat was downloaded from g:Profiler (https://biit.cs.ut.ee/gprofiler/gost), a web server for functional enrichment analysis and gene list conversion (Raudvere et al., 2019). The enrichment score of each gene set was calculated using the full ranking, and the normalized enrichment score (NES) was determined for each gene set. The GSEA results, including rank, expression, and class files, were visualized as a network using Enrichment Map (Merico et al., 2010). For Kyoto Encyclopedia of Genes and Genome (KEGG) pathway enrichment analysis, the KEGG Orthology Database in KOBAS-i was used to predict the putative pathways of DEGs (Bu et al., 2021). The plant transcription factor data were obtained from the Plant Transcription Factor Database (PlantTFDB) (Tian et al., 2020). Protein Basic Local Alignment Search Tool (BLASTP) was used on the peptide sequences of the DEGs, based on the local transcription factor database obtained from PlantTFDB, with an E-value threshold of 1 × 10⁻¹ and sequence identity of >80%. Mev software (http://sourceforge.net/projects/mev-tm4/files/mev-tm4/) was used for k-means clustering of DEGs identified from the GSEA, KEGG pathway, and transcription factor analyses. The results of the GSEA and KEGG pathway analysis were generated using an R script and the ggplot2 R package. Additionally, MapMan was used to identify the pathways of stage-specific genes (Sreenivasulu et al., 2008).

Enzyme activities assays

The crude enzyme was extracted from 100 mg of wheat leaves using a protein extraction buffer containing 50 mM potassium phosphate buffer (pH 7.5). The activities of catalase (CAT), peroxidase (POD), and superoxide dismutase (SOD) and total antioxidant activity (TAC) were measured using commercially available assay kits. Specifically, CAT activity was determined using a catalase microplate assay kit (kit number: MBS8243260; MyBiosource, Inc., San Diego, CA, USA), POD activity was measured using a POD assay kit (kit number: KTB1150; Abbkine, Inc., Wuhan, China), and SOD activity was estimated using a total SOD activity assay kit (WST-1 method) (kit number: MBS2540402; MyBiosource, Inc., San Diego, CA, USA). TAC was assessed using a TAC assay kit (kit number: MAK187; Sigma-Aldrich, St. Louis, MO, USA). The preparation of the reaction mixture and the calculations for each measurement were performed as described in the respective protocol books provided with each assay kit.

Characteristics of the salt-tolerant colored wheat mutant induced via gamma irradiation

Throughout the mutation breeding process, detailed records of agricultural traits, including the flowering time, plant height, and yield, were meticulously collected for the mutant lines. These data allowed for a comprehensive assessment of the phenotypic characteristics of the lines. Evidence supporting the stable phenotype of the mutant lines is provided in Fig. S1, which also presents the field performance of PL1 and PL6. Additionally, the difference in grain color between the colored wheat lines used in this study is illustrated in Fig. S2. To select salt-resistant wheat lines, 50 wheat mutant lines were treated with 150 mM NaCl, and germination rate, shoot length, and root length were subsequently investigated. The results showed that 27 lines exhibited a germination rate of 90%. Among these, six lines demonstrated a growth increase of approximately 20% or more compared to the control group subjected to salt treatment (Fig. S3). Through preliminary salt-tolerance screening, PL6 was selected as the gamma ray-derived mutant line that exhibited favorable salt-tolerance characteristics (Fig. S3). To assess the growth response of the control line (PL1) and PL6 under varying salt concentrations, the seeds were treated with NaCl solutions of 50, 100, 150, 200, 250, 300, and 500 mM, along with distilled water as the control (Choudhary et al., 2021). Overall, high salt concentrations negatively affected seed germination and seedling growth (Figs. 1A and 1B). The germination percentage and seedling growth were reduced with increasing salt concentration in both PL1 and PL6 (Table S2). However, PL6 demonstrated higher germination percentages, particularly at the maximum NaCl concentration, exceeding those of PL1. Remarkably, a maximum increase of 20% in germination was observed for PL6 following treatment with 250 mM NaCl. Moreover, PL6 consistently outperformed PL1 in terms of seedling growth under all salt treatment conditions, as evidenced by its longer shoot and root lengths (Figs. 1C and 1D). The comprehensive data strongly indicates that the gamma ray-derived mutant PL6 exhibits higher resistance to salt stress than PL1.

Assessment of Na⁺, K⁺, and chlorophyll contents under salt stress conditions

Prior to treatment, PL6 had a higher Na⁺ ion content than PL1 (Fig. 2A). However, with increasing time of exposure to salt stress, the Na⁺ ion content markedly increased in both PL1 and PL6. Notably, the rate of increase in Na⁺ ion content was lower in PL6 than in PL1. Conversely, the K⁺ ion content steadily decreased with salt treatment in both PL1 and PL6 (Fig. 2B). To further analyze the ion contents, we calculated the relative ratios of K⁺ and Na⁺ ions in PL1 and PL6, considering their respective contents under control conditions (Figs. 2C and 2D). In PL1, the Na⁺ ion content increased significantly by 47 times from the baseline (0 h) to 48 h following salt treatment. In contrast, PL6 exhibited a milder increase in Na⁺ ion content, approximately 20 times higher at 48 h after salt stress. Consequently, the relative Na⁺ content was more profoundly affected by salt stress in PL1 than in PL6. Interestingly, the chlorophyll concentrations of both PL1 and PL6 remained relatively stable under salt stress (Fig. 2E), indicating that they were not significantly affected by the imposed salinity conditions.

DEGs during salt stress

After treatment with 150 mM NaCl, leaves were harvested from PL1 and PL6 at 0, 3, 24, and 48 h and subjected to RNA sequencing. Following quality evaluation and trimming, an average of 38.1 million trimmed reads and over 22.1 billion bases were generated from each sample under both control and salt stress conditions. The average percentage of Q20 and Q30 bases was found to be 98.4% and 95.5%, respectively, indicating high sequencing quality. Moreover, >96% of the sequenced data exhibited an average mapping rate of 96.16%, successfully aligning to the IWGSC wheat reference sequence (Table S3). During data analysis, a total of 4,017 DEGs were identified with a p-value of <0.05, FDR of <0.05, and absolute fold change value of >2 (Fig. 3A and Table S4). Specifically, in PL1, 872, 1,588, and 1,080 DEGs were identified at 3, 24, and 48 h after salt treatment, respectively, in comparison to the untreated condition (0h) (Fig. 3B). Similarly, for PL6, the number of DEGs was 566, 1,248, and 1,810 at 3, 24, and 48 h after salt treatment, respectively (Fig. 3C). These results highlight the dynamic gene expression changes in PL1 and PL6 under salt stress at different timepoints, contributing to a better understanding of the underlying molecular responses to salt stress in these wheat lines.

Functional analysis of the DEGs during salt stress

Overall, 33 GO terms were identified for each treatment condition (Fig. 3D and Table S5). Notably, several gene sets, including defense response (GO: 0006952), glutathione metabolic process (GO: 0006749), peroxidase activity (GO: 0004601), ROS metabolic process (GO: 0072593), response to biotic stimulus (GO: 0009607), and response to stress (GO: 0006950), were positively correlated with salt stress and PL6, exhibiting a positive NES (Fig. 3D). To visualize the results, all the gene sets from the GSEA were organized into four networks using Enrichment Map (Merico et al., 2010) (Figs. 4A–4E). The expression patterns of each network in Figs 4A–4E for PL1 and PL6 were clustered by expressed patterns (Figs. 4F–4J, respectively). The K-means clustering algorithm in the Mev software was used to identify the clusters of DEGs in each GO term under control and salt stress conditions based on their expression patterns. Most of the expression patterns from the identified clusters did not differ between the control and salt stress conditions. Three clusters that demonstrated different expression patterns for PL1 and PL6, especially those upregulated in PL6, were selected and marked in red boxes in Figs. 4F, 4G, and 4I, and a heatmap of the genes from these clusters was generated (Fig. 4K). Plant hormone-related genes (TRAESCS1B02G145800 and TRAESCS1B02G138100), ROS-related genes (TRAESCS1B02G059100, TRAESCS1B02G095800, TRAESCS1B02G096200, TRAESCS1B02G096900, and TRAESCS1B02G115900), and stress-response genes (TRAESCS5D02G492900, TRAESCS1A02G009900, TRAESCS1B02G023000, and TRAESCS2A02G037400) were highly expressed in PL6 under salt stress conditions. Furthermore, six genes related to chromatin remodeling (TRAESCS1B02G048900, TRAESCS1B02G049100, TRAESCS1D02G286700, TRAESCS1B02G149000, and TRAESCS7D02G246600) showed high expression patterns in PL6 under salt stress conditions [Table 1]. A high number of transcriptomes of MADS-box transcription factors (TRAESCS4A02G002600, and TRAESCS6D02G293200) were also detected in PL6 under salt stress. An auxin-responsive protein (TRAESCS1B02G138100) and probable histone H2A variant 3 (TRAESCS7D02G246600) were also found in cluster 4 [Table 1].

In the case of the differences in the KEGG pathways between PL1 and PL6 under salt stress conditions, the rich factor of “Biosynthesis of secondary metabolites” in PL6 after 3 h of salt stress was ∼0.05, increasing to ∼0.2 after 48 h of salt stress (Fig. 5). Likewise, the rich factors of “Flavonoid biosynthesis” were 0.17 and 0.23, after 24 and 48 h of salt stress, respectively. This was only observed in PL6 during salt stress conditions (Fig. 5 and Table S6).

In addition to GO and KEGG analysis, the role of DEGs as transcription factors was investigated. DEGs at different timepoints under salt stress in PL1 and PL6 were identified using PlantTFDB (http://planttfdb//planttfdb.gao-lab.org). In total, 255 genes were identified with an e-value threshold of 1 × 10⁻¹ and a sequence identity of >80% and further selected to compare the expression patterns between PL1 and PL6 under salt stress conditions. The most abundant type of transcription factor was the ethylene-response factor (ERF) protein family, followed by the basic helix-loop-helix (bHLH) protein family; heat shock transcription factor protein family; myeloblastosis (MYB)-related protein family; and Nam, ATAF, and CUC (NAC) protein family (Fig. 6A). Moreover, 255 putative transcription factors were grouped by expression pattern into six clusters and an unclassified group (Fig. 6B). Overall, 72, 44, and 35 DEGs were annotated by the ERF, bHLH, and MYB (related) protein families, respectively. These three transcription factors accounted for 59% of the total number of transcription factors. The expression patterns of DEGs in clusters 2 and 6 (marked with red boxes in Fig. 6B) were selected and expressed in heatmaps (Fig. 6C) to display differences in the expression patterns of DEGs between PL1 and PL6 under salt stress conditions (Table S7). Notably, PL6 exhibited higher expression of specific transcription factors under salt stress conditions than PL1, as displayed in the heatmap (Fig. 6C).

Lastly, 22 protein kinase genes were identified with significant expression patterns at different timepoints, including two calcineurin B-like (CBL)-interacting protein kinases and one mitogen-activated protein kinase (MAPK) with more than two-fold changes in PL6 under salt stress (Table 2). Additionally, 70 differentially expressed salt stress-responsive genes involved in regulating the circadian clock system, cytoskeleton organization, and cell wall organization were identified using MapMan, with 15 of them showing more than a two-fold change in PL6 (Table 3).

Enzyme activities assays

To investigate the differences in enzyme activities between PL1 and PL6 under salt stress conditions, we measured CAT, POD, and SOD activities and TAC (Fig. 7). Upon subjecting both wheat lines to salt stress, we observed distinct patterns in enzyme activities. In PL6, CAT and POD activities significantly increased after 24 and 48 h of exposure to salt stress (Fig. 7A and 7B). Conversely, in PL1, SOD activity slightly decreased after 24 and 48 h exposure to salt stress (Fig. 7C). Furthermore, the TAC in PL1 was not significantly changed by salt stress (Fig. 7D). Conversely, in PL6, the TAC notably increased after 24 and 48 h of exposure to salt stress. This increase in TAC suggests that PL6 has a higher capacity to counteract oxidative stress and maintain cellular redox balance than PL1, contributing to its enhanced salinity tolerance.

Validation of the DEG results using reverse transcription-quantitative polymerase chain reaction

Supporting the DEG results, 12 genes from the three aforementioned clusters from PL1 and PL6 were selected for RT-qPCR (Fig. 8). All the selected genes were more highly expressed in PL6 than in PL1. peroxidase 2 (TRAESCS1B02G095800), nitrate transporter (TRAESCS1B02G038700), auxin-responsive protein (TRAESCS1B02G138100), and replication protein A (TRAESCS1B02G102200) transcripts in PL6 were highly expressed at 48 h following salt treatment (Fig. 8A). Nuclear transport factor 2- like protein (TRAESCS2A02G046200), histone H2A (TRAESCS1B02G048900), integral membrane protein (TRAESCS1B02G071800), and histone H2A variant 3 (TRAESCS7D02G246600) transcripts in PL6 continuously decreased at 24 h and peaked at 48 h following salt treatment (Fig. 8B). Argonaute 1C-like isoform X2 (TRAESCS6B02G466700), MADS-box (TRAESCS6B02G017900), and aspartokinase 1 (TRAESCS5D02G537600) transcripts in PL6 peaked at 3 h, and all gradually decreased, except for ribosome biogenesis protein NOP53 (TRAESCS1B02G105100) (Fig. 8C). These results are consistent with those of RNA sequencing (RNA-seq).