We identified an alkalinity-tolerant wild soybean variety N24852 from a preliminary screening of 129 soybean varieties (our unpublished data), which showed more alkalinity-tolerant than 95% of the 129 varieties. N24852 is tolerant to biotic and abiotic stresses (Wang et al., 2013) and its alkalinity tolerance was confirmed in this study (Supplementary Figure S1). A widely used salinity-tolerant (but moderately sensitive to alkalinity) cultivated soybean variety Lee 68 (Abel and Mackenzie, 1964), was used as a genotype control to evaluate the alkalinity tolerance of N24852. These two varieties were grown in a growth chamber (E-41HO, Percival, USA) with 60% relative humidity, 15 h light (50000 lx)/9 h darkness photoperiod and corresponding temperature regime of 28°C/24°C. These two varieties showed obvious phenotypic difference when using 90 mM NaHCO₃ treatment, but only showed subtle difference at 80 mM NaHCO₃ and no difference at 100 mM NaHCO₃ treatment in the preliminary experiment (Supplementary Figure S1). Therefore, 90 mM NaHCO₃ was used for RNA-seq study.

Soybean seeds were surface-sterilized with 1% sodium hypochlorite for 30 s, and rinsed three times with deionized water. Twenty seeds were sown in plastic pots (Φ10 × 8 cm) which filled with clean quartz sand. Seven days after germination, the seedlings were thinned to four plants per pot. Four pots were placed in a plastic container (34 cm × 24 cm × 8 cm) containing 1.5 L fresh 1/2 strength Hoagland nutrient solution (pH ≈ 6.5), which could be absorbed through small holes at the bottom, and the solution was changed every 2 days. When the second trifoliolate leaves appeared (V3, approximately 14 days after planting), a treatment solution containing 90 mM NaHCO₃ (with a final pH of 8.5, adjusted by addition of KOH) was applied to induce alkaline stress. The control contained no NaHCO₃. The experiment included three independent biological replications.

Before harvest, roots were dipped into an iso-osmotic solution of 10 mM Ca(NO₃)₂ for 10–20 s to avoid turgor loss and rapid efflux of ions from the apoplast and epidermal cells due to osmotic shock (Munns et al., 2010), and rinsed three times with deionized water. Root tips (3 cm) of NaHCO₃-treated and control plants were harvested after 12 and 24 h treatment, then immediately frozen in liquid nitrogen and stored at -80°C for RNA isolation. The remaining root tissues and leaves from the stressed and control plants were harvested at the same time points for measurement of ion concentration. For determination of ion concentrations, samples were pretreated at 105°C for 30 min to deactivate enzymes and dried at 80°C for 5 days. In order to minimize variation between individual plants, bulked roots and leaves were harvested from four uniform plants in a single pot for each biological replicate.

Methods for ion extraction and measurement followed Munns et al. (2010). Dried samples were ground into powder. Then 0.1 g of the powder was digested in 2 mL nitric acid (HNO₃) in a microwave digestion system (ETHOS T, Milestone, Italy) following the manufacturer’s recommendations. Na⁺, K⁺, and Ca²⁺ concentrations were estimated using an Inductively Coupled Plasma Optical Emission Spectrometer (ICP-OES) Optima 8000 (PerkinElmer, USA) according to the manufacturer’s instructions. The mg/L values were converted to mg/g dry weight (DW) using the formula of (concentration in mg/L as measured × dilution factor × volume of dilute nitrate acid)/(DW of tissue used in extraction × 1000). Differences between two soybean genotypes were compared by t-tests using the SAS proc ttest.

Roots of wild soybean N24852 under 90 mM NaHCO₃ and control at 12 and 24 h were harvested for RNA extraction and RNA-seq. A total of 12 samples (2 time points × 2 treatments × 3 biological replicates) were prepared for RNA-seq. Total RNA was extracted using Trizol^® reagent (Invitrogen, USA). Total concentration and quality of the RNA samples were determined by agarose gel electrophoresis and a NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific, USA). An Agilent 2100 Bioanalyzer RNA Nano chip (Agilent Technologies, USA) was used for accurate quantification.

