De novo transcriptome sequencing and analysis of salt-, alkali-, and drought-responsive genes in Sophora alopecuroides.

BACKGROUND: Salinity, alkalinity, and drought stress are the main abiotic stress factors affecting plant growth and development. Sophora alopecuroides L., a perennial leguminous herb in the genus Sophora, is a highly salt-tolerant sand-fixing pioneer species distributed mostly in Western Asia and northwestern China. Few studies have assessed responses to abiotic stress in S. alopecuroides. The transcriptome of the genes that confer stress-tolerance in this species has not previously been sequenced. Our objective was to sequence and analyze this transcriptome.
RESULTS: Twelve cDNA libraries were constructed in triplicate from mRNA obtained from Sophora alopecuroides for the control and salt, alkali, and drought treatments. Using de novo assembly, 902,812 assembled unigenes were generated, with an average length of 294 bp. Based on similarity searches, 545,615 (60.43%) had at least one significant match in the Nr, Nt, Pfam, KOG/COG, Swiss-Prot, and GO databases. In addition, 1673 differentially expressed genes (DEGs) were obtained from the salt treatment, 8142 from the alkali treatment, and 17,479 from the drought treatment. A total of 11,936 transcription factor genes from 82 transcription factor families were functionally annotated under salt, alkali, and drought stress, these include MYB, bZIP, NAC and WRKY family members. DEGs were involved in the hormone signal transduction pathway, biosynthesis of secondary metabolites and antioxidant enzymes; this suggests that these pathways or processes may be involved in tolerance towards salt, alkali, and drought stress in S. alopecuroides.
CONCLUSION: Our study first reported transcriptome reference sequence data in Sophora alopecuroides, a non-model plant without a reference genome. We determined digital expression profile and discovered a broad survey of unigenes associated with salt, alkali, and drought stress which provide genomic resources available for Sophora alopecuroides.

Background

At present, it is an indisputable fact that global climate has changed, posing a potential threat to the sustainable development of agriculture and food security [1]. Increasing global temperatures cause sea level to rise, which in turn increases the salinity of groundwater in coastal and arid areas [1]. Salt-alkali land is widely distributed around the world, covering about 100 million hectares [2]. According to the Food and Agriculture Organization of the United Nations, more than 400 million hectares of land on the major continents are affected by salt [3]. Rising salinization may reduce agricultural acreage by up to 20% per year by 2050 [4]. The amount of land in India that has been degraded by excess sodicity and salinity is estimated to be about 6.75 million hectares [5]. In China, there are about 100 million hectares of salinized land [6]. Drought, which can cause salinity to increase, has a great impact on crop yield [7, 8]. With the changes in global climate, the frequency and duration of drought events is increasing, with serious impacts on crop yields [9, 10]. To solve the growing global food shortage, it is essential to use saline-alkali land for agriculture. Using effective gene resources to cultivate salt-, alkali-, and drought-resistant crops is the most economical and productive measure to solve this problem.

Plant breeding and gene transformation are important ways to improve crop tolerance to abiotic stress. Stress-regulatory mechanisms in higher plants have been analyzed by researching many of the genes related to abiotic-stress tolerance at the transcriptional level [11, 12]. Various signal transduction pathways are involved in plant responses to abiotic stress; these include the phospholipid signaling pathway [13], calcium-dependent protein kinase pathway, mitogen-activated protein kinase cascade pathway [14], and abscisic acid (ABA) pathway [15]. These pathways form a signal transduction network by which plants respond to abiotic stress [16]. In addition, stress tolerance mechanisms include a series of gene expression and gene product interactions, which enhance plant adaptations to abiotic stress at cellular and molecular levels [17]. Many differentially expressed genes (DEGs), which encode reactive oxygen species scavenging proteins, aquaporins, heat shock proteins, and ion transporters have been identified in stress resistance [18].

Technological developments have provided great convenience in biological science research. Next-generation sequencing of RNA, which can directly determine cDNA sequences, has been widely used to identify plant genes [19, 20]. RNA-seq and DEG-analysis have revealed mechanisms of responses to complex biotic and abiotic stressors in many plant species [21], including Arabidopsis thaliana [22], Vitis vinifera [23], Ammopiptanthus mongolicus [24], Cucumis sativus [25], and Gossypium hirsutum [26]. For instance, in A. thaliana, about 30% of the transcriptome is considered to be involved in abiotic-stress regulation, and 2409 genes have been identified as being of great importance in drought resistance, salt tolerance, and resistance to cold [27, 28].

