Transcriptomic profiling and analysis of differentially expressed genes in asparagus bean (Vigna unguiculata ssp. sesquipedalis) under salt stress

Asparagus bean (Vigna unguiculata ssp. sesquipedalis) is a warm season legume which is widely distributed over subtropical regions and semiarid areas. It is mainly grown as a significant protein source in developing countries. Salinity, as one of the main abiotic stress factors, constrains the normal growth and yield of asparagus bean. This study used two cultivars (a salt-sensitive genotype and a salt-tolerant genotype) under salt stress vs. control to identify salt-stress-induced genes in asparagus bean using RNA sequencing. A total of 692,086,838 high-quality clean reads, assigned to 121,138 unigenes, were obtained from control and salt-treated libraries. Then, 216 root-derived DEGs (differentially expressed genes) and 127 leaf-derived DEGs were identified under salt stress between the two cultivars. Of these DEGs, thirteen were assigned to six transcription factors (TFs), including AP2/EREBP, CCHC(Zn), C2H2, WRKY, WD40-like and LIM. GO analysis indicated four DEGs might take effects on the “oxidation reduction”, “transport” and “signal transduction” process. Moreover, expression of nine randomly-chosen DEGs was verified by quantitative real-time-PCR (qRT-PCR) analysis. Predicted function of the nine tested DEGs was mainly involved in the KEGG pathway of cation transport, response to osmotic stress, and phosphorelay signal transduction system. A salt-stress-related pathway of “SNARE interactions in vesicular transport” was concerned. As byproducts, 15, 321 microsatellite markers were found in all the unigenes, and 17 SNP linked to six salt-stress induced DEGs were revealed. These candidate genes provide novel insights for understanding the salt tolerance mechanism of asparagus bean in the future.

Introduction Soil salinity is among the leading environmental stresses affecting global agriculture, causing major reductions in crop productivity and quality [1]. This abiotic stressor is a growing problem, affecting approximately 20% of irrigated land and leading to continuing loss of arable areas in the world [2,3]. Salt stress inhibits the growth and productivity of crops via several distinct processes including photosynthesis limitation, metabolic dysfunction, and cellular structure damage [4]. Generally, there are at least two primary stresses on plants caused by salinity [5]. The first one is hyperosmotic stress resulting from the reduction of water potential and consequently diminished water availability. The second one is hyperionic stress caused by the toxic effects of accumulated ions. Salt-sensitive plants confine the absorption of salt and synthesize compatible solutes (e.g., proline, glycine betaine, and sugars) to adapt to their osmotic pressure. Salt-tolerant plants uptake and accumulate salt into the cell vacuoles, adjusting the salt concentrations of the cytosol and keeping a high cytosolic K + /Na + concentration ratio in their cells [6][7][8]. Ion exclusion mechanisms could partly improve plant salt tolerance to relatively low concentrations of NaCl, but would not work at high concentrations of salt [9].
Salt tolerance is a complex trait in various plants and is usually regulated by multiple biosynthetic and signaling pathways [10]. Actually, with increasing investigation into plant salt tolerance, three main concepts of salinity tolerance mechanisms in crop plants have been described, including osmotic tolerance, ion exclusion, and tissue tolerance [11,12]. In maintenance of osmotic homeostasis, plants commonly sequester Na + into the vacuole or accumulate compatible solutes by biosynthesis [13,14]. To regulate Na + efflux and vacuolar compartmentalization under salt stress, SOS1, a plasma membrane protein Na + /H + antiporter, was activated by interacting with the calcium sensor protein SOS3 and a Ser/Thr protein kinase SOS2 [13,15,16]. In the cytosol, a high K + /Na + concentration ratio was critical for cytoplasm ion balance during plant growth in high salinity soils [6]. Tissue tolerance is related to the ability of organs to function normally with high Na + and Clconcentration in tissues or cells [12]. A number of component factors contribute to plant tissue tolerance, including transporting of Na + and Cland the maintenance of functional water status in leaves [12,17]. Moreover, plant salt tolerance appears to be a developmentally regulated process that differs depending on plant age and developmental stages [9,18]. Recent advances have investigated salt-tolerancerelated genes, such as OsNHX2, GmsSOS1, and SOS2, that helped to better understand the key components of the plant salt tolerance network [15,16,19]. However, large gaps exist in complete comprehension of the salt tolerance trait [20].
The past decade has seen the application of next-generation DNA sequencing technologies to the analysis of crops and model species under stress [21]. In these sequence-based profiling methods, the RNA-sequencing (RNA-seq) technique has been widely employed to investigate gene expression profiles of plants at different developmental stages in response to salt stress. By comparison with earlier methods of transcriptome sequencing, the RNA-seq technique has dramatically increased the throughput of RNA sequencing and allowed comprehensive measurement of transcript abundance [22].
Asparagus bean [Vigna unguiculata ssp. sesquipedalis (L.) Verdc], a member of the Fabaceae family, is among the top five edible legumes planted worldwide [23]. Currently, soil salinity represents one of the major constraints on the productivity of asparagus bean in agriculture. Morphological feature varied in different genotypes of V. unguiculata under salt stress [24][25][26]. Although great progress has been made on the genome sequencing of cowpea (V. unguiculata) in recent years [27,28], knowledge on the genetic basis of salt tolerance is still limited in asparagus bean (V. unguiculata ssp. sesquipedalis).
Therefore, the aim of this study was 1) to investigate whole-transcriptome expression profiles of transcripts through RNA-seq in asparagus bean under salinity conditions and 2) to explore differentially expressed genes (DEGs) related to salt stress. Two asparagus bean cultivars (one salt-sensitive genotype and one salt-tolerant genotype) were treated with NaCl (125 mM) at the seedling stage. A set of non-redundant transcripts have been generated, and then they have been used for the analysis of gene annotation, functional categorization and identification of candidate NaCl stress-responsive genes.

