Transcriptomic analysis of differentially expressed genes in leaves and roots of two alfalfa (Medicago sativa L.) cultivars with different salt tolerance

Background Alfalfa (Medicago sativa L.) production decreases under salt stress. Identification of genes associated with salt tolerance in alfalfa is essential for the development of molecular markers used for breeding and genetic improvement. Result An RNA-Seq technique was applied to identify the differentially expressed genes (DEGs) associated with salt stress in two alfalfa cultivars: salt tolerant ‘Halo’ and salt intolerant ‘Vernal’. Leaf and root tissues were sampled for RNA extraction at 0 h, 3 h, and 27 h under 12 dS m− 1 salt stress maintained by NaCl. The sequencing generated a total of 381 million clean sequence reads and 84.8% were mapped on to the alfalfa reference genome. A total of 237 DEGs were identified in leaves and 295 DEGs in roots of the two alfalfa cultivars. In leaf tissue, the two cultivars had a similar number of DEGs at 3 h and 27 h of salt stress, with 31 and 49 DEGs for ‘Halo’, 34 and 50 for ‘Vernal’, respectively. In root tissue, ‘Halo’ maintained 55 and 56 DEGs at 3 h and 27 h, respectively, while the number of DEGs decreased from 42 to 10 for ‘Vernal’. This differential expression pattern highlights different genetic responses of the two cultivars to salt stress at different time points. Interestingly, 28 (leaf) and 31 (root) salt responsive candidate genes were highly expressed in ‘Halo’ compared to ‘Vernal’ under salt stress, of which 13 candidate genes were common for leaf and root tissues. About 60% of DEGs were assigned to known gene ontology (GO) categories. The genes were involved in transmembrane protein function, photosynthesis, carbohydrate metabolism, defense against oxidative damage, cell wall modification and protection against lipid peroxidation. Ion binding was found to be a key molecular activity for salt tolerance in alfalfa under salt stress. Conclusion The identified DEGs are significant for understanding the genetic basis of salt tolerance in alfalfa. The generated genomic information is useful for molecular marker development for alfalfa genetic improvement for salt tolerance. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-021-03201-4.


Background
Alfalfa (Medicago sativa L.) is an important forage legume in the world. Cultivated alfalfa is an outcrossing autotetraploid (2n = 4x = 32) with a genome size of 800-1000 Mb [1]. Although alfalfa is regarded as moderately tolerant to salinity [2], alfalfa yield reduces by approximately 6-7% for each dS m − 1 increase above a salinity of 2 dS m − 1 [3]. To stabilize alfalfa production under saline regions, the development of superior salt tolerant cultivars becomes an important breeding goal. Identification of candidate genes for salt tolerance can increase the accuracy of parental selection as this trait has low heritability [4]. Salt tolerance is a complex trait controlled by multiple genes, involving different signaling pathways, osmotic tolerance, ion transport, compartmentalization of salt ions in vacuoles, the synthesis of plant hormones and photosynthesis [5].
Next-generation sequencing technologies have been used to identify candidate genes involved in salt tolerance of alfalfa. Transcriptomic studies in the 1-week old root tissue of alfalfa under salt stress found 1165 DEGs, including 86 transcription factors, which are responsible for stress tolerance, kinase, hydrolase, and oxidoreductase activities [6]. Luo et al. [7] identified 8861 DEGs in 12-day old seedlings of alfalfa under salt stress, which are responsible for ion homeostasis, antiporter, signal perception, signal transduction, transcriptional regulation, and antioxidative defense. Lei et al. [8] revealed 2237 DEGs between salt tolerant and intolerant alfalfa cultivars and found a salt tolerant alfalfa cultivar maintained relatively stable expression of genes responsible for reactive oxygen species and Ca 2+ pathway, phytohormone biosynthesis and Na + /K + transport under stress. Gruber et al. [9], using bulked genotypes as replications, studied transcriptomes in alfalfa and found genes responsible for numerous functions in a salt intolerant alfalfa cultivar. In recent years, genetic modification of certain genes controlling salt tolerance have also been conducted in alfalfa. Overexpression of salt responsive genes or transcription factors had improved salt tolerance in transgenic alfalfa. Such genes include Alfin1 [10], AVP1 [11], GmDREB1 [12], SsNHX1 [13], TaNHX2 [14], GsCBRLK [15], GsZFP1 [16], OsAPX2 [17], SeNHX1 [18], AtNDPK2 [19], AgcodA [20], and GsWRKY20 [21]. These studies have advanced our understanding of the genetic control for salt tolerance in alfalfa. However, most studies mainly focused on single time point sampling of root tissue at the seedling stage after salt stress, limiting the analysis of the temporal expression of genes affecting salt tolerance.
Tissue specific protein induction is regulated during salinity stress and is unique to roots and shoots [22]. Thus, there should be tissue specific transcriptomic responses [23][24][25]. Although the root is the first receptor of salt stress [6,7], leaf tissue is the main energy source for plant growth and stress tolerance during active growth and developmental stages. To advance our knowledge about the temporal gene expression in different tissues for the genetic control under salt stress between tolerant and intolerant cultivars, we conducted a RNA-Seq analysis with the objective to simultaneously analyze gene expressions of leaf and root tissues of two alfalfa cultivars with different tolerance to salinity after exposing them to 12 dS m − 1 of electrical conductivity salt stress for 0 h, 3 h, and 27 h. The analysis was fruitful with the identification of many unique genes conditioning salt tolerance in alfalfa.