Messenger RNA (mRNA) was isolated from the total RNA by magnetic beads coated with oligo (dT) using a Dynabeads mRNA DIRECT Kit (Invitrogen, USA). Approximately 5 μg of mRNA from each sample was prepared to construct a library. The mRNA was fragmented into small pieces with fragmentation buffer at 90°C. Purified mRNA fragments were reverse transcribed (Invitrogen, USA) with random hexamer adaptors to synthesize the first strand cDNA. Second strand cDNA was synthesized with RNaseH (Invitrogen, USA) and DNA polymerase I (Invitrogen, USA). The cDNA was cleaned using Agencourt Ampure XP SPRI beads (Beckman Coulter, USA). The synthesized cDNAs were appended with an ‘A’ base at the 3′-end and ligated with paired-end adaptors. Paired-end cDNA libraries were amplified using PCR, and fragments of 400–500 bp insertions were selected from 2% agarose gel electrophoresis separation and quantified using qPCR. The twelve cDNA libraries were sequenced on an Illumina HiSeq 2500 platform (Illumina Inc. USA) according to the manufacturer’s recommendations at Berry Genomics Company, Beijing, China.

The sequence data (in FastQ format) have been submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) databases with the study accession SRP093892. The raw reads were cleaned by removing adaptor reads and low-quality reads (ambiguous sequences with ‘N’ percentage values ≤3 and the percentage of low-quality bases less than 3 is ≥50%) using an in-house script written in C. The clean data were used for subsequent data analysis. High-quality clean reads were mapped to the latest version (G. max Wm82.a2.v1) of soybean reference genome (Goodstein et al., 2012) downloaded from phytozome¹, using Tophat2 (v2.0.13) software (Kim et al., 2013). Mismatches of less than or equal to 2 bp were allowed in mapping to reference sequences. Subsequently, the sequencing saturation of the library, coverage analysis of clean reads on reference genes, as well as the genomic distributions in CDS (exons), introns, and intergenic regions were analyzed.

Gene expression levels were normalized to Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) by Cufflinks v2.2.2 software (Trapnell et al., 2010). DEGs were identified by comparing the NaHCO₃-treated and control samples at the same time point using the R package DESeq (Anders and Huber, 2010). The significant DEGs were identified using the false discovery rate (FDR) ≤ 0.01 and | log₂FoldChange|≥ 1.

All mapped unigenes and new genes found in this study were compared against seven databases including National Center for Biotechnology Information (NCBI) non-redundant protein database (NR), Gene Ontology (GO), Clusters of Orthologous Groups of proteins (COG) (Tatusov et al., 2000), euKaryotic Ortholog Groups (KOG) (Koonin et al., 2004), reviewed protein sequence database (Swiss-Prot) (Apweiler et al., 2004), and Kyoto Encyclopedia of Genes and Genomes (KEGG) Ontology (KO) (Kanehisa et al., 2004) using BLASTn (v 2.2.26) software (Altschul et al., 1999) with an E-value cutoff at 10^-5, and searched against the Protein family (Pfam) database (Finn et al., 2014) by hmmscan (v 3.0) software (Johnson et al., 2010). All genes were also BLAST against the Plant Transcription Factors Database v3.0 (Jin et al., 2014) with an E-value cutoff at 10^-5. The WEGO software package (Ye et al., 2006) was used for describing GO functional classification of cellular component, molecular function and biological process. GO enrichment analyses of DEGs were performed using Singular Enrichment Analysis (SEA) method with P < 0.01 and FDR < 0.05 by agriGO (Du et al., 2010), and the newest soybean genome Wm82.a2.v1 was set as background. The hypergeometric Fisher exact test (P < 0.01) and Benjamini and Hochberg method (FDR < 0.05) was performed to detect statistically significant enrichment of KEGG pathway and transcription factors, in comparison with the whole soybean transcriptome as the background. DEGs were also blasted against the KOG databases² for functional classification.