Sophora alopecuroides L. (Fabaceae) is a highly stress-tolerant leguminous perennial herb in the genus Sophora, distributed mainly in western Asia and northwestern China [29, 30]. Sophora alopecuroides is an important potential resource for stress resistance genes. However, few studies have focused on finding stress resistance genes in S. alopecuroides using the transcriptome sequencing method. In this study, we perform transcriptome sequencing of S. alopecuroides plants subjected to three stress treatments (salt, alkali, and drought). We use de novo sequence assembly and differential gene expression analysis, and screen many genes related to abiotic resistance. Our study provides new genetic resources for research of abiotic resistance in crop plants, thereby increasing the options for genetic crop improvements.

Results

Sophora alopecuroides is a type of perennial herb and drought tolerant plant, with its drought tolerance closely related to its root system. Using the method of saline, alkaline and drought treatments used on Arabidopsis and soybean, it was found that when concentrations of NaCl, NaHCO3 and PEG were more than 1.2, 1.2 and 8% respectively, the growth of Sophora alopecuroides was inhibited or the plant wilted (Fig. S1). After 72 h, plant growth and physiological indexes were significantly different and relatively stable. Therefore, we used 1.2% NaCl, 1.2% NaHCO3 and 8% PEG-treated Sophora alopecuroides roots as tissues for constructing a cDNA library for transcriptome sequencing, from which differences in gene expression under saltine, alkaline and drought conditions could be explored.

Transcriptome sequencing and assembly

In total, 605,800,814 raw sequencing reads were obtained from the control and treated samples. And 586,189,628 clean reads were used to gather the data (Table 1). 1,382,370 transcripts and 902,812 unigenes were obtained, with an average length of 366 bp and 294 bp, respectively. Six hundred eighty-six thousand one hundred twenty-nine unigenes were 200 to 500 bp long, and 80,452 were > 1000 bp long (Fig. 1).

Table 1

Summary of Sophora alopecuroides sequences analyzed

Sample	Raw Reads	Clean reads	Clean bases	Error (%)	Q20 (%)	Q30 (%)	GC (%)
CK_1	47,691,636	46,202,320	6.93G	0.01	97.79	94.15	44.01
CK_2	41,826,612	40,474,310	6.07G	0.01	97.90	94.42	43.73
CK_3	50,190,638	49,125,200	7.37G	0.01	97.61	93.93	43.91
ST_1	46,150,030	44,279,092	6.64G	0.01	97.32	93.51	43.22
ST_2	50,369,242	48,508,750	7.28G	0.01	97.82	94.37	43.51
ST_3	54,517,670	52,656,970	7.90G	0.01	97.85	94.40	43.48
A_ST_1	48,586,426	47,023,258	7.05G	0.01	97.67	93.94	46.67
A_ST_2	58,315,134	56,447,032	8.47G	0.01	97.63	93.90	46.23
A_ST_3	53,214,668	51,438,942	7.72G	0.01	97.54	93.68	44.09
DT_1	57,566,808	55,738,546	8.36G	0.01	97.76	94.20	43.46
DT_2	49,791,136	48,134,912	7.22G	0.01	97.91	94.42	43.59
DT_3	47,580,814	46,160,236	6.92G	0.01	97.92	94.42	43.62
Summary	605,800,814	586,189,568	87.93G

The numbers 1–3 after CK, ST and A_ST, and DT identify the three independent biological replicates for the control and salt, alkali, and drought treatments, respectively

Q20: The percentage of bases with a Phred value > 20

Q30: The percentage of bases with a Phred value > 30

Fig. 1

De novo assembly length distribution of sequences for Sophora alopecuroides. Transcripts: red; Unigenes: blue

Functional annotation of all non-redundant unigenes

The annotated number unigenes for each database is shown in Table 2. Among the unigenes, 357,522 (39.6%) had significant matches in the Nr database (NCBI redundant protein Sequences), 214,419 (23.75%) in the Nt database (NCBI nucleotide Sequences), and 293,553 (32.51%) in the Swiss-Prot database. Among the 902,812 unigenes, 545,615 (60.43%) had at least one highly match with an identified gene in BLAST searches (Table 2).

Table 2

BLAST analysis of non-redundant unigenes sequenced for Sophora alopecuroides, against public databases

	Number of Unigenes	Percentage (%)
Annotated in Nr	357,522	39.6
Annotated in Nt	214,419	23.75
Annotated in KO	157,394	17.43
Annotated in Swiss-Prot	293,553	32.51
Annotated in PFAM	356,271	39.46
Annotated in GO	366,814	40.63
Annotated in KOG	178,137	19.73
Annotated in all Databases	47,161	5.22
Annotated in at least one Database	545,615	60.43
Total Unigenes	902,812	100

Nr NCBI non-redundant protein Sequences, Nt NCBI nucleotide Sequences, Pfam Protein family, KOG/COG KOG: euKaryotic Ortholog Groups; COG: Clusters of Orthologous Groups of proteins, Swiss-Prot A manually annotated and reviewed protein sequence database, KEGG Kyoto Encyclopedia of Genes and Genomes, and GO Gene Ontology

Functional classification by GO and KOG