Plant materials and growth conditions
Two asparagus bean cultivars with different salt tolerance were used in this study. One is saltsensitive type (V. unguiculata ssp. sesquipedalis var. Zhijiang14, Coded "A10") and the other is salt-tolerance type (V. unguiculata ssp. sesquipedalis var. Yulongteyou, Coded "A33") [29]. Both were collected from the Hubei Province Engineering Research Center for Legume Plants, China.
The seeds of the two cultivars were sterilized with 70% ethanol for 1 minute, followed by washing with sterile water. Filled seeds of uniform size were sown in plastic pots containing the same soil matrix. The plastic pots were placed under long-day (14 h light /10 h dark cycle) conditions at temperatures of 30˚C (light) and 20˚C (dark) in an environmental chamber (HP1000GS-B, Wuhan Ruihua, China) with 60% relative humidity. After two weeks, seedlings at the same growth status were transplanted into standard Hoagland nutrient solution for three days to acclimate at a hydroponic culture condition. On the fourth day, two-week-old seedlings were transferred to the Hoagland solution with NaCl concentration of 125 mM. Meanwhile, control seedlings were treated in the Hoagland solution without NaCl (CK). Due to different salt tolerance mechanisms stimulated by salt stress [18], six time points (0 h, 6 h, 12 h, 24 h, 48 h, and 78 h) were set for the leaf and root sampling under salt stress treatment. The "0 h" point was set for control samples (CK). Three biological replicates of roots and leaves were randomly sampled at the five points after the treatment, and then they were mixed together as one sample. Each sample was immediately frozen in liquid nitrogen and stored at -80˚C for future use. At last, roots and leaves of the two cultivars were harvested from all six time points (0 h, 6 h, 12 h, 24 h, 48 h and 78 h) for treatment and control. As a result, a total of 24 samples were prepared for RNA libraries.