Heat map of the DEGs overlapped between two time points (12 and 24 h), were analyzed using heatmap.2 function of the R/Bioconductor package gplots with default options (Warnes, 2016). The gene expression levels were transformed by log₂ (FPKM+1) using three biological replications.

To validate the RNA-seq gene expression patterns, 14 DEGs were selected and investigated by qRT-PCR using the same RNA samples for the RNA-seq library construction. We used the housekeeping gene GmUKN1 (Glyma.12G020500) as an internal control to normalize the level of gene expression (Hu et al., 2009; Guan et al., 2014). Primer pairs were designed for qRT-PCR using the online website NCBI primer-BLAST³. Primer sequences and gene annotations are listed in Supplementary Table S1. First-strand cDNA was synthesized from 500 ng/10 μl total RNA by using a PrimeScript RT Reagent Kit (TaKaRa, Japan). qRT-PCRs were performed in a LightCycler^® 480 II (Roche, Germany) in a final volume of 25 μl containing 12.5 μl SYBR Premix Ex Taq II (Tli RNaseH Plus) (TaKaRa, Japan), 2 μl (200 ng) of cDNA, 1 μl (10 mM) of the forward and reverse primers, and 8.5 μl of ddH₂O. PCR conditions were set at 95°C pre-denaturation for 3 min and followed by 40 cycles of 95°C denaturation for 15 s, 60°C annealing for 15 s, and 72°C extension for 15 s. Relative gene expression level was calculated using the 2^-ΔΔCT method (Livak and Schmittgen, 2001). All qRT-PCR were performed in three technical replicates. The correlation of RNA-seq data with qRT-PCR analysis was calculated using SAS (v9.3) proc corr based on log₂ fold changes.

Maintenance of a low Na⁺ concentration, and low ratios of Na⁺/K⁺ and Na⁺/Ca²⁺ are widely used as indices of plant salinity tolerance (Puranik et al., 2011; Yong et al., 2014; Tuyen et al., 2016). To evaluate the alkalinity tolerance of wild soybean N24852, the salt tolerant soybean variety Lee 68 was used for comparison. The wild soybean N24852 showed more tolerance to 90 mM NaHCO₃ stress than Lee 68, as indicated by the later appearance of chlorosis and wilting (Supplementary Figure S1). The concentration of Na⁺, ratios of Na⁺/K⁺ and Na⁺/Ca²⁺ in the roots of both soybean varieties increased after NaHCO₃ treatment (Figure 1). But the wild soybean N24852 showed significantly (P < 0.05, student’s t-tests) lower concentration of Na⁺, ratios of Na⁺/K⁺ and Na⁺/Ca²⁺ in the roots than Lee68 under NaHCO₃ stress (Figure 1). And the Na⁺ concentration, ratios of Na⁺/K⁺ and Na⁺/Ca²⁺ in the leaves of N24852 were also lower than those of Lee68 (Figure 1).

Due to the lower Na⁺ concentration in the roots of N24852 than that of Lee 68 (Figure 1) was observed as early as 12 h after alkalinity, RNA-seq analyses of N24852 roots after 12 and 24 h of NaHCO₃ stress in comparison with control were performed on three independent biological replicates. A total of 12 cDNA libraries were constructed and sequenced on Illumina HiSeq 2500, and 283.6 million raw reads were generated (Supplementary Table S2). After removing adaptors and low quality sequences, a total of 270.5 million (95.38% of the raw reads) clean reads (approximately 67.63 Gb clean data) were obtained. On average, 23.6 million clean reads (5.64 Gb clean data) were obtained from each sample (Supplementary Table S2). The percentages of Phred-like quality scores at the Q30 level (an error probability of 1‰) ranged from 87.24 to 94.45% and the average GC content was estimated as 46.43% (Supplementary Table S2). Among the 12 samples, 88.16–91.69% of the clean reads were mapped to the reference genome, and 89.85–96.60% of clean reads were uniquely mapped (Supplementary Table S2). The saturation curves of 12 RNA-seq samples (genes with FPKM ≥ 0.01) estimated that the gene coverage started to show saturation when approximately more than 5 million clean reads were aligned (Supplementary Figure S2A). The average clean reads of our 12 samples were 22.54 million, which is more than the saturation threshold. Gene coverage analysis indicates that the sequencing reads were uniformly distributed from the 5′ to 3′ of genes (Supplementary Figure S2B). On average, more than 80% of the mapped reads were located at the exon region (Supplementary Figure S3). Detailed information of the RNA-seq data was listed in Supplementary Table S2, suggesting that the sequencing quality was high and sequencing depth was sufficient for transcriptome coverage.