The GO (Gene Ontology) classification that we used includes three main classes of ontology. The salt-, alkali-, and drought-treatment samples were examined by GO functional significant enrichment analysis. For the salt-treatment samples, 1178 DEGs were annotated into 47 categories; 5863 DEGs were annotated into 59 categories for the alkali-treatment samples; and 2232 DEGs annotated into 60 categories for the drought-treatment samples. In the biological process, the most enriched categories in salt- and drought-treatment were the biosynthetic, organic substance biosynthetic, and cellular biosynthetic processes. In contrast, in the alkali-treatment samples, the most enrichment occurred in the metabolic process category. In the cellular component, the most enrichment occurred in relation to cellular morphology, cell, and intracellular. In the molecular function, the most enrichment occurred in relation to structural constituents of ribosome, structural molecule activity, and molecular function (Fig. 2).

Fig. 2

Histogram of GO classification for Sophora alopecuroides. The results are summarized in three main categories: Biological Process, Cellular Component, and Molecular Function. The x-axis indicates the subcategories, and the y-axis shows the number of genes associated with the GO terms. The subset “ST&CK” (panel a) indicates the number of DEGs between the salt treatment and control, “A_ST&CK” (panel b) the number of DEGs between the alkali treatment and control, and “DT&CK” (panel c) the number of DEGs between the drought treatment and control

One hundred seventy-eight thousand one hundred thirty-seven unigenes were annotated to 26 groups in KOG database. Among these groups, the largest were those involved in protein turnover, post-translational modification, and chaperones (28271), followed by translation, ribosomal structure, and biogenesis (27448), general function prediction (20248), signal transduction mechanisms (15600). Few unigenes relate to groups involved in cell motility (223) and extracellular structures (245) (Fig. 3).

Fig. 3

Functional classification of the assembled unigenes for Sophora alopecuroides. The y-axis indicates the percentage of genes annotated relative to all the annotated genes

Functional classifications using KEGG pathways under salt stress