High throughput sequencing and assembly
A total of 408 million raw sequence reads were generated using the Illumina HiSeq sequencing platform. The reads were reduced to 93.5% (381 million clean reads) by removing adapter contamination and reads with length lower than 36 bp (Table 1). There were 84.8% of clean reads mapped to the alfalfa reference genome using STAR (v2.6.1a). The samples showed high percentages (78.8-92.4%) of mapping with the alfalfa reference genome except for 'Halo' root tissue sampled at 27 h of salt stress.

Differentially expressed genes (DEGs)
In leaf tissue, there were 237 DEGs identified between the two alfalfa cultivars. Among them, 34 DEGs were expressed at all three time points (0 h, 3 h, and 27 h) and 17 DEGs expressed after exposing to salt stress (3 h and 27 h) (Fig. 1a, b; Additional file 1: Table S1). Of these DEGs, 39, 31, and 49 DEGs were specific to 'Halo', and 34, 34, and 50 DEGs were specific to 'Vernal' at 0 h, 3 h and 27 h of salt stress, respectively (Fig. 1b). The number of DEGs in leaf tissue decreased after 3 h compared to 0 h treatment. Then, the number of DEGs increased from 3 h to 27 h of salt stress for both cultivars (Fig. 1b).
In root tissue, a total of 295 DEGs were identified between the two alfalfa cultivars. There were 33 DEGs expressed at all three time points and 5 DEGs expressed after exposing to salt stress (Fig. 2a, b; Additional file 1: Table S2). Of these DEGs, 68, 55, and 56 DEGs were specific to 'Halo' at 0 h, 3 h and 27 h, whereas 64, 42, and 10 DEGs were specific to 'Vernal', respectively (Fig.  2b). The number of DEGs in root tissue decreased at 3 h as compared to 0 h treatment for both cultivars, but the decrease was greater for 'Vernal' than for 'Halo'. The main difference of DEGs between the two cultivars in the root was from 3 h to 27 h, with a 76% decrease in DEGs in 'Vernal' while there was almost no change for 'Halo' (Fig. 2b). After 27 h of salt stress in root tissue, the number of DEGs in 'Halo' were five times more than that of 'Vernal' (Fig. 2b).
To identify pathways involved in salt tolerance, we carried out Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis of the DEGs. In total, 64 (27%) DEGs from leaf tissue and 86 (29.15%) DEGs from root tissue were assigned to 65 KEGG pathways ( Table  2). In both tissues, the most significant DEGs were represented in the pathways of metabolism and biosynthesis of secondary metabolites. Of these, five pathways were common among different time points and alfalfa tissues. The highest level of enriched DEGs were in 14 pathways in leaf tissue and 6 pathways in root tissue after 27 h of salt stress. Among these pathways, the three highest enriched DEGs were involved in plant hormone signal transduction.