A total of 54,844 unigenes were detected in our transcriptome by cufflinks program, including 54331 genes that were aligned to the soybean reference genome (96.94% of the total 56044 genes in Williams 82) and 513 new genes. Among these, 54342 genes (including 53888 mapped unigenes and 454 new genes) were annotated by at least one of the seven databases (Supplementary Table S3). Taken together, 96.08% (54342) of the expressed genes were successfully annotated in at least one of databases, with 10.51% (5946) of genes annotated in all databases (Supplementary Table S3). The 54342 annotated genes are listed in Supplementary Table S4.

By using the criteria of FDR ≤ 0.01 and | log₂FoldChange|≥ 1, a total of 449 genes were differentially expressed in wild soybean roots under NaHCO₃ stress compared with control (Supplementary Table S5). As shown in the Venn diagram (Figure 2), there were 95 and 140 up-regulated genes and 108 and 135 down-regulated genes after 12 and 24 h of alkaline stress, respectively, with more DEGs detected after 24 h alkaline stress than 12 h, implying more transcriptional changes at 24 h. The numbers of up- and down-regulated genes were similar. Nine genes were up-regulated and 20 genes were down-regulated at both time points (Figures 2 and 3), suggesting their importance in soybean response to NaHCO₃ Stress. The overlapped up-regulated genes (Figure 3) included two aluminum-activated malate transporter (ALMT) genes (Glyma.11G179100 and Glyma.12G094400) and one gene (Glyma.13G363300) encoding a late embryogenesis abundant protein (LEA). ALMT transporters have been reported to play roles in adaptation of plants to abiotic stress (Henderson et al., 2014), and LEA proteins accumulate in response to water deficit caused by drought, heat and salinity (Battaglia and Covarrubias, 2013; Kosova et al., 2014). No DEG showed opposite regulation patterns at two time points (Figure 2; Supplementary Table S5).

By WEGO database, the DEGs were classified into 12, 9, 12 categories for cellular components, molecular functions, and biological processes, respectively, after 12 h of NaHCO₃ stress, and were classified into 12, 13, and 16 categories, respectively, at 24 h of NaHCO₃ treatment (Supplementary Figure S4). By KOG database, there were 97 and 138 DEGs with KOG functional classification information at 12 and 24 h of NaHCO₃ stress, respectively, and they were classified into 17 functional categories (Supplementary Figure S5).

We also performed GO enrichment analyses of DEGs using agriGO (P < 0.01, FDR < 0.05). At 12 h of alkaline stress, genes with transcription factor activity (GO:0003700) were significantly enriched in up-regulated genes (Figure 4A), while genes with peptidase activity (GO:0008233) and related to biological process of proteolysis (GO:0006508) were enriched in down-regulated genes (Supplementary Figure S6A). At 24 h of alkaline stress, genes with the molecular function of transporter activity (GO:0005215) and carbohydrate binding (GO:0030246), and involved in biological process of transport (GO:0006810) were significantly enriched in up-regulated genes (Figure 4B), while GO terms in molecular functions of “hydrolase activity (GO:0016787),” “transferase activity (GO:0016757),” and “oxidoreductase activity (GO:0016705),” biological process of “carbohydrate metabolic process (GO:0005975)” especially “cellular glucan metabolic process (GO:0006073)” were enriched in down-regulated genes (Supplementary Figure S6B). Consistent with the number of DEGs, the numbers of enriched GO terms at 24 h were greater than 12 h after alkaline stress.