RNA extraction, library construction, and sequencing
The total RNAs were extracted from roots and leaves using Colum Plant RNA extraction Kit (BioTeke, Cat. No. RP3202) according to the manufacturer's instructions. The RNA samples were digested using DNase I (RNase-free) (Fermentas, Cat. No. EN0521) at 37˚C for 30 min to remove potential genomic DNA contamination. Then, the RNAs were examined by gel electrophoresis and quantified with NanoDrop (Quawell Q5000, Quawell Technology). Integrity of the quantified RNA samples was analyzed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, U.S.A). High-quality RNAs (RIN > 8) from each time point were used for cDNA preparation and RNA-seq.
The cDNA libraries were constructed using TruSeq RNA Sample Preparation V2 (Illumina, USA) according to the manufacturer's instructions. Poly (A) mRNA was purified from total RNA using magnetic beads with oligo (dT) primers. The mRNA was then fragmented and purified with a high-temperature dilution after fragmentation. These cleaved mRNA fragments were then reverse transcribed into first strand cDNA using reverse transcriptase and random hexamers primers. Subsequently, RNA templates were removed, and the second strands were synthesized. The double-stranded cDNAs (ds cDNAs) were purified using AMPure XP beads. After purification, the ds cDNAs were subjected to end repair, single "A" nucleotides were added, and they were ligated to sequencing adapters. The products were then enriched by PCR with a cocktail primer to construct the cDNA library.
After quality control with Agilent 2100 Bioanalyzer and Qubit (Thermo Fisher Scientific) to detect fragment size and concentration, the libraries were sequenced using the Illumina HiSeq 2500 system at the Institute of Genetics and Developmental Biology, Chinese Academy of Sciences (Beijing, China).

Transcript assembly and annotation
Raw sequencing image data were transformed into raw reads by base calling. The raw reads were filtered for clean reads by Trimmomatic [30]. The criteria was set as follows: adapter mismatch � 2 bases; score of adapter palindrome mode match � 30; score of adapter simple mode match � 3; the adaptor leading and trailing � 2 bases; the sliding window 1:2; the minimized reads length � 50 bases. The subsequent analysis was based on these clean and highquality reads. The clean reads from each library were de novo assembled into contigs with Trinity software [31]. The generated contigs were blasted to the NR database to remove the contaminating bacterial sequences. The clean contigs were blasted to get rid of high-similarity sequences and to obtain sequences that can no longer be extended on either end, which were referred as all-unigene. All-unigene sequences were aligned using BLASTX against the Nr, Swiss-Port, KEGG (Kyoto Encyclopedia of Genes and Genomes), and COG (Clusters of Orthologous) databases to determine or predict protein sequence with lengths � 50 amino acids. All unigenes that could not be aligned to sequences in any of the databases mentioned above were translated into predicted protein with lengths � 100 amino acid using the software Translation. For annotations, all unigenes were searched against the NR database. To obtain the GO (Gene Ontology) terms to describe biological process, molecular function and cellular components, Blast2GO software was used for the GO annotation, and the results were functionally classified by WEGO (Web Gene Ontology) at macromolecular level. BLASTX analysis against the KEGG pathway database was performed to assign putative metabolic pathways to all-unigene [32]. Additionally, the COG database was also used for further analysis of these unigenes.