Candidate genes to enhance salt tolerance in alfalfa
The detected DEGs can be classified into two major groups for the candidate genes responsible for salt tolerance in alfalfa: 1) genes consistently expressed under short-term and long-term salt stress (3 h and 27 h) in 'Halo', and 2) the genes consistently expressed at all three time points in 'Halo'. In the first group, there were 13 genes (11 in leaf; 2 in root) consistently expressed at both 3 h and 27 h of salt stress. While in the second group, there were 46 genes (17 in leaf, 29 in root) consistently expressed at all three time points. Thirteen candidate genes were highly expressed in both leaf and root tissues of 'Halo' as compared to 'Vernal', while 15 and 18 candidate genes revealed tissue specific expression in the leaf and root tissues of 'Halo', respectively (Tables 3,  4, and 5). Among the genes expressed in both tissues, MS.gene029203 (F-box/LRR-repeat protein 4) showed increasing expression with time in both leaf and root tissues of 'Halo', while MS.gene049294 (caffeic acid 3-Omethyltransferase) showed increasing expression with time in leaf tissue and MS.gene01091 (T-complex protein 1 subunit gamma) and MS.gene32989 (hypothetical protein TSUD_06780) showed increasing expression with time only in root tissue. Among the genes with leaf tissue specific expression, MS.gene029201 (replication protein A 70 kDa DNA-binding subunit C), MS.gene029206 (FAD synthetase 1, chloroplastic), and MS.gene24098 (thioredoxin-like protein CDSP32 chloroplastic-like) showed increasing expression with time. Among the genes with root tissue specific expression, MS.gene011517 (14 kDa proline-rich protein DC2.15) and MS.gene013923 (histone lysine Nmethyltransferase, H3 lysine-9 specific SUVH1), had higher and consistent expression with time under salt stress.
In addition, there were also genes consistently expressed under salt stress in leaf (Additional file 1: Table S1) and root (Additional file 1: Table S2)