All DEGs were BLAST against the KEGG Ontology (KO) database. At 12 h of NaHCO₃ stress, 30 DEGs were classified into 29 biological pathways belonging to 13 KEGG categories (Supplementary Table S6). At 24 h after NaHCO₃ stress, 48 DEGs were classified into 37 biological pathways belonging to 12 KEGG categories (Supplementary Table S7).

We also used the Fisher exact test to analyze KEGG enrichment for DEGs (P < 0.01, FDR < 0.05). No enriched biological pathways were found at 12 h after NaHCO₃ stress, but enrichment of “phenylpropanoid biosynthesis” and “phenylalanine metabolism” was found at 24 h after 90 mM NaHCO₃ (pH = 8.5) treatment (Supplementary Figures S7 and S8). The genes involved in both pathways were peroxidase-encoding genes, including three up-regulated genes (Glyma.07G209900, Glyma.17G053000, and Glyma.20G169200) and one down-regulated gene (Glyma.02G171600). As one of the important antioxidant enzymes, peroxidases (PODs) could protect plants from oxidative damage by scavenging of reactive oxygen species (ROS), which is induced by various abiotic and biotic stresses (Miller et al., 2008; Gill and Tuteja, 2010).

To validate the expression data of RNA-seq, we selected 14 DEGs for qRT-PCR analysis, including four POD-encoding genes, three up-regulated genes and two down-regulated genes at both 12 and 24 h, and five genes that were up- or down- regulated at only one time point (Figure 5; Supplementary Table S1). To compare the expression data between RNA-seq and qRT-PCR, the relative expression level was transformed to log₂ fold change. The qRT-PCR results showed a high consistency (linear regression equation y = 0.9142x + 0.0851, R² = 0.9763) with RNA-seq data (Figure 5), indicating the reliability of RNA-seq expression profile in this study.

Transcription factors are essential for regulation of gene expression, by binding to the specific cis-acting elements in the genes that they regulate. A total of 130 and 173 transcription factors representing 36 and 38 different families were differentially expressed at 12 and 24 h, respectively, under NaHCO₃ stress comparing with control in wild soybean roots (Figure 6). Among the differentially expressed transcription factors, high percentages of basic helix-loop-helix (bHLH), ethylene-responsive factor (ERF), Trihelix, and zinc finger (C2H2, C3H) families were found at both time points. Nuclear Factor Y subunit A (NF-YA) family was enriched at 12 h after NaHCO₃ stress. These transcription factors have already been reported to play roles in plant response to abiotic stress (Ge et al., 2010; Kielbowicz-Matuk, 2012; Li et al., 2013; Yong et al., 2014; Dey and Vlot, 2015). No transcription factor family was enriched at 24 h after NaHCO₃ stress.

ABC transporter (Lee et al., 2004; Henderson et al., 2014), ALMT family (Barbier-Brygoo et al., 2011; Henderson et al., 2014), glutamate receptor (GLR) (Demidchik et al., 2002), nitrate transporter (NRT)/proton dependent oligopeptide (POT) family (Li et al., 2010; Barbier-Brygoo et al., 2011; Chen et al., 2012), and S-type anion channel (SLAH) (Negi et al., 2008; Barbier-Brygoo et al., 2011) have been reported to play roles in ion transport, which indicates their importance in maintaining ion homeostasis in plant response to NaHCO₃ stress. Nine transporter genes were differentially expressed at 12 h after NaHCO₃ stress (Table 1). At 24 h after NaHCO₃ treatment, the up-regulated genes were enriched in transport related genes as shown by GO enrichment analysis (Figure 4B). Twenty-four transporter genes were up-regulated and three transporter genes were down-regulated at 24 h after NaHCO₃ stress (Table 1), including the genes encoding ABC transporter, ALMT, ammonium transporter, aquaporin transporter, cation efflux family, GLR, NRT/POT family, SLAH, sulfate transporter, and zinc/iron transporter. Our results indicate the importance of these transporter genes in maintaining ion homeostasis in soybean response to NaHCO₃ stress.