Identification of differentially expressed genes (DEGs)
DEGs under salt stress were analyzed. All the de novo clean reads were firstly assembled together to construct a reference genome, because the whole genome of V. unguiculata ssp. sesquipedalis is unavailable yet. Then each sample reads were aligned to the reference genome using the TopHat software 2.0.12 [33]. Subsequently, the TopHat-generated sequences were put into the software Cufflinks [34] and were run to calculate FPKM to indicate gene expression levels and significant differential expression between samples. Within each cultivar, DEGs were from the comparison between control (CK) and NaCl treatment during different treating time (6 h, 12 h, 24 h, 48 h, and 78 h). Between the two cultivars, DEGs were from the comparison at the same treating time and in the same tissue (root or leaf). The DEGs were defined as genes with FDR � 0.001 and the absolute ratio of log 2 (FPKM sample1/FPKM sample2) � 1 set as threshold values.
The GO terms and the KEGG pathways enriched within DEGs were also analyzed with the Fisher's exact test and the BH multiple testing correction method.
Hierarchical cluster analysis of the expression patterns of DEGs was performed after normalization by MultiExperiment Viewer v4.9 (The TM4 Software Development Team, http:// mev.tm4.org/) [35]. The GO enrichment of DEGs was analyzed by the AgriGO [36]. The GO functional enrichment analysis was tested at a significance cutoff of 0.05 false discovery rate (FDR).

Analysis of transcription factor families
Online software PlantTFcat (http://plantgrn.noble.org/PlantTFcat) [37] was employed to identify the potential genes of transcription factors (TFs) in the assembled unigene sequences derived from the salt-stress-induced RNA-seq data.

qRT-PCR validation of differential expression genes
To validate the DEGs between the two asparagus bean cultivars, candidate genes were selected randomly for qRT-PCR tests (quantitative real-time RT-PCR) run on three independent biological samples of each time point. The total RNA was isolated using the same method described above. Reverse transcription reactions were performed using PrimeScript 1st strand cDNA Synthesis Kit (TAKARA). The primers were designed using an online primer tool (http://sg.idtdna. com/calc/analyzer) to amplify 100-250 bp regions of the selected genes. The following genes and primers were used as references: gapdh (F:5'ATCAGCCAAGGACTGGAGAG3'; R: Real-time PCR was performed on an ABI 7900HT PCR instrument in a 20 μL reaction volume containing 10 μL SYBR Premix Ex Taq II (Tli RNaseH Plus), 0.5 μL primer pairs (internal standard and target genes), and 2 μL (50 ng μL −1 ) cDNA. The reactions were carried out as follows: 95˚C for 10 min followed by 35 cycles of 95˚C for 15 s, 55˚C for 60 s, and 72˚C for 20 s. With the salt-sensitive asparagus bean as a control, the relative expression levels of chosen genes were normalized to the expression levels of reference genes calculated from cycle threshold values using the 2 −ΔΔCt method [38].

Identification of microsatellites and primers calling
As byproducts, the assembled unigenes sequences were analyzed to identify potential microsatellites using software MISA (MIcroSAtellite) [39]. The criteria of microsatellites were set based on different repeat motif type. The minimum number of repeat motifs was set as following: mono-nucleotide repeats for ten, di-nucleotide repeats for six, tri-nucleotide repeats for five, tetra-nucleotide repeats for five, penta-nucleotide repeats for five and hexa-nucleotide repeats for five. Maximal number of bases interrupting two SSRs in a compound microsatellite was less than 100 base pairs. Microsatellite primers were designed using Primer 3 (http:// primer3.sourceforge.net/).

Mapping of the salt-stress-induced DEGs in asparagus bean genome
To identify the genomic location of the salt-stress-induced DEGs in asparagus bean, alignment between the DNA sequences of the DEGs and newly-published genome sequences of Vigna unguiculata v1.0 (GenBank accession: MATU00000000.1, https://legumeinfo.org/blat). Subsequently, SNP loci linked to the DEGs regions were revealed by alignment between the genome sequence Vigna unguiculata v1.0 and a recently published SNP map of V. unguiculata [40]. After that, the genomic position of the DEGs and their adjacent SNP loci could be identified.

Illumina sequencing, quality filtering, and de novo assembly
In this study, we constructed 24 cDNA libraries, including leaf and root samples from the two asparagus bean cultivars and six time points (0 h, 6 h, 12 h, 24 h, 48 h, and 78 h) for treatment and control. These libraries were then sequenced on the Illumina HiSeq 2500 platform. After removing sequencing adaptors and low quality data, 692,086,838 filtered reads were obtained amounting to 86.48 Gb of data (S1 Table). Transcriptome was de novo assembled using the Trinity software. Statistics of all the unigenes was listed in Table 1. 121,138 unigenes were generated, with a total length of 150,960,298 bp. The lengths of all-unigene ranged from 201 to 32,391 bp. The average unigene length was 968 bp with the N50 of 1,720 bp.

Functional annotation and classification of assembled transcripts
After assembly, all the 121,138 unigene sequences were searched against several public databases, including the NR, GO, COG, and KEGG databases using BLAST. 65,322 unigene sequences with each length less than 500 bp had BLASTX hits in the NR database (S2 Table). Among them, 28,019 unigene sequences were mapped to known genes with the best hit (E-value < 1e -50 ) (Fig 1A). Based on similarity distribution analysis, 14,421 unigene sequences matched deposited sequences with a similarity > 80% ( Fig 1B). Moreover, 67% of the annotated unigenes could be assigned to the sequences from the top two legume species: Glycine max (59%) and Medicago truncatula (8%) (Fig 1C). 27,936 unigenes were annotated based on the GO analysis. Notably, biological regulation, cellular processes, metabolic processes, and response to stimulus were significantly over-represented in the 23 biological process GO groups. In the further analysis of GO classification, 27,794 unigenes was assigned to 25 COG categories (S3 Table). "General function prediction only" represented the largest group (5,681), followed by "post-translational modification, protein turnover, chaperones" (2,872).

DEGs of asparagus bean in response to salt stress between the two cultivars
In this study, a set of salt-stress-induced DEGs were obtained (Fig 2; S6 Table). In the roots of the two cultivars, 467 DEGs were revealed including 235 upregulated DEGs and 232 were found in their leaves (Fig 2A and 2C; S6 Table). After removing the repeated unigenes in all the DEGs, a total of 127 leaf-based DEGs and 216 root-based DEGs were identified between "A10" and "A33" (Fig 2B and 2D). Compared with the salt-tolerance cultivar "A33", the salt-sensitive cultivar "A10" triggered much more DEGs both in roots and in leaves. As shown in Fig 2A, regarding the number of up-regulated DEGs in roots, there was forty-eight DEGs in the "A33" and seventy-nine DEGs in the "A10". Similarly, the number of down-regulated DEGs was thirty in the "A33" roots and seventy-eight in the "A10" roots.
Furthermore, gene expression patterns of roots and leaves were analyzed between "A10" and "A33" under salt stress conditions. As shown in Fig 3, a heatmap of the 216 root-based DEGs exhibited temporal and spatial expression patterns between the two cultivars during salt stress treatment. These novelty DEGs would help to reveal the differentiation of salt tolerance occurred between the two cultivars under salt stress. Salt-related DEGs in asparagus bean
Among the thirteen TFs, three TFs ((CCHC(Zn) (c53466_g1_i2), CCHC(Zn) (c53535_g7_i1) and C2H2 (c58717_g1_i2) illustrated opposite regulation under salt stress between the two cultivars. As can be seen from the Fig 5B, c53535_g7_i1 and c58717_g1_i2 showed up-regulation in both roots and leaves of "A10". On contrary, both of them were down-regulated in both roots and leaves of "A33". Similarly, c53466_g1_i2 was down-regulated in "A10", but it was up-regulated in "A33". Thus, these TF-related DEGs would play diverse roles in salt stress response and adaptation in asparagus bean.

Validation of salt-stress-induced DEGs by qRT-PCR between the two cultivars
To validate the RNA sequencing results and expression profiling, nine salt-stress-induced DEGs were randomly selected for qRT-PCR analysis ( Table 2). The root tissues of the two cultivars ("A10" and "A33") were sampled for 0, 12, 24, and 48 h under salt stress and control. The qRT-PCR results demonstrated that the nine DEGs were salinity responsive. And most of them exhibited different responses to salt stress in both cultivars "A10" and "A33" (Fig 6).

Fig 4. The GO enrichment of the four DEGs between "A10" and "A33" at CK (0h) and at the salt stress treatment stages (6 h, 12 h, and 24 h). (A)
The GO enrichment of the DEGs between the "A10" root and the "A33" root at CK (0 h) and at the salt stress treatment stages (6 h, 12 h, and 24 h). (B) The GO enrichment of the DEGs between the "A10" leaf and the "A33" leaf at CK (0 h) and at the salt stress treatment stages (6 h, 12 h, and 24 h). (C) and (D) The heatmaps of the DEGs of enriched GO terms between the "A10" root and the "A33" root at CK (0 h) and at the salt stress treatment stages (6 h, 12 h, and 24 h). (E) The heatmaps of DEGs of enriched GO terms between the "A10" leaf and the "A33" leaf at CK (0 h) and at the salt stress treatment stages (6 h, 12 h, and 24 h). https://doi.org/10.1371/journal.pone.0219799.g004

The vesicular-transport-related KEGG pathway under salt stress in asparagus bean
Further analysis of these DEGs was conducted to investigate the biological functions in the 125 KEGG pathways (S7 Table). In the present study, notably, the pathway "SNARE interactions in vesicular transport" (ko:vvi04130), which was a salt-related process, was enriched. Six genes (c92920_g1_i1, c64313_g1_i1, c45795_g1_i2, c49659_g1_i1, c98985_g1_i1, c33236_g2_i1) were involved in "SNARE interactions in vesicular transport" (Fig 7). Thus, the "SNARE interactions in vesicular transport" pathway might be affected under salt stress conditions, implying a potential contribution to the regulation of salt stress in asparagus bean.
The position of these DEGs and their SNPs was newly found, and they were not reported before. So, they provide new sight into asparagus bean genome concerning salt-stress-induced genes.

Functional annotation, classification of assembled unigenes
In the present study, our results increased the number of transcript sequences available that could be integrated to public genome resources of asparagus bean [27,28,[41][42][43].We Salt-related DEGs in asparagus bean identified 121,138 unigenes with an average length of 968 bp and N50 length of 1,720 bp. Based on GO analysis, there was 23 biological process groups. The dominant GO terms of the unigenes under salt stress concurred with previous similar studies, comprising of "cellular process" and "metabolic process" [44,45]. A set of 343 DEGs (127 leaf-based DEGs and 216 root-  based DEGs) were identified in asparagus bean seedlings between control and salt-treated conditions.
The salt-stress-induced DEGs in the processes of "oxidation reduction", "transport" and "signal transduction" In higher plants, intracellular and extracellular antioxidants could give rise to complicated networks to protect against oxidation activities under abiotic stresses [46]. In this study, four novel DEGs were uncovered by the GO analysis. Three of them (c58509_g1_i1, c88702_g1_i1, and c57788_g1_i2) were in the root, while one of them (c58967_g9_i1) was in the leaf. They were predicted function in "oxidation reduction" (c58509_g1_i1, c88702_g1_i1) ( Fig 4C) "transport" (c57788_g1_i2) (Fig 4D) and "signal transduction" (c58967_g9_i1) (Fig 4E). In the salt-tolerance genotype "A33", the two genes (c57788_g1_i2, c88702_g1_i1) were downregulated.
Generally, root is the plant organ that first senses and reacts to environmental stresses [47], and salt stress suppressed root growth in most plants [48]. Notably, in the oxidation reduction, gibberellin 2-oxidase (c58509_g1_i1) gene of the gibberellin catabolic process was up-regulated in the roots of "A33" under salt stress in this study. Gibberellin 2-oxidase gene encodes 2-oxoglutarate-dependent dioxygenase (2ODD), which changes bioactive and intermediate forms of GA to inactive forms [49]. Endogenous GA levels in Arabidopsis was decreased by the induction of GA 2-oxidase, causing that growth is repressed for high-salinity stress [50]. Recent results suggest that GA 2-oxidase is crucial for controlling Arabidopsis root meristem cell number and suppresses IAA-directed primary root and root hair growth in response to salt stress [51,52]. Similarly, the candidate gene of gibberellin 2-oxidase (c58509_g1_i1) possibly stimulate the oxidation reduction of asparagus bean in response to salt stress.
One DEG "c57788_g1_i2" was a candidate gene of transmembrane transporters NRT2 (nitrate transporter 2). The NRT2 family, which contains seven genes in Arabidopsis, encodes high-affinity NO 3transporters in roots [53]. Except for nitrate availability, the NRT2 genes could be regulated at the transcriptional level by other stress factors [54][55][56]. In our study, the DEG "c57788_g1_i2" was down-regulated in the salt-tolerance genotype "A33". As a result, its low expression level possibly affected root NO 3− uptake in asparagus bean, and finally caused growth reduction under low and limiting NO 3− availability.
In addition, a potential candidate gene (c58967_g9_i1), encoding protein suppressor of NPR1 (nonexpressor of pathogenesis-related protein 1), was related to signal transduction. Actually, the function of NPR1 is a crucial modulator of the plant UPR (unfolded protein response) that accompanies with ER stress-induced reduction of the cytosol and translocation of NPR1 to the nucleus [57]. Under salt stress, NPR1 might regulate SA signaling to mediate Salt-related DEGs in asparagus bean Na + entry into roots and its subsequent transport into shoots [58]. So, this might be the reason that the DEG "c58967_g9_i1" showed up-regulation in the salt-sensitive genotype "A10" under salt stress. Undoubtedly, signal transduction system is pivotal for asparagus bean in defense responses against salinity.

Novel TF-related DEGs under salt stress revealed in asparagus bean
Transcription factors regulated expression of salt-related genes and ultimately determined the level of salt tolerance of plants [3]. Many of these transcription factors relating to abiotic stress were predominant, including "WD40-like" [59], "CCHC(Zn)" [60], "AP2/EREBP" [61], "basic helix-loop-helix (bHLH)" [62], "WRKY" [63], and "bZIP" [64,65]. Similar transcription factors were found and they exhibited different expression regulation in asparagus bean in this study. Four TFs (AP2/EREBP, CCHC(Zn), C2H2 and WRKY) in the leaf-based DEGs and five TFs (AP2/EREBP, C2H2, CCHC(Zn), WD40-like and LIM) were in the root-based DEGs. It is noticed that gene expression were up-regulated in four TFs (AP2/EREBP, CCHC(Zn), C2H2, and LIM) in the roots of salt-sensitive genotype "A10". They may play an important role in coping with salt stress in asparagus bean. As for AP2/EREBP, the increased level of expression of AP2/EREBP improves salt tolerance in Oryza sativa and Arabidopsis [66]. C2H2-type zinc fingers contribute to the potato response to salt and dehydration stresses through an ABA-dependent pathway [67]. C2H2 zinc-finger genes probably responded positively to salt stress and were upregulated in leaves and/or roots in poplar [68,69]. In Arabidopsis C2H2-type zinc-finger proteins functioned as transcription repressors under salt stress conditions [70,71]. During the salt stress process in Arabidopsis, overexpression of the CSDP1 gene (containing seven CCHC-type zinc fingers) suppressed seed germination, while overexpression of the CSDP2 (containing two CCHC-type zinc fingers) gene promoted seed germination [72]. In tomato, a novel regulatory gene, SlZF3, encodes a C2H2 zinc-finger protein and enhances plant salt-stress tolerance by interacting with CSN5B [73]. In wheat, two proteins (TaRZ1 and TaRZ2) containing a CCHC-type zinc finger were inhibited under salt stress conditions [74]. These novel TFs provided potential candidate genes for future understanding of the molecular mechanism of salt tolerance in asparagus bean.

DEGs involved in ion transport, hormones, and signal transduction
The processes of ion transport, hormones, and signal transduction have important roles in response to salinity stress in asparagus bean. The nine tested DEGs are mainly involved in those processes. Several DEGs had potential for maintaining homeostasis in the ion transport process. Their function was associated with ion transmembrane transport activities, cation transport (c58342_g1_i1), substrate-specific transmembrane transporter activity (c33878_g1_i2), response to osmotic stress (c17757_g1_i1), and cytoplasmic localization (c59189_g3_i2, c59535_g2_i1). Another DEG (c58509_g1_i1) was related to the gibberellin catabolic process in asparagus bean. Previous studies showed that GA enhanced by salt stress, resulting in the irregular root growth [75,76]. Therefore, to cope with salt stress, it is essential to activate plant hormone metabolism pathways in asparagus bean. Additionally, one of the DEGs (c87471_g1_i1), down-regulated in "A33", was probably involved in the phosphorelay signal transduction system in asparagus bean under salt stress. Salt stress could stimulate the phosphorylation of the Arabidopsis vacuolar K + Channel TPK1 by calcium-dependent protein kinases [77]. These novel DEGs might be the reason for the discrepancy of salt tolerance in asparagus bean.

Mechanism of the vesicular-transport-related KEGG pathway under salt stress
Intracellular membrane dynamics play a significant role in plant salt tolerance [78]. In the KEGG pathway analysis, it is noticed that several DEGs were assigned to the pathway "SNARE interactions in vesicular transport" which mediates vesicular trafficking. Many researchers have studied SNARE complexes in plant tolerance to salt stresses [79][80][81][82][83][84][85]. The AtVAMP7C family of vesicle soluble N-ethylmaleimide-sensitive factor attachment protein receptors (v-SNAREs) mediates fusion of H 2 O 2 -containing vesicles with the tonoplast. Suppression of the AtVAMP7C genes that showed downregulation effect could improve plant salt tolerance [86]. Similarly, overexpression of GsVAMP72, which is a member of R-SNARE family in Glycine soja, resulted in the reduction of plant tolerance to salt stress [83]. Additionally, a Qc-SNARE protein (AtSFT12) took part in salt response via sequestration of Na + in vacuoles [84]. In a word, as an important strategy of osmotic tolerance, this "SNARE interactions in vesicular transport" pathway influenced ion influx or efflux. It should play an important role in cellular metabolic processes, ultimately leading to the salt tolerance disparity among various asparagus bean cultivars.

Molecular markers and mapping of the salt-stress-induced DEGs in asparagus bean
Molecular markers have been playing an increasing role in plant genetics studies. As byproducts in this study, a great number of microsatellites were screened out from sequences of all the unigene. Further, 15, 321 novelty microsatellite primers were obtained. Once the marker validation of these microsatellites is verified, they could be useful tools for future studies on genome mapping, marker-assisted selection, and population genetic analysis in asparagus bean. In addition, mapping of six salt-stress-induced DEGs exhibited their distribution among six chromosomes in asparagus bean genome. 17 SNP loci are linked to the six salt-stressinduced DEGs, which are different from the seven SNP markers related to salt tolerance at seedling stage in V. unguiculata [87]. These salt-stress-related markers (both microsatellites and SNP markers) would offer paths to shed light on the molecular genetic basis of salt tolerance in asparagus bean.

Conclusion
In this study, a salt-stress-induced transcriptome of asparagus bean was obtained and annotated. A set of 343 DEGs (127 leaf-based DEGs and 216 root-based DEGs) were identified at seedlings of asparagus bean between control and salt-treated conditions. Four DEGs might be involved in the biological processes of "oxidation reduction", "transport" and "signal transduction". Seven leaf-based DEGs were assigned to four TF families (AP2/EREBP, CCHC(Zn), C2H2, and WRKY), while six root-based DEGs were found in five TF families (AP2/EREBP, C2H2, CCHC(Zn), WD40-like and LIM). Expression profiles of nine salt-stress-induced DEGs which involved in ion transport, hormones, and signal transduction were validated by qRT-PCR. In the KEGG pathway analysis, six genes were enriched in the salt-stress-related pathway "SNARE interactions in vesicular transport". In addition, 15, 321 microsatellite markers were found in all the unigenes, while 17 SNP linked to six salt-stress induced DEGs were revealed. These candidate genes provide novel insights for understanding the salt tolerance mechanism of asparagus bean in the future.
Supporting information S1