Discussion
This study generated a unique set of differentially expressed genes associated with salt tolerance in alfalfa. This finding is not only significant for understanding the temporal expression of genes conditioning salt tolerance in alfalfa, but also can be used to characterize alfalfa breeding material and develop molecular markers for salt tolerance selection. First, 84.8% of sequence reads were mapped onto the alfalfa reference genome. Secondly, 237 DEGs in leaf and 295 DEGs in root tissues were identified between the two alfalfa cultivars. Third, this study was able to determine candidate genes   Due to polyploidy and its out-crossing nature, alfalfa has encountered many challenges in genomic studies [26] as compared to self-pollinated crops such as wheat [27] and soybean (Glycine max) [28]. Thus, this transcriptomic study had considered several factors to overcome certain technical difficulties. First, identical clones were sampled at different time points from different alfalfa tissues. Unlike previous alfalfa transcriptomic studies, two technical replicates for each treatment were included to minimize technical errors. Furthermore, this study focused on both leaf and root tissues of alfalfa cultivars to capture tissue specific gene expression. These considerations seemed to be effective in capturing about 381 million high quality reads, which likely represents most of the genome of M. sativa. The raw reads showed a high percentage of mapping with the reference genome. These outputs should have enhanced our detection  [29]. One of the main differences between the two cultivars was the number of DEGs in roots. In the root of salt tolerant alfalfa, the number of DEGs was similar between 3 h and 27 h of salt stress, but a sharp decrease was observed between 3 h and 27 h in the intolerant cultivar 'Vernal'. We speculate that such earlier activation of salt responsive genes and maintenance of a large number of DEGs might be a key characteristic for salt tolerance in alfalfa, suggesting alfalfa tolerance is associated with upregulation of key genes from short term salt stress. About 60% of DEGs were assigned to GO categories, while KEGG pathways for less than 30% DEGs were identified in this study. The DEGs were mainly involved in metabolic pathways as revealed by KEGG pathway analysis. Although certain pathways involved in salt tolerance may be conserved in plant species such as in halophytes, there was still variation among plant species, cultivars, and tissues [5]. This study demonstrated that transcriptional variation in adaptation to salt stress exists not only among the alfalfa cultivars but also between the different tissues. 'Ion binding' (GO:0043167) was significantly enriched in both leaf and root tissues of 'Halo', but not in 'Vernal' under salt stress. This suggested that the genes responsible for 'ion binding' should be unique for salt tolerance of 'Halo' alfalfa. Therefore, the tissueand genotype-specific salt responsive genes might be useful in identification of salt tolerant genotypes in the future.
Among 13 candidate genes expressed in leaf and root tissues of 'Halo' under salt stress (Table 3) in this study, two genes (MS.gene013222 and MS.gene52595) are responsible for transmembrane protein function. These transmembrane proteins control gateways and selective transport of salt ions to facilitate salt tolerance in plants. Likewise, MS.gene013211, a homologous gene to ribonuclease TUDOR1, is involved in stress adaptation and highly expressed in leaf and root tissues of 'Halo' in our study [30]. MS.gene93979, a homologous gene to NF-X1-type zinc finger protein, is part of mechanisms that regulate growth under salt stress and was highly expressed in leaf and root tissues of 'Halo' in our study [31]. In addition, MS.gene029202 (E3 ubiquitin-protein ligase CIP8), MS.gene029203 (F-box/ LRR-repeat protein 4), MS.gene36780 and MS.gene36960 (elongation factor 1-alpha) were highly expressed in leaf and root tissues of salt tolerant alfalfa in our study. These genes are involved in regulation of a number of biological processes including biotic and abiotic stress tolerances [32][33][34]. For example, MS.gene049294, which is a homologous gene of Omethyltransferase, was found to improve salt tolerance in transgenic Arabidopsis [35]. MS.gene01091, a homologous gene to the T-complex protein 1 subunit gamma, showed high expression in both root and leaf tissue and is involved in intracellular assembly and folding of various proteins [36]. MS.gene029200, a homologous gene to replication factor A protein, was highly expressed in both leaf and root tissues of 'Halo' in our study, which might play a role in binding, replication, repair, and recombination of DNA under stress conditions [37].
In this study, we found 15 and 18 candidate genes specific to leaf and root tissues of salt tolerance 'Halo' alfalfa (Tables 4,5). In leaf tissue, nine genes showed consistent expression under salt stress, while six of them were expressed at all three time points. In our study, salt tolerant alfalfa showed an increased expression of MS.gene024018 and MS.gene24098 with putative functions of chloroplastic glutaredoxin and thioredoxin-like protein CDSP32, respectively. The two genes (MS.gene024018, MS.gene24098) were found to be important for defense against protein oxidative damage in other studies [38,39]. This is important because salt stress results in the formation of reactive oxygen species, which damage protein, membrane lipids, and nucleic acids [40]. MS.gene63155, a homologous gene to receptor-like kinases (RLKs), are a family of transmembrane proteins, showed lowered expression with time under salt stress. This gene is involved in plant growth as well as stress response [41]. Two nucleoporin proteins (MS.gene037960 and MS.gene39381) were expressed consistently under salt stress in leaf tissue of 'Halo'. These proteins connect cytoplasm and nucleoplasm, and are involved in abiotic stress tolerance [42]. MS.gene038586 is a homologous gene to kinesin super family proteins which plays a significant role in intracellular transport and are critical for cellular functioning and survival [43]. MS.gene029206, a homologous gene to FAD synthetase 1, is a co-factor for various enzymes that participate in numerous metabolic processes like photosynthesis, electron transport, fatty acid oxidation and biosynthesis of secondary metabolites [44]. MS.gene36621, a homologous gene to stem 28 kDa glycoprotein, which is known as a vegetative storage protein, was highly expressed under salt stress in our study. This protein plays a certain role as a somatic storage protein during early seedling development [45]. Salt tolerant alfalfa showed a high expression of MS.gene07287 in leaf, a homologous gene to calvin cycle protein CP12-2. This gene is involved in photosynthesis and improved plant growth [46]. The Calvin-Benson cycle is the primary pathway of carbon fixation, producing carbon compounds. CP12 facilitates the formation of a complex between glyceraldehyde-3-phosphate dehydrogenase and phosphoribulokinase, thereby increasing the photosynthetic capacity of the plant [46]. MS.gene99197, a homologous gene to zeaxanthin epoxidase (ZEP), was highly expressed at all three time points in leaf tissue of 'Halo'.
It is an important enzyme in ABA biosynthesis and plays an important role in osmotic tolerance [47].
In root tissue, one (MS.gene61130) of the 18 genes detected was consistently expressed under salt stress, while the rest were expressed at all three time points in our study (Table 5). MS.gene002389, a homologous gene to secretory carrier-associated membrane proteins, is involved in membrane trafficking and found to influence accumulation of secondary cell wall components in Poplus [48]. MS.gene011517, a homologous gene to 14-kDa proline-rich protein DC2.15, is involved in cell wall modification and organization [49]. The plant cell wall is not only a physical barrier between the plant and the environment but also is a responsive part of the plant to biotic and abiotic stresses. The finding of tissue specific salt tolerant candidate genes responsible for the plant cell wall is promising and underlines the need for further research on its role in response to salt stress. In addition, salt stress causes lipid peroxidation, resulting in damage of membrane lipids and eventual cell leakage. This study showed salt tolerant alfalfa had an increased expression of MS.gene049130, a homologous gene to aldehyde dehydrogenase, responsible for oxidation of aldehydes produced during lipid peroxidation thereby detoxifying cells [50]. MS.gene95536 is a homologous gene to acyl-CoAbinding domain-containing protein 6, which is associated with phospholipid metabolism. This gene also was shown to play a role in the freezing tolerance of Arabidopsis [51]. MS.gene070486, a homologous gene to phosphatidylinositol transfer proteins, plays an important role in signal transduction and facilitates lipid transfer between membranes [52]. MS.gene056386, a homologous gene to fructokinases, are important enzymes catalyzing fructose phosphorylation and are involved in plant growth and development [53]. MS.gene058673, a homologous gene to heavy-metalassociated domain-containing protein conferring tolerance to abiotic stress [54]. MS.gene073760, a homologous gene to probable E3 ubiquitin-protein ligase LOG2, which induces amino acid secretion. This is the main form of organic nitrogen in the plant [55]. MS.gene02427, a homologous gene to soluble inorganic pyrophosphatase, is tightly linked with carbohydrate metabolism. It plays an important role in stress adaptive responses [56]. Carbohydrate metabolism produces soluble carbohydrates that are important for salt tolerance because of its osmotic adjustment function in the root.