All unigenes were annotated and mapped to KEGG database (http://www.genome.jp/kegg/). Of the 1673 DEGs sequenced from the salt-treatment samples, 616 were annotated to 61 metabolic pathways (Table S1). Of these, 64 up-regulated DEGs were annotated in 25 of the metabolic pathways, and 552 down-regulated genes were annotated in 55 of the metabolic pathways. In these 61 metabolic pathways, only 6 were annotated to the up-regulation of DEGs, and 19 metabolic pathways annotated both up-regulated and down-regulated genes. Only 36 metabolic pathways annotated down-regulated genes.

In the process of plant secondary metabolite synthesis, the coenzyme A gene, SaCoA, involved in the phenylpropanoid biosynthesis pathway (ko00940) and phenylalanine metabolism (ko00360), corresponds to the upregulation of differential genes under salt stress. Coenzyme A is an important cofactor in many biosynthesis, degradation, and energy generation pathways [31]. Previous studies have shown that the coenzyme A biosynthesis enzyme phosphoryltransferase participate in plant growth, salt resistance, and osmotic stress [32]. The coenzyme A biosynthetic pathway corresponds to regulating plant salt tolerance [33]. In the study of Zygophyllum spp., it was found that under salt stress conditions, the CoA contents of the salt-tolerant varieties in the control group and the salt-treated group did not differ significantly, whereas the CoA contents of the salt-sensitive varieties decreased significantly [33]. In the ABA signaling pathway, SaPYL4–1, SaPYL4–2, SaPYL4–3, SaPYL4–4, and SaPYL5–1 were found to be related to ABA receptors (Table S2). These five genes were down-regulated both under salt and alkali treatments. Four upregulated DEGs, SaPP2C8, SaPP2C16, SaPP2C37, and SaPP2C53, were related to protein phosphatase 2C (Table S2). Moreover, SaPP2C8 and SaPP2C53 also showed up-regulation under alkali stress. This result is consistent with the confirmed relationship between PYL and PP2C in the ABA signaling pathway [34]. In the signaling pathway of the plant hormone brassinolide, a gene SaTCH4 related to xyloglucosyl transferase TCH4 was identified under salt stress. This gene is involved in the regulation of cell elongation, which may be related to the suppression of plant growth under salt stress conditions [35].

Functional classifications using KEGG pathways under alkali stress

Of the 8142 DEGs sequenced from the alkali-treatment samples, 2644 DEGs were annotated to 118 metabolic pathways in the KEGG database. Of these, 2056 up-regulated DEGs were annotated in 117 metabolic pathways, and 588 down-regulated DEGs were annotated in 64 metabolic pathways. Among these 118 metabolic pathways, there were only 54 metabolic pathways that contained up-regulated DEGs. There were 63 metabolic pathways that annotated both up-regulated and down-regulated DEGs, whereas only one metabolic pathway annotated down-regulated genes.

Under alkali stress, the positive regulatory gene SaNPR1 was obtained in the signal transduction pathway. The SaNPR1 gene was annotated as a regulatory protein NPR1-like, which is a positive regulatory gene in the salicylic acid signal pathway. Studies have shown that NRP1 participates in abiotic stresses including low temperature and salt stress through the salicylic acid signaling pathway [36-40], and that salicylic acid, as a plant stress signal, plays an important role under high pH stress conditions [41]. In the ethylene signaling pathway, 17 genes related to serine/threonine-protein kinase CTR1 were upregulated. These included SashkB, SashkC1, SashkC2, SashkC3, SashkC4, SakinX, SadrkD, SagefX, SaDDB1, SaDDB2, SaDDB3, SaDDB4, SaDDB5, SaMIMI1, SaMIMI2, Sapats1–1, and Sapats1–2 (Table S2). Studies have confirmed that CTR1 was a positive factor regulating abiotic stress [42]. Six negative regulatory genes were obtained in the auxin signaling pathway, including the auxin influx carrier (AUX1 LAX family) gene SaAUX1, and auxin-responsive protein IAA genes SaIAA8–1, SaIAA8–2, SaIAA14, SaIAA26, and SaIAA27 (Table S2). These genes are mainly involved in cell enlargement and plant growth. The histidine kinase 4 (cytokinin receptor) gene SaAHK4, which is involved in cytokine signaling pathway, also exhibits negative regulation under drought stress. In the gibberellin signaling pathway, the negatively regulated DELLA protein GAI-like gene SaGAI was obtained, which was mainly involved in plant stem growth and induction of germination. In the brassinolide signaling pathway, the negatively regulated D3-type cyclin isoform 1 gene, SaCYCD3, was obtained, which was mainly involved in the cell division process. The negative regulatory genes we obtained are mainly involved in the growth and development of plants [43-52]. Previous studies in Arabidopsis have found that plants can further improve their resistance to stress by slowing down growth and promoting leaf senescence [34]. Furthermore, we obtained two up-regulated genes in the plant secondary metabolite synthesis pathway, namely SaGCDH and SaOMT6. The caffeoyl-CoA O-methyltransferase gene is related to the Citrus reticulata flavonoid biosynthesis process, and flavonoids, as the primary secondary metabolites, play a crucial part in plant stress resistance [42]. The up-regulated differential genes related to antioxidant enzymes mainly include SaPXR1, superoxide dismutase genes (SaSOD1, SaSOD2–1, and SaSOD2–2), and the putative peroxisomal-coenzyme A synthetase gene SaHACL1. These genes are mainly involved in the removal of ROS responding to stress, and thus reduce the damage of plant cells and tissues by reactive oxygen species [53, 54].

Functional classifications using KEGG pathways under drought stress

Of the 17,479 DEGs sequenced from the drought-treatment samples, 6546 DEGs were assigned to 121 metabolic pathways; of these, 576 up-regulated DEG annotations were in 101 metabolic pathways, and 5970 down-regulated DEG annotations were in 120 metabolic pathways. Of the 121 metabolic pathways, one was annotated to up-regulated DEGs, and 100 contained both up-regulated and down-regulated genes; the other 20 annotated only down-regulated genes.

Under drought stress, multiple positively regulated expression genes were obtained in the phytohormone signal transduction pathway and found to participate in the ABA signaling pathway. Four genes, namely SaPYL4–1, SaPYL4–2, SaPYL5, and SaPYL9, were found to be related to the abscisic acid receptor PYR/PYL family (Table S2). In a study on Arabidopsis abiotic stress, it was found that ABA receptor-related genes are accompanied by up-regulation of ABA to activate the Arabidopsis stress resistance system, which is a positive regulator of Arabidopsis adaptation to abiotic stress [34]. In this process, PYL-related genes promote the expression of serine/threonine-protein kinase gene by inhibiting the expression of PP2C-related genes [34]. Further analysis revealed four serine/threonine-protein kinase genes (SaSRK2e-1, SaSRK2e-2, SaSRK2e-3, and SaSAPK1) and one ABA response element-binding protein 1 gene (SaABF1). These genes are consistent with the expression pattern of Arabidopsis ABA in response to stress [42, 55]. Therefore, we presumed that the up-regulated expression of these genes may responsible for the drought tolerance of Sophora alopecuroides. In addition, we obtained four negative regulatory genes involved in plant hormone signal transduction pathways, including SAUR-like auxin-responsive family protein (SaSAUR32), sensory histidine protein kinase (SaAHK2 and SaAHK4), and histidine-containing phosphotransfer protein 1 (SaAHP1), which are mainly involved in cell growth, cell division, bud germination, and plant growth. Studies have confirmed that under drought stress, Arabidopsis can respond to adversity stress by weakening its growth [42, 55, 56]. We speculate that there is a similar mechanism in Sophora alopecuroides, and the downregulated expression of genes SaSAUR32, SaAHK2, SaAHK4, and SaAHP1 is a stress response. In secondary metabolite synthesis-related pathways, we identified upregulated genes, including shikimate O-hydroxycinnamoyl transferase (SaHST), catalase isozyme 1 (SaCAT1–1, SaCAT1–2, SaCAT1–3, SaCAT1–4, and SaCAT1–5), and peroxisome biogenesis protein 5 (SaPEX5) (Table S2). These genes are mainly involved in the active oxygen scavenging mechanism. After the plant is subjected to adversity stress, it will be accompanied by secondary stress damage, such as that by reactive oxygen species. In response to oxidative stress, plants form peroxidase, superoxide dismutase, and catalase, which are used to remove active oxygen species and reduce the damage caused by them to plant cells [53, 54].

Analysis of differential gene expression

The number of DEGs under saline, alkaline and drought conditions was quite different (Fig. 4). The number of DEGs in the drought treatment was 17,479, of which 2036 were upregulated and 15,443 were downregulated. In the alkaline treatment 8142 DEGs were obtained, of which 6271 were upregulated and 1871 were downregulated. There were only 1673 DEGs under saline stress conditions, with 159 upregulated and 1514 downregulated. Of all the DEGs, four were upregulated in all 3 treatments and 899 were downregulated.

Fig. 4

Venn diagram of DEGs sequenced for Sophora alopecuroides. The sum of the numbers in each large circle represents the total number of DEGs between combinations. The overlapping part of the circles represents DEGs for the treatment combinations. “ST&CK”: number of DEGs between salt treatment and control; “A_ST&CK”: number of DEGs between alkali treatment and control; “DT&CK”: number of DEG between drought treatment and control. a: up-regulation of DEGs; b: down-regulation of DEGs

Furthermore, 11,936 DEGs were annotated into 82 transcription factor families, including 814 MYB transcription factors, 472 bZIP transcription factors, 166 WRKY transcription factors, and 123 NAC transcription factors (Table S3). There were 42 MYB transcription factors corresponding to 6 KEGG metabolic pathways, and 25 bZIP transcription factors corresponding to 7 KEGG metabolic pathways. However, one NAC transcription factor corresponded to only one KEGG metabolic pathway, and 25 WRKY transcription factors corresponded to 6 KEGG metabolic pathways.

Among the differential genes, three MYB transcription factors, SaMYB1, SaMYB5, and SaMYB14, were identified up-regulating under salt and alkali stress (Table S2). SaMYB1 belongs to LHY (K12133) of the circadian rhythm – plant pathway (Ko04712) in the KEGG pathway analysis. SaMYB5 corresponds to CDC5 (K12850) of the spliceosome pathway (Ko03040). In addition, we screened the SaMYBG gene, which showed negative regulation under drought and alkali conditions. Previous results suggested that overexpression of soybean MYB transcription factors GmMYBJ1, GmMYBJ2, GmMYB12B2, GmMYB3a, and GmMYB68 increased the tolerance of Arabidopsis to salt, alkali, and drought [57-61]. AtMYB11, AtMYB12, and AtMYB111 of Arabidopsis thaliana have different regulatory effects on flavonoid biosynthesis [62-64]. We speculate that SaMYB14, SaMYB1, SaMYB5, and SaMYBG might also have similar regulatory functions under salt, alkali, and drought stress.

Many bZIP transcription factor genes have been identified in different plant species, such as cucumbers, corn, and legumes [65-67]. In plants, bZIP transcription factors are essential for various biological processes, such as seed maturation, stress signal transduction, and flower development [68]. Arabidopsis AtbZIP17 and AtbZIP28 are involved in regulating root elongation during stress response [69]. In this study, we screened 3 bZIP transcription factor genes: SabZIP1, upregulated under alkali stress; SabZIP3, upregulated under salt stress; and SabZIP5, downregulated under both salt and drought stress (Table S2). Furthermore, SabZIP1 was involved in ko00280, ko00270, ko00260 and ko00250 pathway.

WRKY transcription factors response to stress by regulating plant secondary metabolites, such as alkaloids, terpenes, and other subclasses [70, 71]. In this study, three WRKY transcription factor genes were obtained: SaWRKY27, SaWRKY33, and SaWRKY38 (Table S2). Both SaWRKY27 and SaWRKY33 were upregulated under alkali stress and participated in the plant-pathogen interaction (Ko04626) pathway. SaWRKY38 was downregulated under both alkali and drought stress conditions. In common bean, 8 up-regulated WRKY transcription factors were confirmed under drought stress [71]. Most WRKY transcription factors obtained in Dunaliella bardawil can bind to the W-box and response to n abiotic stress [72]. In peanut, 73 differentially expressed AhWRKY genes affected by drought stress were obtained through expression pattern analysis of the WRKY gene family [73]. SbWRKY30 enhanced the drought tolerance of sorghum by regulating the drought stress response gene SbRD19 [74]. The overexpression of the wheat WRKY transcription factor gene TaWRKY13 can significantly improve salt tolerance in rice, indicating that the WRKY transcription factor was a positive factor responding to salt stress [75].

We analyzed NAC transcription factors and found that SaNAC2, SaNAC5, and SaNAC25 were upregulated under alkali stress and drought stress (Table S2). Studies in celery (Apium graveolens L.) have shown that AgNAC63 (a homologous gene of ANAC072/RD26) is highly induced under the conditions of heat, cold, and salt [76]. Through transcriptome analysis of 4 cotton varieties, it was determined that 120 NAC transcription factor genes may be involved in response to salt, alkali and drought stress [77]. MfNACs in Medicago falcate maintain glutathione levels by regulating the gene expression of glyoxalase 1, which is a defense response against drought stress [78].

Idenfication of the differentially expressed genes (DEGs)

To confirm that the DEGs obtained by Illumina Hiseq 2000 and 2500 platform sequencing were credible, 24 DEGs (13 up-regulated and 11 down-regulated) were randomly chosen for qRT-PCR analysis (Fig. 5 and Table S4). The results showed that expressions of these DEGs were similar to those obtained by RNA-seq. The results indicate that the method used to confirm DEG in this experiment is feasible.

Fig. 5

The relative expression levels of representative DEGs between control and stress-treatment samples, sequenced for Sophora alopecuroides. Actin was used as internal reference. Relative transcription levels were calculated using the 2−ΔΔCt method. Data represent means ± standard deviation (SD) from three biological replicates and three technical replicates. RNA-seq value is base on fold chang of upregulation or downregulation DEGs. The relative expression level indicates the fold change obtained by quantitative RT-PCR. Values are the means ±SD. Means were generated from three independent replicates. Statistical comparisons (one-way ANOVA) are presented for each variable (** P < 0.01; * P < 0.05)

Discussion

Sophora alopecuroides (known as ku dou zi in Chinese medicine) is famous for its tolerance to high salinity, alkalinity and drought. The genome of S. alopecuroides has not been sequenced previously. We therefore analyzed the transcriptome sequence of S. alopecuroides under different stress treatments, explored the changes in DEGs and metabolic pathways under salt, alkali, and drought stress, and searched for genes related to abiotic-stress resistance. This study focuses on the expression of genes in the roots of Sophora alopecuroides, which has some limitations. Many genes (DEGs) were expressed differentially between the treatments and control. We found many transcription factor genes, including MYB, bZIP, NAC, ERF and WRKY family members, which were functionally annotated in all three treatments. Transcription factor families (MYB, WRKY, NAM, NAC, ATAF, AP2-EREBP, bHLH, C2H2, HSF, ERF, and bZIP) have been found responsible for abiotic stress in rice [79-84], maize [82], A. thaliana [83], soybean [84] and tobacco; overexpression of these TFs enhanced stress tolerance in transgenic plants [79, 83]. In this study, 4 MYB TFs (SaMYB14, SaMYB1, SaMYB5, and SaMYBG), 3 bZIP TFs (SabZIP1, SabZIP3, and SabZIP5), 3 WRKY TFs (SaWRKY27, SaWRKY33, and SaWRKY38), and 3 NAC TFs (SaNAC2, SaNAC5, and SaNAC25) were screened out as potential candidate genes that respond to salt, alkali, and drought stress (Table S2).

Significant enrichment analysis showed that the DEGs, mainly those related to plant membrane binding and signal transduction, were basically consistent with those of cotton [26, 85, 86], Caragana korshinskii [87], and banana [88, 89]. In the salt and drought treatments, the most abundant DEG types were those involved in biosynthesis. In the alkali treatment, the most abundant DEGs were those involved in metabolic processes. The most abundant DEGs involved in molecular functions and producing cellular components were the same under all three treatments.

By cross-analysis of DEGs obtained under the three treatments, we found that only four of the up-regulated DEGs were shared among the three treatments. Among the down-regulated DEGs, those that were common to the three treatments were mainly related to plant signal transduction and metabolism, including the ABA signaling pathway, peroxisome and flavonoid biosynthesis. ABA, as an important plant hormone, plays a key role in plant growth and development, including seed germination, reproduction, and seed dormancy. ABA also plays an important role in plant responses to biological and abiotic stress, particularly salt, drought, and cold stress [26, 85–91]. Transcription factor ATHB6 is involved in down-regulating ABA responses, for instance during seed germination and stomatal closure, when plants are sensitive to ABA [92]. However, the up-regulated DEG SaATHB annotation obtained in this study is related to the homeobox-leucine zipper protein, reflecting differences in DEG responses to stress under different treatments. The results need to be validated prospectively.

The results of GO enrichment analysis showed that there were differences in DEGs and metabolic pathways between the salt-, alkali-, and drought-stress samples, indicating that the response mechanisms of plants to salt, alkali, and drought stress were not identical under different treatments. Hormone signal transduction pathway is a KEGG pathway directly related to plant metabolism. This pathway involves zeatin biosynthesis, tryptophan metabolism, diterpenoid biosynthesis, carotenoid biosynthesis. Nine DEGs were involved in carotenoid biosynthesis under salt treatment. Among these, five downregulated DEGs (SaPYL4–1, SaPYL4–2, SaPYL4–3, SaPYL4–4, and SaPYL5–1) were related to the pyrabactin resistance 1-like (PYR/PYL) ABA receptor, and four upregulated DEGs (SaPP2C8, SaPP2C16, SaPP2C37, and SaPP2C53) were related to protein phosphatase 2C (PP2C). Furthermore, five SaPYLs and two SaPP2Cs (SaPP2C8 and SaPP2C53) showed the same expression pattern between salt and alkali stress. This result indicates that there is a negative regulation between PYR/PYL and PP2C, which is same as previous studies that showed that PYR/PYL inhibits PP2C [34, 93]. Under drought treatment the expression of PYL-related genes in the ABA signaling pathway was different from that under saline-alkali stress. Four PYL-related genes, SaPYL4–1, SaPYL4–2, SaPYL5, and SaPYL9 were up-regulated under drought stress, but were down-regulated under salt and alkali. The finding indicated that the mechanisms of drought stress and saline-alkali stress responses is different in Sophora alopecuroides.

In addition, under salt treatment conditions, a down-regulated gene SaTCH4 was identified in the brassinosteroid biosynthesis pathway, specifically in xyloglucosyl transferase TCH4. TCH is related to plant perception of external stimuli and morphogenesis during adaptation. TCH is up-regulated when plants are stimulated by environmental stimuli [94]. However, the differential genes that correspond to TCH were down-regulated under salt stress in this experiment.

CRE1 has been identified as a cytokinin receptor in A. thaliana [95]. Cytokinin is an important plant hormone regulating cell division and differentiation [96]. In response to abiotic stress, plants may repair tissue damage, which may require cytokinins. SaAHK4 was found to be expressed both under alkali and drought stress, corresponding to the CRE1 metabolic pathway and the expression level was similar between alkali and drought treatments. CTR1 encodes a member of the protein kinase Raf family, and down-regulates the ethylene response pathway in A. thaliana [97]. There have been some reports about CTR1 involvement in stress responses. For instance, the expression of GmCTR1 in soybean infected with scab was up-regulated, peaked 48 h after infection, and subsequently decreased slightly; this indicates that CTR1 might have a role in stress-responses [98]. CTR1 overexpression in tobacco can make transgenic tobacco more sensitive to salt stress, suggesting that CTR1 may down-regulate responses to salt stress [99]. The study of Arabidopsis ctr1 mutants showed that Ctr1 may reduce the salt tolerance of Arabidopsis through negative regulation of ethylene signal [100-104]. Our results showed that nine genes corresponding to CTR1 were down-regulated under drought stress, which is consistent with previous findings about responses to abiotic stress processes. However, under alkali stress, we found that 17 up-regulated DEGs corresponded to CTR1. This indicates that S. alopecuroides uses different processes to cope with alkali and drought stress. Further, this analysis provides a new approach for studying the role of the CTR1 gene in alkali-stress responses in S. alopecuroides.

Under stress, plant tissues produce ROS, which can destroy membrane structure, damage biological macromolecules, cause metabolic disorders, and in severe cases cause plant death [105-110]. Superoxide dismutase, peroxidase, and catalase are protective enzymes in plants. Their main functions are to scavenge ROS free radicals, prevent their excessive accumulation, avoid or mitigate ROS attack on membranes, and prevent membrane damage [106, 111, 112]. In the KEGG enrichment pathway sequences, our salt, alkali, and drought treatments all reflected the presence of superoxide dismutase, peroxidase, and catalase. Through screening and comparison, we found some DEGs, including SaPXR1, SaSOD1, SaSOD2–1, SaSOD2–2, and SaHACL1, which code for superoxide dismutase, peroxidase, and catalase. This suggests that the physiological and biochemical processes and regulatory mechanisms differ depending on the type of abiotic stress.

Conclusions

Our study is the first to analyze DEG expression using no-reference transcriptome sequencing under salt and alkali drought stress in S. alopecuroides. We have identified many candidate functional genes that are involved in multiple salt, alkali, and drought stress tolerance mechanisms, and these should be studied further. We found that key pathways, such as those response to plant hormone signal transduction, biosynthesis of secondary metabolites, metabolic processes, peroxisome production, and cellular morphology, were involved in abiotic-stress tolerance in this species. Many genes (DEGs) were expressed differentially between the treatments and the control. Our findings may assist in molecular design breeding about molecular adaptations to extreme environment and will provide resources for the study of plant stress resistance.

Methods

Plant materials and treatment

The 10 g seeds of S. alopecuroides were collected from Korla region, Xinjiang Autonomous Region, China (41° 45′ N, 86° 8′ E). After treatment using 98% H2SO4, the seeds were dipped in sterile water for 12 h and planted in pots filled with soil and vermiculite (1:1), and kept in day/night of 16 h/8 h, 25 °C /22 °C. Four-week-old seedlings were divided into 4 groups: control (CK), salt treatment (S_T), alkali treatment (AS_T) and drought treatment (DT), and each treatment was replicated three times. For the salt, alkali, and drought treatments, 1.2% NaCl, 1.2% NaHCO3 and 8% PEG6000, respectively, were added to 1/8th strength Hoagland’s nutrient solution. After 72 h of treatment, the root tissues of the control group and the treatment group were taken, and stored at − 80 °C.

RNA extraction, quantification, and qualification

The extraction and isolation of RNA used the previous methods [57], with 3 replicates for each group. Then use the 2100 bioanalytical instrument to evaluate the quality of the obtained RNA.

cDNA library preparation for transcription sequencing

Following the manufacturer’s recommendations of RNA Library Preparation Kit (New England Biolab, Massachusetts, USA), a sequencing library was constructed. Magnetic beads with poly-T oligonucleotides attached were used to enrich the mRNA. At a high temperature, divalent cations were used to carry out the cleavage. First strand of cDNA was synthesized using random hexamer primers. cDNA library was prepared followed the previous research [113]. The library was estimated and sequenced on the Illumina HiSeq platform (New England BioLabs, MA, USA) [113].

Quality control, assembly, and annotation of unigenes

Adapter, ply-N and reads with low quality were removed to generated clean data, which was performed in downstream analysis. Transcriptome was assembed using Trinity by Beijing Novogene Company with reference to previous research methods [114]. Annotation of gene function using the BLASTX algorithm (E value < 1.0 E− 5) is based on the following databases: Nt, Nr, GO, Pfam, Swiss-Prot, KOG/COG, and KO. To obtain effective unigenes, all assembled unigenes were screened against multiple databases. Using Nr annotation, unigenes were assigned to BP, MF, and CC ontologies according to the GO framework using Blast2GO software [115]. Pathway distributions were perfomed on the basis of KEGG database using the BLASTX algorithm (E value < 1.0 E− 5).

Analysis of differentially expressed genes (DEGs)

The transcriptome obtained using the Trinity platform provided a reference sequence (hereafter “ref”). RSEM software was used to map clean reads for every sample on ref. [116]. The results were calculated using RSEM, and the read count of each gene was procured. FPKM was used to determine the number and expression level for each unigene. The input data of DEG was the read count data obtained from the analysis of gene expression levels [117]. We used the DESeq method to analyze three biological replicate samples, and the screening threshold was padj < 0.05 [118]. Negative binomial models were used to calculate the p-value. GO and KEGG pathway enrichment analysis were used to analyze the DEGs. GOseq5 GO enrichment analysis method is based on Wallenius’ noncentral hypergeometric distribution [115]. KOBAS software was used to examin the effective enrichment of KEGG pathway [119].

Verification of sequencing results by real-time quantitative RT-PCR

To verify the reliability of the sequencing results, 24 unigenes were randomly selected for expression analysis. Total RNA (5 μg) from each of three replicates was extracted and the primers were listed in Table S4. To perform qRT-PCR, 20 μL of reaction mixture was placed in a 96-well plate in the StepOnePlus™ Real-Time PCR System using SYBR Premix Ex Taq™ II kit. Each sample had three independent biological and technical replications. The expression level of the actin gene (J01298) was used to normalize target gene expression. Relative transcription levels were calculated using the double delta CT (2−ΔΔCt) method.

Additional file 1: Figure S1. Phenotypes of Sophora alopecuroides at different concentration in salt, alkali and PEG. This figure provides the phenotype under salt, alkali and drought treatment (The figures were taken by authors themselves).

Additional file 2: Table S1. Statistics of the number of DEGs corresponding to KEGG metabolic pathways. This table provides the number of DEGs corresponding to KEGG metabolic pathways.

Additional file 3: Table S2. Statistics of DEGs corresponding to KEGG metabolic pathways. This table provides the potential genes corresponding to KEGG metabolic pathways.

Additional file 4: Table S3. Statistics of the number of transcription factors corresponding to each transcription factor family. This table provides the number of transcription factors corresponding to each transcription factor family.

Additional file 5: Table S4. Primer sequences validated by qRT-PCR in this study. This table provides the primer sequences that were validated by qRT-PCR in this study.