Introduction
In the Middle East, especially in Iran, wheat (Triticum aestivum L.) is the most important food crop. However, in a span of 30 years (i.e. 1960–1990), salt-affected lands have increased by more than 48% (15.5–23 Mha) in Iran, and the latest estimates indicate that this figure has reached more than 25 million hectares (FAO, 2019). It was estimated that over 12% of irrigated agricultural land out of about 34 million ha of the cultivable area in Iran was affected by salt, resulting in significant yield losses and total annual damage of more than 1 billion USD (Qadir et al., Reference Qadir, Qureshi and Cheraghi2008). Even though wheat is moderately tolerant to salinity, yield losses have been reported up to 40% (Qadir et al., Reference Qadir, Quillérou, Nangia, Murtaza, Singh, Thomas, Drechsel and Noble2014). Salinity refers to the amount of accumulated and dissolved salts in soil and water (Grattan, Reference Grattan2002). High salt concentrations affect the growth and development of plants (Yadav et al., Reference Yadav, Bharadwaj, Nayak, Mahto, Singh and Prasad2019). Therefore, plants need to adapt to stress-inducing conditions. It is becoming increasingly difficult to ignore salinity as a major threat to global food security, causing significant reductions in plant productivity and wheat production worldwide (Oyiga et al., Reference Oyiga, Sharma, Baum, Ogbonnaya, Léon and Ballvora2018). The response of plants to salt stress is highly complex, and it depends on plant species, growth stage, salt concentration, types of ions, and various environmental factors (Ashraf, Reference Ashraf2004; Sairam and Tyagi, Reference Sairam and Tyagi2004; Munns and Tester, Reference Munns and Tester2008; Rao et al., Reference Rao, Din, Akhtar, Sarwar, Ahmed, Rashid, Azmat Ullah Khan, Qaisar, Shahid, Ahmad Nasir, Husnain and Abdurakhmonov2016; Yadav et al., Reference Yadav, Bharadwaj, Nayak, Mahto, Singh and Prasad2019). These responses may include tolerance to osmotic stress, Na+ exclusion, and tissue tolerance, which is well demonstrated in Munns and Tester (Reference Munns and Tester2008). In addition, many genes are involved in each mechanism, some of which are identified and some still remain unknown. In addition, by the invention of next-generation sequencing (NGS), sequencing can be performed faster, cheaper, more accurately, and with much fewer limitations than previous sequencing technologies, such as Sanger. As a result, many new genes have been identified in different species. Over the past decade, there have been increasingly rapid advances, and more researchers are interested in sequencing technology. NGS has been used as a highly accurate and sensitive tool for determining gene expression patterns across different studies (Cho et al., Reference Cho, Garvin and Muehlbauer2006; Ge et al., Reference Ge, Li, Zhu, Bai, Lv, Guo, Ji and Cai2010; Guimarães et al., Reference Guimarães, Brasileiro, Morgante, Martins, Pappas, Silva, Togawa, Leal-Bertioli, Araujo, Moretzsohn and Bertioli2012; Zhang et al., Reference Zhang, Liao, Chang and Zhang2014; Bahieldin et al., Reference Bahieldin, Atef, Sabir, Gadalla, Edris, Alzohairy, Radhwan, Baeshen, Ramadan, Eissa, Hassan, Baeshen, Abuzinadah, Al-Kordy, El-Domyati and Jansen2015; Li et al., Reference Li, Li, Chen, Tang, Li and Huang2016; Zhou et al., Reference Zhou, Yang, Cui, Zhang, Luo and Xie2016; Moazzzam Jazi et al., Reference Moazzzam Jazi, Seyedi, Ebrahimie, Ebrahimi, De Moro and Botanga2017; Wang et al., Reference Wang, Geng, Zhu, Li, Tong, Wang, Chen, Tang and He2017; Wu et al., Reference Wu, Hu, Huo, Zhang, Chen and Zhang2017; Mansouri et al., Reference Mansouri, Naghavi, Alizadeh, Mohammadi-Nejad, Mousavi, Hosseini Salekdeh and Tada2018; Wu et al., Reference Wu, Zhao, Wu, Yuan, Ma, Lin, Pan, Li and Sun2019), especially in response to abiotic stress conditions, such as salinity, drought, cold and heat (Qin et al., Reference Qin, Wu, Peng, Yao, Ni, Li, Zhou and Sun2008; Winfield et al., Reference Winfield, Lu, Wilson, Coghill and Edwards2010; Takahashi et al., Reference Takahashi, Tilbrook, Trittermann, Berger, Roy, Seki, Shinozaki and Tester2015; Goyal et al., Reference Goyal, Amit, Singh, Mahato, Chand and Kanika2016; Wani et al., Reference Wani, Tripathi, Zaid, Challa, Kumar, Kumar, Upadhyay, Joshi and Bhatt2018; Amirbakhtiar et al., Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019; Chaichi et al., Reference Chaichi, Sanjarian, Razavi and Gonzalez-Hernandez2019). Moreover, this method provides an opportunity to study the relationship between wheat varieties in different regions and to evaluate the mechanisms associated with the adaptability of the plants, leading to the possibility of crop production under various environmental conditions (Shavrukov et al., Reference Shavrukov, Suchecki, Eliby, Abugalieva, Kenebayev and Langridge2014). Since its appearance, this robust technology has considerably progressed in terms of read length and quality, speed, throughput and reduction in time and cost (da Fonseca et al., Reference Da Fonseca, Albrechtsen, Themudo, Ramos-Madrigal, Sibbesen, Maretty, Zepeda-Mendoza, Campos, Heller and Pereira2016). In a recent experiment focused on salt tolerance by Amirbakhtiar et al. (Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019), RNA sequencing was utilized, resulting in over 113 million reads, around 104,013 genes, 5128 differentially expressed genes (DEGs) and 227 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. However, our study focused on leaf transcriptome analysis of a famous local salt-tolerant wheat landrace (i.e. Roshan) under severe prolonged stress (i.e. 250 mM of NaCl for 48 h). The study resulted in 137.5 million clean reads, 214 KEGG pathways, 17,897 up-regulated and down-regulated DEGs, and a large number of genes, enzymes and transcription factors (TFs) involved in the response to salinity, which will be explained below. We aimed to evaluate the leaf transcriptome of a local salt-tolerant wheat landrace (i.e. Roshan, a local tall short-awn salt and drought-tolerant landrace) under 250 mM salt stress by applying the RNA deep sequencing paired-end technology and poly-A capturing method for analysing pathways to reveal potential regulatory genes and molecular mechanisms.
Materials and methods
Plant materials and salt stress conditions
In the present study, Roshan seeds (i.e. a local salt-tolerant wheat landrace) (Poustini and Siosemardeh, Reference Poustini and Siosemardeh2004; Dehdari et al., Reference Dehdari, Rezai and Maibody2005) were received from the Seed and Plant Improvement Institute (SPII), Karaj, Iran. First, they were germinated in Petri dishes with filter paper for 3 days, and then, the uniform seedlings were transformed into 2-cm holes made in Styrofoam sheets, which were on the top of pots with a diameter of 20 cm and a depth of 25 cm, using Hoagland solution (Hoagland and Arnon, Reference Hogland and Arnon1950). The experiment was a completely randomized design with three replicates. The growth conditions of the seedlings at phytotron (BINDER: KBWF 720(E 5.2)) were as follows: hydroponic Hoagland solution (Hoagland nutrient solution was increased gradually to ½ and full strength), 16 h of light and 8 h of dark, at 25 and 20°C temperatures, respectively, and with an air humidity of 60%. The pH was kept at 5.8, and the nutrient solution was renewed daily. After 14 days, 250 mM salt (NaCl) stress was imposed for 48 h. Finally, the total plant leaves were harvested 48 h after salt stress treatment. The leaves collected under both control and salinity stressed conditions were immediately transferred to liquid nitrogen, and then, they were kept in the freezer at −80°C for further analysis.
Total RNA isolation and mRNA extraction
Total RNA was extracted from the three biological replicates of normal and salt-stressed samples (after 48 h of exposure to salt stress) using RNeasy Plant Mini Kit QIAGEN based on the relevant instructions. The total RNA quality was determined by 1.2% formaldehyde agarose gel electrophoresis and the NanoDrop 2000c spectrophotometer. The concentration of total RNA, 28S/18S and 23S/16S test, and RNA integrity number (RIN) of samples were analysed using the RNA 6000 nano Reagents kit by Agilent 2100 Bioanalyser. Finally, samples with an OD260 nm/OD280 nm ratio of ~2 and a RIN >7 were chosen for RNA sequencing. The RNA samples were then sent to Beijing Genomics Institute (BGI), China (https://www.bgi.com) for library construction (according to the poly-A capturing method, i.e. the most well-known method used for RNA sequencing, based on the identification of poly-A tails at the 3′ end of the RNA (Guo et al., Reference Guo, Zhao, Sheng, Guo, Lehmann, Pietenpol, Samuels and Shyr2015)) and sequencing. The samples were sequenced using Illumina HiSeq 2000 with a read length of 150 bp.
RNA-seq data analysis
Raw data were pre-processed (the percentages of Q20 for all reads were more than 97%. Table 1 shows the statistical results after data treatment) using Trimmomatic (v0.35) (Bolger et al., Reference Bolger, Lohse and Usadel2014) to remove very short reads and low-quality sequences, with parameters including TRAILING:20, MAXINFO:120:0.9 and MINLEN:120. Data quality control was conducted using the FastQC (v0.11.5) toolkit (Andrews, Reference Andrews2017). After trimming, to map the NGS reads, clean reads were aligned against the IWGSC wheat reference genome v1 (https://plants.ensembl.org/Triticum_aestivum/Info/Index) using STAR v.2.5.3a (Dobin et al., Reference Dobin, Davis, Schlesinger, Drenkow, Zaleski, Jha, Batut, Chaisson and Gingeras2012) and Hisat2 v.2.1.0 (Kim et al., Reference Kim, Langmead and Salzberg2015) software applications. The outputs generated by STAR and Hisat2 were separately used as input for the HTSeq-count v 2.7.3 software (Anders et al., Reference Anders, Pyl and Huber2015) to quantify the gene expression levels. Differential expression analysis was performed using DESeq2 v 1.16.1 (Love et al., Reference Love, Huber and Anders2014) as follows: the genes obtained from Hisat2 and STAR were analysed separately using the DESeq2 package. The adjusted P-value < 0.05 (false discovery rate (FDR)) and fold-changes > 2(Log2(treated/control) ⩾ 1) or −2 > fold-change > 2(−1 > Log2(treated/control) > 1) were used as thresholds for DEGs. Finally, DEGs identified by the two methods (i.e. STAR + DESeq2 and Hisat2 + DESeq2) were employed for further analysis. The NGS raw data of the present study were submitted to SRA (sequence read achieve) of NCBI under access number PRJNA589280 (SRR10914949, SRR10914948, SRR10914947, SRR10914946, SRR10914945 and SRR10914944).
Statistical results (number of genes, raw and clean reads, bases and reads quality) after data treatment in exposure to salinity stress in local salt-tolerant wheat landrace (Roshan)
Functional enrichment analysis
The GO-seq R package (Young et al., Reference Young, Wakefield, Smyth and Oshlack2010) was used to classify DEGs according to gene ontology (GO) terms based on molecular function, biological processes and cellular components. A hypergeometric test was employed, and GO terms with FDR < 0.05 were selected as significant. The Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) was applied to perform pathway enrichment analysis of the DEGs. The same statistical test was used for KEGG pathway enrichment. We performed the TF family analyses using ‘PlantTFDB 4.0’ (http://planttfdb.cbi.pku.edu.cn) (Jin et al., Reference Jin, Tian, Yang, Meng, Kong, Luo and Gao2017).
Furthermore, MapMan (Version 3.6.0RC1; https://mapman.gabipd.org/web/guest/mapman-version-3.6.0) (Thimm et al., Reference Thimm, Blasing, Gibon, Nagel, Meyer, Kruger, Selbig, Muller, Rhee and Stitt2004) was employed with a P-value cut-off of ⩽ 0.05 to investigate changes caused by salinity stress in the pathways of the DEGs.
Validation of RNA deep sequencing data using quantitative real-time polymerase chain reaction (qRT-PCR)
To validate the RNA-Seq results, six DEGs were randomly selected for qRT-PCR analysis with a Thermal Cycler® Biorad real-time PCR Instrument (Corbbet-Rotor gene 6000) using REALQ PLUS 2X MASTER MIX SYBR Green (AMPLIQON) according to the manufacturer's protocol. The qRT-PCR amplification steps included pre-denaturation at 95°C for 30 s, followed by 35 cycles of 95°C for 5 s and 72–95°C for 30 s (TM depends on the genes). The gene-specific primers for RT-qPCR were designed using the Primer3 Plus online software (Untergasser et al., Reference Untergasser, Nijveen, Rao, Bisseling, Geurts and Leunissen2007) (www.bioinformatics.nl/primers3plus), and their sequences are presented in online Supplementary Table S1. First-strand cDNA synthesis was performed (1 μg of mRNA was used) using the SINACLON kit oligo(dT)18 (SINACLON, CAT. NO: RTS202) according to the manufacturer's instructions. PCR conditions included 1 μl cDNA, 5 μl REALQ PLUS 2X MASTER MIX SYBR GREEN (AMPLIQON), 2 mM of MgCl2 and 10 μl of primer forward and primer reverse each, in a final volume of 10 μl. Similar to previous studies (Amirbakhtiar et al., Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019; Chaichi et al., Reference Chaichi, Sanjarian, Razavi and Gonzalez-Hernandez2019), the actin-2 housekeeping gene was used as an internal control to normalize the gene expression levels. Samples were run in duplicate. The relative quantitative procedure (2−ΔΔCt) (Livak and Schmittgen, Reference Livak and Schmittgen2001) was used to calculate quantitative gene expression levels. The same RNA samples as those described for sequencing were used for RT-qPCR.
Results
Data production statistics
We performed transcriptome profiling in order to investigate molecular responses and identify the genes involved in the salt stress response in the leaves of a local salt-tolerant wheat landrace (i.e. Roshan). During the experiment, 14-day-old seedlings were exposed to 250 mM of NaCl for 48 h. Then, leaf samples were used for RNA extraction, and finally, six RNA samples were sent to BGI for high-throughput RNA-seq. The samples were sequenced using Illumina HiSeq 2000 with a read length of 150 bp. In total, 137,509,859 raw reads and 137,508,542 clean reads (~99% of the raw reads) were obtained.
All the clean reads were mapped against the Chinese Spring reference genome IWGSC_V1. Overall, the alignment rate showed that on average, about 89.60 and 93.17% of clean reads were mapped to the wheat reference genome by STAR and HISAT2, respectively (online Supplementary Table S2). Of these, ~82 and 77% were uniquely mapped via STAR and HISAT2, respectively. Moreover, 5.46% (HISAT2) and 7.54% (STAR) of the reads showed multi-position matches.
Analysis of DEGs
Differential expression analysis of the transcripts was performed using DESeq2 (online Supplementary Table S3). Only significant DEGs identified by both approaches (STAR + DESeq2 and Hisat2 + DESeq2) were subjected to further analysis. The number of common DEGs in both approaches for up-regulated and down-regulated genes was 7871 and 10,026, respectively (online Supplementary Table S3, Figs 1 and 2).
DEGs (up and down regulated) between control and salt stress in bread wheat by using HISAT2 + DESeq2 and STAR + DESeq2.
Fold change distribution of RNA-Seq data Processing and DEGs analysis.
Identification of TFs
In the present study, a total of 1857 DEGs, including 797 and 1060 up-regulated and down-regulated genes, respectively, representing 1.41% of the transcriptome, were classified into 87 TFs families (online Supplementary Table S4). Among them, GRAS (269 genes or 14.48%), C2H2 (186 genes or 10.01%) and CAMTA (176 genes or 9.47%) were the most abundant. In the present study, other large TFs families, such as bHLH, AP2/ERF, bZIP, NAC, MYB and WRKY were also affected by the salt stress conditions (i.e. 250 mM of NaCl). In addition, 110 bHLH protein-encoding genes were identified, of which 37 genes were significantly up-regulated and 73 were significantly down-regulated.
Moreover, 117, 66 and 99 DEGs, belonging to AP2/ERF, bZIP and WRKY families, respectively, were identified in the present study. Moreover, among the significant DEGs, 94 NAC genes were discovered, of which 61 NAC stress-responsive genes were up-regulated and 33 were down-regulated after exposure to salt stress. Finally, 97 wheat MYB proteins were found to be up-regulated under salt stress conditions.
KEGG pathway enrichment analysis
To better understand the functions of DEGs under salt stress, pathway network analysis was conducted using KEGG (KEGG, http://www.genome.jp/kegg/) (Kanehisa and Goto, Reference Kanehisa and Goto2000) (Fig. 3). According to the KEGG results, 7679 genes were divided into 214 KEGG pathways (online Supplementary Table S5), as 104 and 110 pathways belonged to up-regulated and down-regulated genes, respectively. Out of 214 pathways, only 62 pathways (including 5634 genes) were significant at corrected P-value ⩽ 0.05 (online Supplementary Tables S5 and S6, and Fig. S1). The most significant pathways included metabolic pathways, biosynthesis of secondary metabolites, peroxisome, fatty acid metabolism, glycerophospholipid metabolism, carbon metabolism, glutathione metabolism and galactose metabolism, respectively.
KEGG classification of identified DEGs was divided into five main categories (cellular processes; environmental information processing; genetic information processing; metabolism; organismal systems) in Triticum aestivum Landrace (Roshan) in response to salt stress conditions.
Overview of metabolism under salinity stress
In order to complete the analysis of genomic data and to better understand the plant's response to salinity stress, metabolic pathways were analysed using MapMan 3.6.0. According to the results presented in Fig. 4, 14,625 genes out of 17,897 DEGs were mapped, of which 3207 genes were involved in these pathways. However, some genes may also be mapped under multiple bincodes. Accordingly, photosynthesis (photophosphorylation, Calvin cycle, photorespiration), amino acid metabolism, cell wall organization (hemicellulose), lipid metabolism and protein homeostasis pathways have significantly been enriched in the studied genotype. The overview of the photosynthesis pathways revealed that genes encoding transketolase, glycolate oxidase and ELIP LHC-related protein were the most enriched. The overview of amino acid biosynthesis pathways indicated that genes encoding aspartate and pyruvate (valine/leucine/isoleucine aminotransferase) families were specifically enriched under salt stress. Based on the cell wall organization overview, the genes encoding mannan biosynthesis, modification and degradation were notably enriched. The overview of lipid metabolism pathways showed that genes encoding endoplasmic reticulum-plasma membrane tethering protein, steroleosin, caleosin, obtusifoliol 14-alpha demethylase and NAD-dependent glycerol-3-phosphate dehydrogenase were increasingly enriched after salt exposure. Moreover, genes associated with the proteolysis process, including FtsH plastidial protease complexes and asparaginyl endopeptidase (Legumain), were significantly enriched in the protein homeostasis overview (online Supplementary Table S7). In addition, we found that the applied salinity stress enriched the genes involved in hormone signalling pathways, including ABA, IAA, Jasmonate, ethylene, BA, GA and cytokinin (online Supplementary Table S8). Furthermore, the regulation overview (Fig. 5) demonstrated that genes contributing to the transcription regulation, including the sucrose nonfermenting-1 (SNF1)-related protein kinase (SnRK1) complex, were enriched under salinity stress (online Supplementary Table S8).
Metabolism overview of DEGs in local salt-tolerant wheat landrace (Roshan) under salt stress using Mapman 3.6.0. (green: up and red: down-regulated genes).
Regulation overview of DEGs in local salt-tolerant wheat landrace (Roshan) under salt stress using Mapman 3.6.0. (blue: up and red: down-regulated genes).
GO enrichment analysis in DEGs
To interpret the functions of DEGs, GO enrichment analysis was conducted using the GO-seq tool (Young et al., Reference Young, Wakefield, Smyth and Oshlack2010). Annotation of DEGs revealed a total of 226 GO terms (104 up-regulated and 122 down-regulated GO terms), of which 108 GO terms (including 3325 genes) were significant at a P-value ⩽ 0.05 (online Supplementary Table S9). Overall, in the ‘biological process’, ‘molecular function’ and ‘cellular component’ categories, most genes were involved in the ‘metabolic process’, ‘catalytic activity’ and ‘chloroplast’ GO groups, respectively.
Among significant GO terms, most up-regulated DEGs were classified in several biological processes, including ‘metabolic process’ and ‘response to stress’. Interestingly, the most significant enriched GO term was ‘response to stress’ according to the over-represented FDR for up-regulated DEGs. In addition, other biological processes, such as ‘metal ion transport’, ‘lipid transport’ and ‘ion transport’, were enriched. The significantly-enriched GO terms, including ‘catalytic activity’ and ‘oxidoreductase activity’, were the most abundant in the ‘molecular function’ category, while ‘vacuolar’ and ‘vacuolar membrane’ were the most abundant GO groups in the ‘cellular component’ category. These functions were more efficacious under salt stress.
For down-regulated DEGs, the GO-enriched biological processes, such as ‘metabolic process’ and ‘protein peptidyl-prolyl isomerization’, were the most represented groups. Moreover, GO enrichment analysis indicated that GO terms including ‘catalytic activity’ and ‘structural constituent of ribosome’ were the largest groups in the ‘molecular function’ category. In addition, several enzyme activity-related GO terms, such as ‘monooxygenase activity’ and ‘hydrolase activity, hydrolysing O-glycosyl compounds’, were also enriched in this category. In the category of ‘cellular component’, ‘chloroplast’ and ‘cell’ were the most abundant GO terms.
Validation of DEGs using qRT-PCR
To confirm the results of RNA sequencing, a total of six randomly-selected genes, i.e. pyrroline-5-carboxylate synthetase (TraesCS3A02G363700), serine–threonine protein kinase (TraesCS2D02G493700), aldehyde dehydrogenase family-3 member_1, chloroplastic-like (TraesCS5B02G210100), proline dehydrogenase_2_mitochondrial-like (TraesCS1D02G212400), putative calcium-binding protein_CML29 (TraesCS7B02G335200) and GDSL_esterase/lipase_At2g40250-like (TraesCS7D02G094900) were examined for RT-qPCR validation. The results of the qRT-PCR analysis showed that the expression levels of the control and salt-treated samples (up and down-regulated) were in agreement with those in the RNA-Seq, indicating that the two results validated one another (online Supplementary Fig. S2). Generally, the significant positive correlation (R 2 = 0.99) between the expression levels of selected DEGs in qRT-PCR and the results of the RNA sequencing analysis confirms the validity of the transcriptomic profiling data.
Discussion
Several studies have investigated the transcriptomic profiling of wheat under salt stress conditions imposed on the roots (Goyal et al., Reference Goyal, Amit, Singh, Mahato, Chand and Kanika2016; Han et al., Reference Han, Wang, Wei, Liang, Dai, Xia and Liu2018; Amirbakhtiar et al., Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019). However, to the best of our knowledge, only a few studies have been performed on transcriptomic profiling of wheat leaves under salt stress conditions (Luo et al., Reference Luo, Teng, Fang, Li, Li, Chu, Li and Zheng2019).
Assembly of raw sequence data yielded 137,508,542 clean reads after quality filtering. DEGs analysis is a statistical analysis based on negative binomial (NB) distributions, in which the expression level of a gene is determined based on the abundance of its transcript. The results showed that the number of common DEGs in the two approaches (i.e. STAR + DESeq2 and Hisat2 + DESeq2) was 17,897 (online Supplementary Table S3 and Fig. 1).
Previous studies have reported a different number of significant DEGs or differentially expressed unigenes in wheat, e.g. 2495 (Goyal et al., Reference Goyal, Amit, Singh, Mahato, Chand and Kanika2016), 5128 (Amirbakhtiar et al., Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019) and 4506 (Mansouri et al., Reference Mansouri, Naghavi, Alizadeh, Mohammadi-Nejad, Mousavi, Hosseini Salekdeh and Tada2018). Observing a greater number of DEGs in the present study probably because of the sudden severe salt stress (online Supplementary Table S3).
Investigating the genes associated with salt in the local salt-tolerant wheat landrace (i.e. Roshan) using the high-throughput NGS method in this study indicated that a portion of the DEGs was engaged in the synthesis of protein kinases, as noted by Ma et al. (Reference Ma, Gu, Liang, Zhang, Jin, Wang, Shen and Huang2016). Protein kinase genes that are widely present in plants play important roles in signal transduction and the response to stress (Ma et al., Reference Ma, Gu, Liang, Zhang, Jin, Wang, Shen and Huang2016; Amirbakhtiar et al., Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019; Chakradhar et al., Reference Chakradhar, Reddy, Chandrasekhar, Khan, Reddy, Ferrante and Khan2019). Key salinity-tolerant genes, such as polyamine oxidase and hormone-related genes, which were previously reported by Xiong et al. (Reference Xiong, Guo, Xie, Zhao, Gu, Zhao, Li and Liu2017), were also identified in our findings (online Supplementary Table S3). In addition, numerous genes have been identified in this study, some of which may be considered as novel genes that can be implicated in tolerating salinity in wheat. Asparagine synthetase is one of the genes identified among the up-regulated DEGs. Although it has been described as an unnecessary amino acid in animal studies, it converts aspartate and glutamine to asparagine and glutamate by consuming ATP, respectively. Asparagine is found in most organs, and it is very effective in the cellular response to stress conditions (Lomelino et al., Reference Lomelino, Andring, McKenna and Kilberg2017). The expression pattern (log2FoldChange = 7.24, Padj = 0) and the presence of this gene in two enriched pathways (i.e. the biosynthesis of secondary metabolites and alanine, aspartate and glutamate metabolism) suggest that it may be considered as one of the key genes for salt tolerance in wheat (online Supplementary Table S3).
Serpins, which are involved in the proteolysis process, have a significant effect on plant growth and development, stress responses, and defence mechanisms against pests and pathogens. In addition, plant serpins are more potent inhibitors than those found in mammals (Roberts and Hejgaard, Reference Roberts and Hejgaard2008). Five serpin genes have been found to be up-regulated in our work (online Supplementary Table S3). Moreover, actin-depolymerizing factors are a family of microfilament proteins, which are less frequent in mammalian genomes compared to plant genomes. They are responsible for plant growth and development through organizing filamentous actin, while they make the plant tolerant or sensitive to pathogens, depending on the types of these pathogens (Inada, Reference Inada2017). We identified three up-regulated genes, encoding actin-depolymerizing factors, which appear to play a role in salt tolerance. Organelle division is an important process in plants, where vital compounds, including dynamin-related proteins (DRPs), are involved. Two types of these proteins, i.e. DRP3A and DRP3B, have been implicated in the peroxisomal and mitochondrial division in Arabidopsis; however, DRP3A is the main component (Aung and Hu, Reference Aung and Hu2012). Based on the obtained results, the DRP3 gene was involved in the enriched pathway of amino acid metabolism (valine, leucine and isoleucine biosynthesis) for up-regulated genes, which is an essential metabolic pathway in the salinity resistance of wheat. Mother of FT and TFL1 (MFT) proteins, encoded by the MFT gene in wheat, takes part in different functions, including regulating germination and seed dormancy (Xi and Yu, Reference Xi and Yu2010; Chono et al., Reference Chono, Matsunaka, Seki, Fujita, Kiribuchi-Otobe, Oda, Kojima and Nakamura2015; Nakamura et al., Reference Nakamura, Makiko, Stehno, Holubec, Morishige, Pourkheirandish, Kanamori, Wu, Matsumoto and Komatsuda2015), and it interferes with the ABA signalling pathway in Arabidopsis (Xi and Yu, Reference Xi and Yu2010). Four genes related to this protein were observed in the results, and their expression levels were found to be up-regulated under the applied salinity stress.
GO and KEGG enrichment analyses
Approximately 30% of the significantly up-regulated DEGs in the ‘biological processes’ category were assigned into four function terms; namely, ‘metabolic process’, ‘response to stress’, ‘metal ion transport’, and ‘lipid transport’ (online Supplementary Table S9). Moreover, the GO analysis based on three main classes of ‘molecular functions’, ‘biological processes’ and ‘cellular component’ showed that the highest percentage of genes belonged to ‘molecular function’, while ~48% of the significantly up-regulated DEGs in this category belonged to ‘catalytic activity’ and ‘oxidoreductase activity’ (Ma et al., Reference Ma, Gu, Liang, Zhang, Jin, Wang, Shen and Huang2016; Zhang et al., Reference Zhang, Wang, Zhang, Dong, Chen and Cui2016). The results of the GO analysis in the present study showed that 12 and 28 GO biological process terms were significantly enriched in up-regulated and down-regulated genes, respectively (online Supplementary Table S9).
The results of KEGG enrichment analysis showed that some pathways, including ‘metabolic pathways’, ‘biosynthesis of secondary metabolites’, ‘carbon metabolism’, ‘plant hormone signal transduction’ and ‘starch and sucrose metabolism’, were the most significant pathways, which were associated with the most DEGs (online Supplementary Table S5). Plant hormones, known as growth regulators, actually stimulate the physiological response of the plants to adverse environmental conditions and stresses, such as salinity and drought (Kaya et al., Reference Kaya, Tuna, Yokaş, Ashraf, Ozturk and Athar2009). Therefore, the plant hormone signal transduction pathway, as a significantly enriched pathway of up-regulated genes, can be one of the most important to enhance salt tolerance in wheat. Signal transduction under stress conditions begins with understanding the signal, then, through the second messenger, it continues with phosphorylation of the proteins that are directly involved in the process of cell protection or TFs that control specific stress genes (Xiong et al., Reference Xiong, Schumaker and Zhu2002). Evaluating the analyses and their correlation with those of other studies focusing on the transcriptome analysis of the shoot using two resistant and susceptible wheat cultivars by Xiong et al. (Reference Xiong, Guo, Xie, Zhao, Gu, Zhao, Li and Liu2017) showed that ‘butanoate metabolism’ was identified as a new pathway in response to salinity stress, which was also found in transcriptome analysis of the local salt-tolerant wheat landrace (i.e. Roshan). In addition, their analysis showed that some pathways, such as ‘metabolic pathways’, ‘carbon metabolism’ and ‘plant hormone signal transduction’, were enriched. These three pathways were significantly enriched in the present study (online Supplementary Table S5). The last pathway was also identified by Amirbakhtiar et al. (Reference Amirbakhtiar, Ismaili, Ghaffari, Nazarian Firouzabadi and Shobbar2019), and it seems to have played an important role under salt stress.
Moreover, the results showed that ~88% of significant pathway genes (4963 out of 5634) were involved in metabolism pathways. This suggests that salt stress affects different aspects of the plant's metabolism, such as the metabolism of carbohydrates, lipids, secondary metabolites and amino acids. Furthermore, in the metabolism class, different sub-categories including carbohydrate metabolism (627), amino acid metabolism (251), biosynthesis of other secondary metabolites (191), metabolism of other amino acids (159), lipid metabolism (151), metabolism of terpenoids and polyketides (66), energy metabolism (55), nucleotide metabolism (48) and metabolism of cofactors and vitamins (48) were the most abundant. The sub-category of ‘carbohydrate metabolism’ was overrepresented since several up-regulated genes, involved in ‘galactose metabolism’, ‘pyruvate metabolism’, ‘starch and sucrose metabolism’, ‘glycolysis/Gluconeogenesis’, ‘glyoxylate and dicarboxylate metabolism’, ‘amino sugar and nucleotide sugar metabolism’, ‘ascorbate and aldarate metabolism’, ‘pentose and glucuronate interconversions’, ‘fructose and mannose metabolism’ and ‘inositol phosphate metabolism’, were significantly enriched. It should be noted that a similar conclusion was reached by Yong et al. (Reference Yong, Zou, Kok, Kwan, Chow, Nasu, Nanzyo, Kitashiba and Nishio2014) (online Supplementary Table S10). Overall, these findings are in accordance with those reported by Zhang et al. (Reference Zhang, Liu, Han, Song, Bai, Gao, Zhang, Yang, Li, Gao and Li2015).
The MapMan tool was used to enrich DEGs in the pathways related to salt stress conditions (Figs 4 and 5). In photosynthesis and the Calvin cycle, the adsorption and conversion of atmospheric carbon dioxide into biomass is especially important. Moreover, transketolase (TK) is one of the enzymes that have been reported to play an effective and key role in this mechanism (Hu et al., Reference Hu, Chang, Chen, Kuo-Huang, Liao, Huang and Juan2011). Although TK is believed to be the main enzyme in the Calvin cycle, it is also associated with the plant's resistance to abiotic stresses, including its role in cold tolerance by increasing carbon uptake capacity, reducing oxidative damage and stabilizing the cell structure (Bi et al., Reference Bi, Li, Wang and Ai2018). Previous studies on glycolate oxidase (GLO), which is known as a major enzyme involved in the photorespiratory metabolism, have also shown that the decreased activity of this enzyme reduces plant growth and inhibits photosynthesis in rice (Cui et al., Reference Cui, Ys, Li, Yang and Peng2016). Early light-inducible proteins (ELIPs), produced in the presence of light, are located in the thylakoid membrane, and they are responsible for protecting the plant's photosynthetic system against various environmental stresses, while they are also known for converting chloroplast to chromoplast (Timerbaev and Dolgov, Reference Timerbaev and Dolgov2019; Lee et al., Reference Lee, Lee, Han and Kim2020). Amino acids, which are the main constituents of proteins, affect plant growth and the production of metabolic energy (Goldfarb, Reference Goldfarb2020). As they contain nitrogen storage, they are essentially important in biological activities (O'Neill and Lee, Reference O'Neill and Lee2020). In general, free amino acids (e.g. valine, leucine, isoleucine, etc.) are known as major factors in protein metabolism, and their concentration is low compared to other proteins (Christensen, Reference Christensen1964). Mannan is involved in the early stages of plant growth and seed development as a source of carbohydrates, and it is synthesized in Arabidopsis by cellulose synthase-like (CSL) proteins as well (Verhertbruggen et al., Reference Verhertbruggen, Yin, Oikawa and Scheller2011; Wang et al., Reference Wang, Mortimer, Davis, Dupree and Keegstra2012). With regard to tethering proteins and their role in plants, several functions have been identified, including plant growth, collaboration in plant adaptation, acting as a connection mediator in donor and receiver membranes and indirect response to abiotic stresses; however, some of them are direct targets of pathogens in response to biotic stress (Ravikumar et al., Reference Ravikumar, Steiner and Assaad2017). Structural proteins are particularly important for plants since they are a crucial source of energy and carbon for seed germination and plant growth and development. Moreover, they are also responsible for stress responses, lipid metabolism and hormone signalling. Caleosins and steroleosins belong to the group of structural proteins. An important role of caleosins involves organelle dynamics and stabilization, while steroleosins are involved in metabolism and signalling of steroid hormone (Shao et al., Reference Shao, Liu, Su, Ma and Wang2019). Obtusifoliol 14-alpha demethylase is a member of the cytochrome P450 gene family, which catalyses the formation of steroid hormones involved in lipid signalling in plants (O'Brien et al., Reference O'Brien, Chantha, Rahier and Matton2005). Glycerol-3-phosphate dehydrogenase (GPDH) plays an important role in glycerolipid metabolism, plant growth and the response to abiotic stresses. Increased levels of GPDH accumulation in roots and early stages of Zea mays seed growth have also been demonstrated in response to abiotic stresses, including salinity (Zhao et al., Reference Zhao, Li, Wang, Zhao, Gao, Zhao, He, Li and Xu2018). Asparaginyl endopeptidases (AEPs) catalyse the biosynthesis of peptides, while they also contribute to plant defence and grain protein storage (Jackson et al., Reference Jackson, Gilding, Shafee, Harris, Kaas, Poon, Yap, Jia, Guarino, Chan, Durek, Anderson and Craik2018; Du et al., Reference Du, Yap, Chan, Rehm, Looi, Poth, Gilding, Kaas, Durek and Craik2020). Elimination of damaged proteins and repairing plant damage, caused by exposure to light levels higher than required for photosynthesis, are the functions of the protein quality control mechanisms in plants (Kato and Sakamoto, Reference Kato and Sakamoto2018). FtsH protease complexes play crucial roles in thylakoid membrane biogenesis and protein quality control mechanisms in photosystem II (Kato and Sakamoto, Reference Kato and Sakamoto2018). The role of the SnRK1 subfamily has been demonstrated in a number of processes, such as regulating plant growth and development by the interaction between protein kinase 1 (SnRK1) and abscisic acid (ABA); controlling metabolisms such as carbohydrate metabolism, regulation of cellular energy (balance and homeostasis), inhibition of energy consumption in anabolic pathways, and enhanced catabolism; contributing to stress tolerance through direct phosphorylation of metabolic enzymes; and responding to stresses that deplete plant ATP (Mohannath et al., Reference Mohannath, Jackel, Lee, Buchmann, Wang, Patil, Adams and Bisaro2014; Crepin and Rolland, Reference Crepin and Rolland2019; Krasnoperova et al., Reference Krasnoperova, Buy, Goriunova, Isayenkov, Karpov, Blume and Yemets2019; Belda-Palazón et al., Reference Belda-Palazón, Adamo, Valerio, Ferreira, Confraria, Reis-Barata, Rodrigues, Meyer, Rodriguez and Baena-González2020).
Transcription factors
TFs play major roles in the response to different abiotic stresses. They are master regulators of abiotic stress responses in plants (Lindemose et al., Reference Lindemose, O'Shea, Jensen and Skriver2013). A large number of studies have been published on the roles of TFs families in plant responses to different stress conditions, including GRAS (Bolle, Reference Bolle2004), C2H2 zinc finger (Wang et al., Reference Wang, Ding, Cai, Chen and Zhu2018a, Reference Wang, Liu, Wu, Li, Wang, Cui and Zhuang2018b), CAMTAs (Finkler et al., Reference Finkler, Ashery-Padan and Fromm2007), basic helix-loop-helix (bHLH) (Kim and Kim, Reference Kim and Kim2006), the APETALA2/ethylene-responsive factor, AP2/ERF (Mizoi et al., Reference Mizoi, Shinozaki and Yamaguchi-Shinozaki2012), basic leucine zipper (bZIP) (Zhu et al., Reference Zhu, Meng, Cai, Li, Dong and Li2018), calmodulin-binding transcription activators, NAM/ATAF1/CUC2 and NAC (Nuruzzaman et al., Reference Nuruzzaman, Sharoni and Kikuchi2013), and MYB (Dubos et al., Reference Dubos, Stracke, Grotewold, Weisshaar, Martin and Lepiniec2010). These TFs families play significant roles in translating abiotic stress signals into changes in gene expression (Lindemose et al., Reference Lindemose, O'Shea, Jensen and Skriver2013). A GRAS protein was the most abundant TF in the present study. Different roles of GRAS proteins, such as plant signal transduction, development (Bolle, Reference Bolle2004; Hirsch and Oldroyd, Reference Hirsch and Oldroyd2009; Sun et al., Reference Sun, Jones and Rikkerink2012), plant growth regulation, responses to multiple stresses (Wang et al., Reference Wang, Ding, Cai, Chen and Zhu2018a, Reference Wang, Liu, Wu, Li, Wang, Cui and Zhuang2018b; Chen et al., Reference Chen, Zhu, Wu, Lu, Sun, Cao, Li and Xu2019; Wang et al., Reference Wang, Wang, Liu, Zhang, Wang, Li and Gao2019; Li et al., Reference Li, Dong, Liu, Bo, Miao, Beckles, Zhang and Gu2020) and adaptation to unfavourable environmental conditions (Chen et al., Reference Chen, Zhu, Wu, Lu, Sun, Cao, Li and Xu2019), have been identified. Furthermore, transcriptome analysis of Roshan revealed 111 down-regulated and 65 up-regulated genes for the CAMTA TF family, which is one of the fast response stress proteins (Büyük et al., Reference Büyük, İlhan, Şener, Özsoy and Aras2019). Plant CAMTAs (calmodulin-binding transcription activators) integrate stress and growth signals, such as increased Ca2+ concentration (Finkler et al., Reference Finkler, Ashery-Padan and Fromm2007). Bailey et al. (Reference Bailey, Martin, Toledo-Ortiz, Quail, Huq, Heim, Jakoby, Werber and Weisshaar2003) identified 162 BHLH genes in the Arabidopsis thaliana genome. They report that this family is one of the largest transcription factor gene families in A. thaliana. Another family of TF, which was identified in the present study, was AP2/ERF. This family has four major subfamilies (Sakuma et al., Reference Sakuma, Maruyama, Osakabe, Qin, Seki, Shinozaki and Shinozaki2006), i.e. DREBs and ERFs, which are involved in plant abiotic stress responses (Mizoi et al., Reference Mizoi, Shinozaki and Yamaguchi-Shinozaki2012); AP2, which is involved in floral organ identity (Dinh et al., Reference Dinh, Girke, Liu, Yant, Schmid and Chen2012), abiotic stress and plant development; and RAV, which is involved in leaf senescence (Woo et al., Reference Woo, Kim, Kim, Kim, Lee, Song, Kim, Lee, Nam and Lim2010) and seed germination (Feng et al., Reference Feng, Chen, Wang, Kong, Wu and Chen2014). As noted earlier, two classes of transcription factors include basic leucine zipper (bZIP) and WRKY. The bZIP genes act as crucial regulators in ABA-mediated stress response in plants (Zhu et al., Reference Zhu, Meng, Cai, Li, Dong and Li2018), while the ClWRKY genes play significant roles in the growth and development of various tissue in Citrullus lanatus (Yang et al., Reference Yang, Li, Yang, Wang, Mo, Zhang, Zhang, Ma, Wei and Zhang2018), signalling pathways and regulatory networks (Chen et al., Reference Chen, Song, Li, Zhang, Zou and Yu2012) and plant pathogen responses (Lindemose et al., Reference Lindemose, O'Shea, Jensen and Skriver2013). Recently, researchers have shown the effects of abiotic stress on the WRKY TFs family. For instance, among WRKY family members in rice, WRKY50 and WRKY72 were reported to be differentially regulated under salt stress (Çelik et al., Reference Çelik, Meriç, Ayan and Atak2019). Moreover, 171 TaWRKY TFs were identified from the entire wheat genome, most of which were enriched on four chromosomes, especially on chromosome 3B (Ning et al., Reference Ning, Liu, Kang and Lv2017). Another TF family identified in the present study was NAC. NAC proteins play important roles in development, abiotic and biotic stress responses and biosynthesis (Welner et al., Reference Welner, Deeba, Leggio, Skriver and Gonzalez2016).
DEGs involved in transportation
Based on many studies, different kinds of genes are involved in the response to salt stress. The analysis showed that there were 615 significant DEGs (including up-regulated and down-regulated genes) involved in different transportation mechanisms. The plasma membrane Na+/H+ antiporter (also called exchanger) Salt Overly Sensitive1 (SOS1) plays a major role by controlling Na+ homeostasis, and possibly contributing to the sensing of sodicity stress (EL Mahi et al., Reference EL Mahi, Hormaeche, Luca, Villalta, Espartero, Arjona, Fernandez, Bundo, Mendoza, Mieulet, Lalanne, Lee, Yun, Guiderdoni, Aguilar, Leidi, Pardo and Quintero2019). In the present study, this gene was identified among the significant up-regulated DEGs. Another plasma membrane transporter for removing calcium (Ca2+) from the cell is calcium-transporting ATPases. In the present study, 13 DEGs for calcium-transporting ATPases, including 12 up-regulated and 1 down-regulated genes, were detected. We also identified other significant DEGs involved in cell transportation, such as the ABC transporter (115 DEGs including 68 up and 47 down-regulated), HKT genes, potassium transporter 10-like, sucrose transporter, cation-chloride cotransporter, cationic amino acid transporter, phospholipid-transporting ATPase, probable copper-transporting ATPase, manganese-transporting ATPase, putative inactive cadmium/zinc-transporting ATPase, copper-transporting ATPase, chloroplastic-like, Annexin (10 DEGs including 4 up and 6 down-regulated), and many others. The most significant up-regulated DEGs (TraesCS3A02G217900) with padj = 7.42 × 10−106 belonged to the ABC transporter family members. In the case of Na+ transporter genes, nine HKT genes (i.e. HKT1, HKT2, HKT4, HKT7 and HKT8) were identified in the local salt-tolerant wheat landrace (i.e. Roshan), of which HKT8 was the most significant up-regulated DEG (online Supplementary Table S11).
DEGs related to the synthesis of late embryogenesis-abundant (LEA) and heat-shock proteins (Hsps)
Another important gene family involved in the reduction of salinity stress includes the LEA proteins. Liu et al. (Reference Liu, Xing, Yang, Mu, Wang, Lu, Wang and Zhang2019) identified a total of 179 LEA genes in T. aestivum under different abiotic stresses. Comparing DEGs in the present study with the NCBI databases and the previously published report (Liu et al., Reference Liu, Xing, Yang, Mu, Wang, Lu, Wang and Zhang2019) led to the identification of 99 LEA genes (online Supplementary Table S12). Moreover, we discovered 51 Hsp genes (35 up and 16 down-regulated), including Hsp40, Hsp70, Hsp90 and Hsp101, under salinity stress in wheat (online Supplementary Table S13). Several studies have investigated the effects of abiotic stress on the synthesis of LEA and Hsps proteins (Wang et al., Reference Wang, Vinocur, Shoseyov and Altman2004; Al-Whaibi, Reference Al-Whaibi M2011; Afzal et al., Reference Afzal, Howton, Sun and Mukhtar2016; Liu et al., Reference Liu, Xing, Yang, Mu, Wang, Lu, Wang and Zhang2019).
Conclusion
In conclusion, the leaf transcriptome analysis of an Iranian local salt-tolerant wheat landrace (i.e. Roshan) indicates that different genes and pathways are involved in achieving tolerance against salt stress. Our findings will enrich existing genomic resources and provide incentives for research into improving salinity tolerance in important wheat species. In addition to being helpful in genomic analysis, the data generated by our study can be effective in identifying genes and analysing their expression, as well as initiating functional and comparative genomic studies. These results will also be a rich and valuable resource for further analysis of salt tolerance and breeding wheat cultivars through the use of genes related to salt stress.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S1479262121000319
Acknowledgements
The authors wish to thank Mrs Azarivand (NIGEB Genome Service Center Laboratory) for carrying out the RNA extractions, RNeasy Plant Mini Kit QIAGEN processing. We would also like to express our appreciation to Mrs Loni for her valuable suggestions.
Financial support
This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.
Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Ethical standards
This article does not contain any studies with human participants or animals performed by any of the authors.
Availability of data and material
The datasets supporting this article have been uploaded as part of the supporting information.