Conclusion
Our study generated a unique set of DEGs for alfalfa salt tolerance studies and breeding efforts. The information is useful for better understanding of temporal expression of genes in response to salt stress. Furthermore, GO annotation and KEGG pathway analysis of the DEGs provided insights to the different molecular and biological processes between salt tolerant and intolerant alfalfa cultivars. In particular, 'ion binding activity' was found as a key molecular activity specific to salt tolerant alfalfa cultivar 'Halo'. Based on this finding, salt tolerance in alfalfa appears to be associated with consistent expression of genes for selective transport of salt ions and compounds, increasing photosynthetic capacity as well as carbohydrate metabolism, enhancing defense against oxidative damage, modification of root cell wall and protection against lipid peroxidation. The SNPs discovered in this study will be valuable for molecular marker-assisted breeding for the development of salt tolerant alfalfa.

Plant material and salt treatment
Two alfalfa cultivars, 'Halo' (obtained from Agriculture and Agri-Food Canada, Swift Current Research and Development Centre) and 'Vernal' (sourced from Dr. Biligetu's lab, Crop Development Centre, University of Saskatchewan) were chosen for the study. Cultivar 'Halo' was selected for improved salinity tolerance for germination, seedling growth, and mature plant regrowth at 100 mM NaCl in the greenhouse conditions [57], and cultivar 'Vernal' was considered as a salinity intolerant cultivar [58,59]. Four genotypes (biological replicates) of each cultivar were grown from seeds in the College of Agriculture and Bioresources greenhouse at the University of Saskatchewan (45 Innovation Blvd., Saskatoon, SK) for 12 weeks. Six identical clones of each biological replicate were produced by stem cuttings. Salt stress of 120 mM NaCl approximately corresponding to 12 dS m − 1 electrical conductivity was applied on 4 week old seedlings. Salt stress of 12 dS m − 1 was selected from our earlier greenhouse study where alfalfa was grown at various gradients of salt stress and alfalfa cultivars showed variation in response to salt stress at 12 dS m − 1 , with increase in salt stress from 12 dS m − 1 all alfalfa cultivars showed very high mortality (Bhattarai et al., unpublished). Leaf and root samples were collected immediately before salt treatment (control, 0 h), and at 3 h and 27 h of salt treatments. The samples were immediately frozen in liquid nitrogen and then stored at − 80°C for 2 weeks until total RNA extraction carried out.

Tissue sample and RNA isolation
About 100 mg of tissue samples were disrupted using TissueLyser II and total RNA was extracted with RLT buffer using the Qiagen RNeasy Plant Mini Kit (Qiagen Inc., Mississauga, ON, Canada) according to the manufacturer's protocol. DNase treatment was performed using the Ambion DNA-free DNase treatment and removal reagents (Life Technologies, Carlsbad, CA, USA) to remove contaminant genomic DNA from the isolated total RNA. Nanodrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA) was used to measure the total RNA concentration. RNA integrity number was evaluated for 12 samples using RNA 6000 Nano labchip on 2100 Agilent Bioanalyzer (Agilent Technologies, Waldbronn, Germany) (Additional file 1: Table S3; Additional file 2: Fig. S1).

Library preparation and sequencing
Poly (A) RNA was purified from total RNA using Magnosphere MS150 OligodT beads according to the manufacturer's protocol. Reference-based mapping, differential gene expression analysis and annotation The quality of the raw sequence was assessed using the FastQC software [60]. The raw reads were cleaned by removing adapters and low-quality sequences using Trimmomatic v.0.36 based on the default setting of pairedend mode, phred 33 and threads 6 [61]. The trimmed high-quality reads of samples from the two technical replicates were merged and mapped with the alfalfa reference genome (https://doi.org/10.6084/m9.figshare. 12327602.v3) [62,63] using STAR (v2.6.1a) [64] with "quantMode" as "GeneCounts". The obtained "ReadsPer-Gene" of each sample were extracted as count matrix and the differentially expressed genes were analyzed using DeSeq2 package [65] where data were normalized by the median of the ratios. The threshold of padj < 0.001 and the Log fold change (Log2FC) > 2 were used to determine the significance of gene expression differences. The functional annotation of the DEGs were also extracted via searches of NR databases as available in "query.blastp.db.out" and gene ontologies were obtained via searches of the GO databases as available in "Msa.-GO.list.up" likewise Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog (KO) were obtained via search of KO databases as available in "query.ko" from Zeng [63]. Gene ontology analysis of the DEGs was done for biological process, cellular components, and molecular function by AgriGO v2.0 software [66]. Venn diagrams were produced using the Venny tool [67].

Identification of single nucleotide polymorphisms (SNPs)
SNPs calling was done using freebayes software using the bam file generated in the mapping process where at least 5 supporting observations were required to be