Enhanced reactive oxygen detoxification occurs in salt-stressed soybean roots expressing GmSALT3.

Soybean (Glycine max) is an important crop globally for food and edible oil production. Soybean plants are sensitive to salinity (NaCl), with significant yield decreases reported under saline conditions. GmSALT3 is the dominant gene underlying a major QTL for salt tolerance in soybean. GmSALT3 encodes a transmembrane protein belonging to the plant cation/proton exchanger (CHX) family, and is predominately expressed in root phloem and xylem associated cells under both saline and non-saline conditions. It is currently unknown through which molecular mechanism(s) the ER-localised GmSALT3 contributes to salinity tolerance, as its localisation excludes direct involvement in ion exclusion. In order to gain insights into potential molecular mechanism(s), we used RNA-seq analysis of roots from two soybean NILs (near isogenic lines); NIL-S (salt-sensitive, Gmsalt3), and NIL-T (salt-tolerant, GmSALT3), grown under control and saline conditions (200 mM NaCl) at three time points (0 h, 6 h, and 3 days). Gene ontology (GO) analysis showed that NIL-T has greater responses aligned to oxidation reduction. ROS were less abundant and scavenging enzyme activity was greater in NIL-T, consistent with the RNA-seq data. Further analysis indicated that genes related to calcium signalling, vesicle trafficking and Casparian strip (CS) development were upregulated in NIL-T following salt treatment. We propose that GmSALT3 improves the ability of NIL-T to cope with saline stress through preventing ROS overaccumulation in roots, and potentially modulating Ca2+ signalling, vesicle trafficking and formation of diffusion barriers.

INTRODUCTION

Plants use reactive oxygen species (ROS) as signalling molecules at low concentrations to regulate various biological processes, such as growth, programmed cell death, hormone signalling and development (Foreman et al., 2003; Mittler, 2002; Neill et al., 2002; Overmyer et al., 2003; Pei et al., 2000). However, ROS are also a toxic by‐product of many reactions, which subsequently need to be eliminated by ROS scavenging enzymes. Under non‐stressed conditions the production of ROS and the capacity of the cell to scavenge them are generally at equilibrium (Foyer & Noctor, 2005). However, under biotic and abiotic stresses, such as salinity, drought, temperature extremes, flooding, heavy metals, nutrient deprivation, and pathogen attack, intracellular ROS levels can increase dramatically. ROS at high concentrations damage plant cells through lipid peroxidation, DNA damage, protein denaturation, carbohydrate oxidation, and enzymatic activity impairment leading to significant damage to cellular functions and even cell death (Foyer & Noctor, 2005; Gill & Tuteja, 2010; Mittler et al., 2004; Noctor & Foyer, 1998). ROS damage due to environmental stress is therefore a leading factor contributing to reduced global crop production. ROS have, for example, detrimental effects on membrane integrity under drought stress and cause early senescence, lowering crop yield (Mittler, 2002; Sharma et al., 2017). In plants, ROS are eliminated through enzymatic and non‐enzymatic pathways. Enzymatic antioxidant systems include scavenger enzymes such as SOD (superoxide dismutase), APX (ascorbate peroxidase), GPX (glutathione peroxidases), PrxR (proxiredoxin), GST (glutathione‐S‐transferase), and CAT (catalase); while non‐enzymatic pathways include low molecular metabolites, such as ASH (ascorbate), GSH (glutathione), proline, α‐tocopherol (vitamin E), carotenoids and flavonoids (Choudhury et al., 2017; Mittler et al., 2004).

Among the above‐mentioned stresses, salinity is one of the most prominent factors restricting crop production and agricultural economic growth worldwide. More than US$12 billion in revenue is estimated to be lost annually because of saline‐affected agricultural land areas (Bose et al., 2014; Flowers et al., 2010; Gilliham et al., 2017). Soybean (Glycine max [L.] Merrill) is an important legume crop that contributes 30% of the total vegetable oil consumed and 69% of human food and animal feed protein‐rich supplements (Lam et al., 2010; Prakash et al., 2001). The yield of soybean can be significantly reduced by salinity stress, especially during the early vegetative growth stage (Pi et al., 2016).

In soybean, GmSALT3 (also known as GmCHX1 and GmNcl) has been identified as a dominant gene that confers improved salinity tolerance (Do et al., 2016; Guan et al., 2014; Qi et al., 2014). It is mainly expressed in root phloem‐ and xylem‐associated cells and the protein is localised in the ER (endoplasmic reticulum), not the plasma membrane (PM) (Guan et al., 2014). NILs (near isogeneic lines), differing through a naturaly occuring transposon insertion in the GmSALT3 gene, have been used to show that GmSALT3 is involved in shoot exclusion of both Na+ (sodium) and Cl− (chloride). In NILs with a functional GmSALT3, leaf Cl− exclusion can be observed prior to Na+ exclusion (Liu et al., 2016), suggesting two distinct mechanisms for the exclusion of the two ions. Reciprocal grafting of NILs showed that shoot Na+ exclusion occurs via a root xylem‐based mechanism and shoot Cl− exclusion likely depends upon novel phloem‐based Cl− recirculation (Qu et al., 2021). Additionally, GmSALT3 is also related to K+ homeostasis (Do et al., 2016; Qu et al., 2021). Despite the known importance of GmSALT3 for soybean salinity tolerance, the molecular mechanism through which the ER‐localised GmSALT3 contributes to improved salinity stress tolerance has not yet been revealed.

High‐throughput “‐omics” technologies, including transcriptomics, proteomics, and metabolomics have been applied in an attempt to understand soybean root responses to salinity stress (Aghaei et al., 2009; Ge et al., 2010; Lu et al., 2013; Pi et al., 2016; Qin et al., 2013; Toorchi et al., 2009). While these studies provide a basis for examining the general responses of soybean roots to salt stress, they compare genetically diverse soybeans that differ in many aspects. For instance, proteomic studies utilising the cultivar Wenfeng07 (relatively salt‐tolerant) and Union85140 (relatively salt‐sensitive), suggested that Wenfeng07's tolerance to salinity stress was associated with flavonoid accumulation. Wenfeng07 contains the GmSALT3 salt tolerant allele and it is absent in Union85140, and flavonoids belong to the metabolites that scavenge ROS. Several key enzymes in flavonoid synthesis including chalcone synthase (CHS), chalcone isomerase (CHI) and cytochrome P450 monooxygenase (CPM) were differently expressed in Wenfeng07 (Pi et al., 2016), suggesting a possible link between ROS detoxification and salinity tolerance in these two diverse soybean cultivars. However, less genetically diverse material is needed to demonstrate an explicit and direct link between ROS detoxification and GmSALT3.

We took advantage of our previously isolated NILs, and performed RNA‐sequencing of NIL‐T (salt‐tolerant, GmSALT3) and NIL‐S (salt‐sensitive, Gmsalt3) roots, to further investigate the mechanism by which GmSALT3 confers salinity tolerance. Genes connected to ROS signalling were significantly differently regulated in NIL‐T compared to NIL‐S under saline conditions, suggesting a direct connection of ROS production to GmSALT3. Enzymatic assays further revealed that presence of GmSALT3 is essential for maintaining ROS homeostasis under saline conditions. Differential gene expression under control conditions suggested that NIL‐S might be primed for responding to biotic stress in comparison to NIL‐T.

MATERIALS AND METHODS

Plant growth conditions and stress treatments

NIL‐T (salt‐tolerant, GmSALT3) and NIL‐S (salt‐sensitive, Gmsalt3) plants were grown in a growth chamber (RXZ‐500D; Ningbo Jiangnan Instrument), with a day length of 16 h (with a light‐emitting diode light source at 400 μmol m−2 s−1) at 28°C, and 8 h dark at 25°C, with 60% relative humidity throughout. Soybean seedlings were treated with 200 mM NaCl (salt treatment) or water (control) at 10 days after sowing (DAS). 200 mM NaCl or water were applied again at 12 DAS.

Total RNA extraction and RNA‐seq library construction

Total RNA was extracted from soybean roots using TRIZOL reagent (Ambion, http://www.ambion.com). Whole roots were harvested at three time points, 0 h, 6 h, and 3 days of a 200 mM salt‐treatment with the corresponding non‐treatment controls. All the harvested soybean roots did not have nitrogen‐fixing nodules and were not inoculated with symbiotic nitrogen fixing bacteria. To remove the residual DNA, the extracted RNA was treated with RNase‐free DNase I (New England Biolabs, https://www.neb.com) for 30 min at 37°C. Thirty RNA libraries were generated for paired‐end reads using an Illumina HiSeq 2500 sequencer, consisting of three biological samples per time point per genotype.

RNA‐seq data analysis and assembly

Raw reads were generated and mapped to the latest soybean genome sequence Gmax_275 Wm82.a2.v1 (Glyma 2.0) using TopHat2 (Kim et al., 2013). Clean mapped reads were obtained by removing low quality (Q < 30) sequences, adapter fragments and barcode sequences. A quality control test by fastQC (Andrews, 2010) revealed the quality of RNA‐seq libraries construction and sequence alignment were sufficient for further analysis.

Gene expression level and DEG analysis

Gene expression was compared in FPKMs (Fragments Per Kilobase of transcript per Million) calculated using the Cufflinks functions cuffquant and cuffnorm (Trapnell et al., 2012). Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). The resulting p‐values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. The cutoffs for differentially expressed genes (DEGs) were Log2FC (fold change) ≥ 1, FDR (false discover rate) < 0.01.

GO enrichment and KEGG pathway analysis

GO enrichment analysis of DEGs was implemented using AgriGO (Du et al., 2010). GO terms with corrected p values <0.05 were deemed significantly enriched. KEGG (Kanehisa et al., 2007) is a database resource for understanding high‐level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular‐level information, especially large‐scale molecular datasets generated by genome sequencing and other high‐throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS (Mao et al., 2005) software to test the statistical enrichment of differential expression genes in KEGG pathways.

Weighted gene co‐expression network analysis (WGCNA)

There were 54,175 genes annotated in all the RNA‐seq libraries. After filtering out non‐expressed genes, low‐expressed genes (FPKM < 2) and genes that possessed large expression variation across replicates, 28,223 genes were selected for WGCNA analysis. An unsigned gene co‐expression network was built using the WGCNA R‐package (Langfelder & Horvath, 2008). The topological overlap measure (TOM) was calculated using the adjacency matrix between all genes with a power of 18. A gene dendrogram was created by using the dissimilarity TOM. Branches of the dendrogram (Figure S6) indicated modules (clustered highly co‐expressed genes) through using the dynamic tree cut algorithm (Langfelder et al., 2008). Each module has an assigned colour. Similar modules were merged by calculating module eigengenes. Using the module eigengenes, the module–trait relationships (MTRs) were plotted by calculating the Pearson's correlations between the module eigengene and the samples (traits). A correlation of 0.75 with at least one of the selected traits was the criteria used to generate modules. Clustered genes in modules of interest were then used for further analysis. Genes annotation and GO term analysis were performed in Soybase (Grant et al., 2010) and AgriGO (Du et al., 2010), respectively.

ROS production and contents measurement

Fourteen‐day‐old NIL‐T and NIL‐S soybean seedlings were used to measure the generation of ROS in roots. Roots were labelled with 30 μM dichlorofluorescein diacetate (DCFDA; Sigma‐Aldrich), by incubating washed soybean roots for 10 min at room temperature in the dark with 30 μM DCFDA. ROS sensitive fluorescent dyes were imaged using a Nikon SMZ25 stereo microscope (Nikon). The excitation and emission wavelengths were 488 nm and 500/535 nm, respectively. Fluorescence intensity was measured by ImageJ/FIJI (Schindelin et al., 2012), signal intensity values were normalised by setting the mean value for NIL‐T control to 1. Superoxide (O2 −) anions were also measured followed by Li et al. (2019), using 0.2% NBT (nitroblue tetrazolium) in 50 mM sodium phosphate buffer (pH 7.5). Soybean root samples were submerged in NBT staining solution for 10 min, which was the time it took samples to develop the dark blue colour.

ROS contents in soybean root samples were measured using Amplex® UltraRed reagent (Invitrogen, USA). Sodium phosphate buffer (0.5 ml 50 mM pH 7.4) was added to 0.1 g soybean powder, after mixing and solubilisation, samples were set on ice for 5 min. Then tubes were centrifuged at 12,000g for 20 min. 500 μl supernatant was transferred into new tubes. An equal volume of 2:1 (vol/vol) chloroform:methanol was added and centrifuged at 12,000g for 5 min. 50 μl aqueous phase was taken from each sample and added into each well of 96‐well microplate. 50 μl working solution (freshly made) was added into each well. Plates were incubated at 25°C for 30 min (protected from light). Fluorescence was measured at 540/590 nm.

Scavenging activity of the superoxide anion (O −) assay

The scavenging activity assay was adapted from Pi et al. (2016) with slight modifications. Less root homogenate (0.1 g) was used for measurement. Antioxidant enzymes were extracted with 2 ml of 0.05 M phosphate buffer (pH 5.5) from 0.1 g root homogenate. The extract was then centrifuged at 12,000g (4°C) for 10 min. Supernatant (40 μl) was added into 160 μl reaction buffer, which contained 80 μl phosphate buffer, 40 μl 0.05 M guaiacol (Sigma, USA), 40 μl 2% hydrogen peroxide (H2O2). The increased absorbance at 470 nm in 96‐well plates due to the enzyme‐dependent guaiacol oxidation was recorded every 30 s until 4 min of reaction.

Real‐time quantitative reverse‐transcription polymerase chain reaction (RT‐qPCR) validation for RNA‐seq results

Total RNA was extracted from soybean root tissues using the TRIZOL reagent (Ambion, http://www.ambion.com). To remove residual DNA, the extracted RNA was treated with RNase‐free DNase I (New England Biolabs, https://www.neb.com) for 30 min at 37°C. For gene expression, first‐strand cDNA synthesis was performed with a PrimeScript RT Reagent Kit (TaKaRa, http://www.takara.co.jp/english). Real‐time PCR was performed using SYBR Premix Ex Taq II (TliRNaseH Plus) (TaKaRa). The level of GmSALT3 transcript was normalised using the control gene GmUKN1 (Hu et al., 2009).

RESULTS

RNA‐sequencing preparation and profiles

GmSALT3 has previously been shown to be expressed in roots, and grafting experiments have shown that presence of GmSALT3 expression in roots is sufficient to confer shoot Na+ and Cl− exclusion (Guan et al., 2014; Qu et al., 2021). Total root RNA was extracted from NIL‐T and NIL‐S, and transcriptomes obtained using RNA‐seq in a time course from plants grown under control and saline conditions. To investigate short‐ and long‐term responses, root samples were harvested at three time points, 0 h (control), 6 h, and 3 days with 200 mM salt‐treatment and controls (Figure 1A and B) and in biological triplicates. Thirty RNA libraries were generated for paired‐end reads using an Illumina HiSeq 1500 sequencer. In total, 1.6 billion paired 100 bp raw reads were generated and mapped to the soybean genome sequence Gmax_275 Wm82.a2.v1 (Glyma 2.0) using TopHat2 (Kim et al., 2013). The average mapping percentage was 81.25%. After trimming of low quality (Q < 30), adapter fragments and barcode sequences, a total of 804 million clean mapped reads were identified. Combined with a quality control test using fastQC (Andrews, 2010), the quality of RNA‐seq libraries construction and sequence alignment was deemed sufficient for further analysis. A summary of mapped reads and quality of sequencing is shown in Table S1. Transcripts corresponding to GmSALT3/Gmsalt3 in the NIL‐T and NIL‐S in all samples were confirmed. NIL‐T transcriptomes contained reads across the entire full‐length transcript for GmSALT3 were obtained, while in NIL‐S reads only mapped to the predicted truncated version Gmsalt3 (Figure 1C). The phenotype of plants grown at the same time but not harvested and retained in the salt treatment after 10 days are shown in Figure 1D, to show the differential downstream impact of salt on NIL‐T and NIL‐S. A flowchart of data analysis is shown in Figure 1E.

FIGURE 1

Treatment, sampling strategy and data analysis. (A) Soybean seedlings were treated with 200 mM NaCl (salt treatment) or water (control) at 10 days after sowing (DAS). Arrows mean the sampling time points (0 h, 6 h, 3 days), triangle indicated the time point when water or NaCl solution was applied. (B) Pictures of soybean that had been treated with NaCl solution for 3 days. (C) Mapped reads to GmSALT3 (bottom) and Gmsalt3 (above). (D) The phenotype of NIL‐T and NIL‐S after 11 days salt treatment. (E) Flowchart of data analysis. NIL‐T, NIL with GmSALT3 allele; NIL‐S, NIL with Gmsalt3 allele; DEGs, differentially expressed genes; GO, gene ontology

Overview of DEGs between NIL‐T and NIL‐S under control conditions

Previous work had shown that GmSALT3 is expressed under both control and saline conditions, suggesting that there might already be a difference in other gene expression between NIL‐T and NIL‐S under control conditions. We therefore first compared the transcriptomes of NIL‐T and NIL‐S (indicated as “T” and “S”, respectively) under control conditions, before analysing changes due to salinity treatment. A PCA plot (for the first two principal components) of 10 grouped samples shows a good separation between all comparisons, confirming replicates clustered together consistently. Those comparisons are: Control 0 h T versus Control 0 h S (grey), Control 3 days T versus Control 3 days S (purple), and Control 6 h T versus Control 6 h S (brown) (Figure 2A). Transcript abundance was compared in FPKMs calculated using the Cufflinks functions cuffquant and cuffnorm (Trapnell et al., 2012). DEGs were classified as genes with a Log2FC (fold change) ≥ 1 and a FDR (false discover rate) < 0.01 when comparing between two conditions.

FIGURE 2

Overview of DEGs between NIL‐T salt‐tolerant and NIL‐S salt‐sensitive soybean samples at the 0 h, 6 h, and 3 days timepoints under non‐saline conditions. (A) PCA plot of 10 groups (g1 − g10). Different coloured dots represent replicates in each group, different coloured circles represent comparisons between groups (g1 vs. g2, g3 vs. g4, g5 vs. g6). (B) Venn diagram of up‐ and downregulated DEGs in NIL‐S soybeans at 0 h, 6 h, and 3 days under control conditions. (C) Three up‐regulated DEGs in NIL‐S soybeans at all 0 h, 6 h, and 3 days under control conditions, and their expression profile is shown in (D). The cutoffs for DEGs are Log2FC (fold change) ≥ 1, FDR (false discover rate) < 0.01

There are 5 upregulated differentially expressed genes (DEGs) at 0 h, 9 upregulated DEGs at 6 h, and 6 upregulated DEGs at 3 days in NIL‐S (compared to NIL‐T), and a Venn diagram demonstrates 3 of these genes are common across three time points. There were 1, 5, and 5 downregulated DEGs at 0 h, 6 h, and 3 days, respectively, but the identity of the downregulated DEGs were different between each of the time point (Figure 2B). The three consistently upregulated DEGs in NIL‐S at all time points under control conditions were Glyma.07G196800 (Linoleate 13S‐lipoxygenase 3–1), Glyma.10G143600 (uncharacterized protein), and Glyma.20G105500 (3‐hydroxybenzoate 6‐hydroxylase 1‐like) (Figure 2C). Figure 2D shows the expression profile (in FPKM, fragments per kilobase of transcript per million) of these three genes at all time points of NIL‐T and NIL‐S under control conditions.

Overview of DEGs in salt treated NIL‐T and NIL‐S when compared to their respective controls

Gene expression differences in response to salt of each genotype was then examined. By comparing control to salt treated NIL‐T, we identified 1816 DEGs that were differently expressed in NIL‐T roots in response to salt (1263 upregulated and 553 downregulated, 6 h T). We then compared expression of control and salt treated NIL‐S, and after 6 h salt treatment NIL‐S roots (when compared to control NIL‐S) showed a greater initial response in term of DEGs than we found in NIL‐T, with a total of 3054 DEGs (1911 upregulated and 1143 downregulated, 6 h S) (Figure 3B). After 3 days of salt treatment, the increased number of DEGs in NIL‐S was not observed anymore. On the contrary, the NIL‐T now showed a slightly higher number of DEGs compared to control conditions. NIL‐T had 2844 DEGs (1333 upregulated and 1511 downregulated, 3 days T) while there were 2573 DEGs (1318 upregulated and 1255 downregulated, 3 days S) in NIL‐S roots (Figure 3B).

FIGURE 3

Overview of DEGs (differentially expressed genes) between salt‐stressed and control samples in NIL‐T and NIL‐S soybeans at 6 h and 3 days. (A) PCA plot of 10 groups (g1–g10). Different coloured dots represent replicates in each group, different coloured circles represent comparisons between groups (g5 vs. g9, g6 vs. g10, g3 vs. g7, g4 vs. g8). (B) Numbers of DEGs (up and downregulated) between groups. Clustering of DEGs in log2FPKMs (fragments per Kilobase of transcript per million) of groups 3–10 (C–F). The false colour scale from green through to red indicates increasing log2FPKM. (G) Venn diagram of upregulated DEGs in NIL‐T and NIL‐S soybeans at 6 h and 3 days under NaCl treatment. (H) Venn diagram of downregulated DEGs in NIL‐T and NIL‐S soybeans at 6 h and 3 days under NaCl treatment. The cutoffs for DEGs are Log2FC (fold change) ≥ 1, FDR (false discover rate) < 0.01

Hierarchical clustering of the DEGs of the three replicates in different grouped samples is shown in Figure 3C–F. These DEGs represent the genes in each soybean NIL that are up‐ or downregulated in response to salt. We used the hierarchical clustering to identify DEGs, which either similarly change in expression in the two genotypes in response to salt, or show a different change in expression between the two genotypes. Analyses of the similarities and differences between the tolerant and sensitive lines indicated that 341 and 989 DEGs are uniquely upregulated after 6 h salt treatment in NIL‐T or NIL‐S, respectively; and 608 and 593 in the 3 days salt treatment (Figure 3G). We found similar values for uniquely downregulated DEGs (Figure 3H).

GO analysis of DEGs between control and salt‐treated samples of NIL‐T and NIL‐S

To gain insight into the putative functions of the uniquely up‐ and downregulated DEGs at 6 h salt treatment, the identified genes were subjected to GO (Gene Ontology) enrichment analysis, the representative GO terms are shown in Table 1, and detailed GO terms are shown in Figures S1–S4.

TABLE 1

GO term analysis of uniquely up and downregulated genes under salt treatment in soybean roots at 6 h

GO term	Ontology	Description	Number in input list	Number in BG/ref	p value	FDR
(a)
GO:0009607	BP	Response to biotic stimulus	7	64	9.3e−08	2.4e−05
GO:0006952	BP	Defence response	7	118	4.4e−06	0.00056
GO:0055114	BP	Oxidation reduction	28	2408	6.1e−05	0.0051
GO:0005507	MF	Copper ion binding	9	219	3.4e−06	0.00071
GO:0016491	MF	Oxidoreductase activity	32	2744	1.5e−05	0.0016
GO:0004866	MF	Endopeptidase inhibitor activity	5	92	0.00016	0.0083
(b)
GO:0006355	BP	Regulation of transcription, DNA‐dependent	69	1874	2e−05	0.0016
(c)
GO:0016021	CC	Integral to membrane	11	1419	0.00065	0.0058
(d)
GO:0020037	MF	Heme binding	23	728	8e−05	0.012
GO:0009055	MF	Electron carrier activity	21	702	0.00033	0.018
GO:0004190	MF	Aspartic‐type endopeptidase activity	9	169	0.00042	0.018
GO:0015035	MF	Protein disulphide oxidoreductase activity	7	96	0.00031	0.018

Note: (a) Upregulated genes in NIL‐T soybean roots. (b) Upregulated genes in NIL‐S soybean roots. (c) Downregulated genes in NIL‐T soybean roots. (d) Downregulated genes in NIL‐S soybean roots.

Abbreviations: CC, cellular component; BP, biological process; MF, molecular function.

The enriched GO terms for the 6 h T (NIL‐T) upregulated DEGs were all related to stress response such as “oxidation reduction”, “defence response” and “response to biotic stimulus” (Table 1a). Downregulated GO terms in NIL‐T were less specific and mainly associated with “integral to membrane” (Table 1c). For 6 h S (NIL‐S), a large number of upregulated genes were associated with the very broad GO term “regulation of transcription, DNA‐dependent (GO:0006355)”, and no further downstream classification was considerably enriched (Table 1b). Downregulated GO terms in NIL‐S included “heme binding”, “electron carrier activity,” “aspartic‐type endopeptidase activity,” and “protein disulphide oxidoreductase activity” (Table 1d), relating to protein types instead of specific biological process GO terms.

After 3 days of salt treatment, GO analysis of uniquely up‐ and downregulated DEGs showed that only “oxidation reduction” and “oxidoreductase activity” were upregulated in NIL‐T (Table 2a). “Oxidation reduction” and “oxidoreductase activity” were consistently upregulated in NIL‐T at both time points. Table 3 shows all the 53 upregulated genes in “oxidation reduction” and “oxidoreductase activity” GO terms after 3 days treatment in NIL‐T roots. A group of Cytochrome P450 enzymes‐encoding genes were significantly more highly expressed in NIL‐T, especially Glyma.13G173500, which had a FPKM of 279 compared to 71 under control conditions (Table 3). Other genes encoding oxidoreductase enzymes were also included such as peroxidases and dehydrogenases (Table 3). All the gene annotations in each GO are shown in Table S1.

TABLE 2

GO term analysis of uniquely up and downregulated genes under salt treatment in soybean roots at 3 days

GO term	Ontology	Description	Number in input list	Number in BG/ref	p value	FDR
(a)
GO:0055114	BP	Oxidation reduction	53	2408	8.9e−06	0.0037
GO:0016491	MF	Oxidoreductase activity	57	2744	2e−05	0.0073
(b)
GO:0006355	BP	Regulation of transcription, DNA‐dependent	51	1874	9.8e−08	3.6e−06
GO:0006468	BP	Protein amino acid phosphorylation	45	2356	0.002	0.023
GO:0048544	BP	Recognition of pollen	14	135	3.8e−09	3e−07
GO:0015743	BP	Malate transport	5	34	0.0001	0.0013
(c)
GO:0007018	BP	Microtubule‐based movement	17	155	9.8e−09	5.2e−06
GO:0003777	MF	Microtubule motor activity	17	155	9.8e−09	2.4e−06
GO:0043531	MF	ADP binding	27	467	2.1e−07	2e−05
GO:0004565	MF	Beta‐galactosidase activity	6	37	9.5e−05	0.0031
GO:0017171	MF	Serine hydrolase activity	17	330	0.00014	0.0041
GO:0009341	CC	Beta‐galactosidase complex	6	37	9.5e−05	0.009
(d)
GO:0055085	BP	Transmembrane transport	29	1167	1.9e−05	0.0075
GO:0022857	MF	Transmembrane transporter activity	24	1020	0.00022	0.036

Note: (a) Upregulated genes in NIL‐T soybean roots. (b) Upregulated genes in NIL‐S soybean roots. (c) Downregulated genes in NIL‐T soybean roots. (d) Downregulated genes in NIL‐S soybean roots.

Abbreviations: CC, cellular component; BP, biological process; MF, molecular function.

TABLE 3

Upregulated genes in oxidation reduction and oxidoreductase activity GO terms of 3 days salt‐treated NIL‐T roots

Note: Annotations are achieved from https://www.soybase.org/genomeannotation/.

KEGG analysis of DEGs between control and salt‐treated samples

To understand what pathways were differently altered in NIL‐T and NIL‐S under salt stress, KEGG (Kyoto Encyclopaedia of Genes and Genomes) analysis was performed. This analysis can reveal enrichments of biological pathways, similar to GO term analysis but KEGG takes fold changes of DEGs into account, while the GO term analysis does not.

In our analysis, 1816 and 3054 DEGs in NIL‐T and NIL‐S, respectively, were enriched in 58 pathways (corrected p‐value <0.05) after 6 h salt‐treatment (Figure 4A). Most of the 58 pathways in both NIL‐T and NIL‐S were within the metabolism category with “phenylpropanoid biosynthesis,” “phenylalanine metabolism,” and “starch and sucrose metabolism” as the top three pathways (Figure 4A). Overall, NIL‐S has more DEGs compared to NIL‐T after salt‐treatment, across all significantly changed pathways.

FIGURE 4

Significantly enriched KEGG (Kyoto encyclopaedia of genes and genomes) of all DEGs between salt‐stressed and control samples in NIL‐T and NIL‐S soybeans. (A) 6 h KEGG enrichment. (B) 3 days KEGG enrichment. KEGG enrichment p‐value cutoff ≤0.05

After 3 days of salt stress (3 days), 2844 and 2573 DEGs in NIL‐T and NIL‐S, respectively, were mapped to 57 pathways (corrected p‐value <0.05) (Figure 4B). The metabolism category was also the most enriched after 3 days with similar pathways involved compared to 6 h. After 3 days salt‐treatment, most of the significantly changed pathways included more DEGs in NIL‐T compared to NIL‐S, in particular in the pathway “Protein processing in endoplasmic reticulum” (16–11 DEGs), the cell organelle where GmSALT3 is localised (Guan et al., 2014). DEG enrichment in NIL‐T was also higher in the pathways “ubiquitin mediated proteolysis” (10–3 DEGs), “phenylpropanoid biosynthesis” (71–60 DEGs), “phenylalanine metabolism” (52–42 DEGs), “starch and sucrose metabolism” (43–35 DEGs), and “plant‐pathogen interaction” (33–20 DEGs). It is important to note that all plants (NIL‐T and NIL‐S) were grown together and no plants showed signs of infection, suggesting that the genes categorised into this pathway also have a function in abiotic stress response.

In the “plant‐pathogen interaction” pathway, the putative CNGC (cyclic nucleotide‐gated ion channel) 15‐like gene, Glyma.13G141000, was significantly downregulated in 3 days T (Figure S5a and S5b), and the CNGC 20‐like gene, Glyma.09G168700, was upregulated in 3 days S. Several CaM (Calmodulin) and CML (Calmodulin‐like) genes were downregulated in 3 days T. CDPK (calcium‐dependent protein kinase 3, Glyma.08G019700) and Rboh (respiratory burst oxidase homologue protein F‐like isoform 1, Glyma.01G222700) were upregulated in 3 days S. However, at 6 h, another CNGC 20‐like gene (Glyma.16G218300), homologue to the Arabidopsis putative Ca2+ channel that is proposed to regulate cytosolic Ca2+ activity and expression of Rhob and CaM/CML (Demidchik et al., 2018), was significantly upregulated in NIL‐T (Figure S5d).

WGCNA identifies genes that have strongly correlated expression profiles, and those genes work cooperatively in related pathways to contribute to corresponding phenotypes. It is based on the concept that genes with strongly correlated expression profiles are clustered because those genes work associatively in related pathways, contributing to the resulting phenotype; in our case, establishing salinity tolerance after salt treatment. We applied the WGCNA to the normalised FPKM data (FPKM > 2) from all 30 RNA‐seq libraries.

A co‐expression network was constructed, and highly co‐expressed genes were assigned to different colour‐coded modules (Figure S6). Clusters were summarised by correlating the module's eigengene (or an intromodular hub gene) to the sample traits (Langfelder & Horvath, 2008). Module–trait relationships (MTRs) were computed based on Pearson's correlation between module and phenotypes, and then be used to screen modules for downstream analysis. Thirty clusters of highly co‐expressed genes (modules) were detected and are shown in differently coloured module (Figure 5). To investigate what genes show strongly correlated co‐expression in NIL‐T under saline conditions (6 h and 3 days), two modules of interest were selected for functional annotation based on their correlation value and corresponding p value. One module contains 79 genes that show strongly correlated co‐expression in NIL‐T after 6 h salt treatment; this model is depicted here in a shade of pink (MEdeeppink). The second selected module, colour‐coded here in a shade of white (MEantiquewhite2), contains 406 genes which expression profiles are highly correlated in NIL‐T after 3 days salt treatment.

FIGURE 5

Module–trait relationships (MTRs) and corresponding p‐values (within brackets) between the detected modules on the y‐axis and groups (trait) on the x‐axis. The MTRs are coloured based on their correlation: from red (strong positive correlation) to blue (strong negative correlation). Each group contains three biological replicates. Pentagram represents modules of interest

GO analysis of 79 genes in module MEdeeppink is shown in Table 4a; GO analysis of 406 genes in module MEantiquewhite2 is shown in Table 4b. Within the highly co‐expressed 79 genes in NIL‐T after 6 h salt treatment (MEdeeppink module), Gene Ontology (GO) analysis indicates 12 of them are integral to plasma membrane (GO:0016021; Table 4a), and they encode a group of Casparian strip membrane proteins (Table S1). Highly co‐expressed 406 genes in NIL‐T after 3 days salt treatment (MEantiquewhite2 module) are enriched in GO:0016021 (integral to membrane; 26 genes), which are related to intracellular vesicle trafficking and ion transport (Table 4b). Within those GOs, there are genes relevant to vesicle transport, such as genes encoding for coatomer subunits and transmembrane protein transporters; three potassium transporters including Glyma.05g213100 (TRH1; tiny root hair 1), Glyma.08g288500 (KUP6; potassium ion uptake transporter 6), and Glyma.13g170700 (Table S1).

TABLE 4

Significant gene ontology (GO) terms in MEdeeppink and MEantiquewhite2 modules

GO term	Ontology	Description	Number in input list	Number in BG/ref	p value	FDR
(a)
GO:0005886	C	Plasma membrane	12	81	5.00e−20	1.70e−18
GO:0031224	C	Intrinsic to membrane	12	1453	4.40e−06	4.90e−05
GO:0016021	C	Integral to membrane	12	1419	3.50e−06	4.90e−05
GO:0044425	C	Membrane part	13	1801	7.00e−06	5.80e−05
GO:0016020	C	Membrane	15	3518	0.00053	0.0035
(b)
GO:0009207	P	Purine ribonucleoside triphosphate catabolic process	6	40	0.0000092	0.00044
GO:0046907	P	Intracellular transport	14	297	0.0000073	0.00044
GO:0016567	P	Protein ubiquitination	9	130	0.000019	0.00076
GO:0046034	P	ATP metabolic process	5	46	0.00021	0.0046
GO:0046578	P	Regulation of Ras protein signal transduction	6	91	0.00061	0.011
GO:0004674	F	Protein serine/threonine kinase activity	14	39	3.1e−16	1.2E‐13
GO:0046872	F	Metal ion binding	65	3600	0.000024	0.0021
GO:0042623	F	ATPase activity, coupled	12	343	0.00047	0.021
GO:0030120	C	Vesicle coat	5	41	0.00013	0.0013
GO:0031981	C	Nuclear lumen	5	105	0.0066	0.039
GO:0016021	C	Integral to membrane	26	1419	0.0068	0.039

Note: (a) Enriched GO terms in MEdeeppink module (NIL‐T after 6 h salt treatment); (b) enriched GO terms in MEantiquewhite2 module (NIL‐T after 3 days salt treatment).

Abbreviations: C, cellular component; P, biological process; F, molecular function.

RT‐qPCR validation for RNA‐seq results

In order to confirm the RNA‐seq results, 10 DEGs were selected for RT‐qPCR validation based on their RPKM transcript abundance. Importantly, this occurred on independently grown plants not used for RNAseq analysis but grown in identical conditions. RT‐qPCR indicated that relative expression values (relative to housekeeping gene GmUKN1) of the selected DEGs were significantly correlated with their FPKM values (Figure S7).

ROS production and scavenger enzymes activity

Due to the enrichment of the “Oxidation reduction” gene ontology category in NIL‐T plants under salt (Table 3), we decided to test if this translated into a difference in ROS generation or detoxification in NIL‐T compared to NIL‐S, again on independently grown plants. The ROS activity was measured in the roots of NIL‐T and NIL‐S after a 3‐day with or without salt treatment, using a fluorescent dye (Figure 6A and B). We could detect a lower ROS activity in NIL‐T compared to NIL‐S roots, already in the control conditions, consistent with the expression of GmSALT3 under control conditions. This difference strongly increased when comparing plants from saline conditions; ROS were strongly produced in NIL‐S under salt treatment. The quantified fluorescence intensity of staining confirmed that under both control and salt treatment, NIL‐T could maintain a significantly lower ROS content in the roots (Figure 6C).

FIGURE 6

ROS production and enzyme activity of the superoxide anion (O2 −) scavenger assay in soybean roots. NIL‐T and NIL‐S were labelled with 30 μM DCFDA for 10 min under control conditions (A) and with 200 mM NaCl salt treatment for 3 days (B). Bar = 1000 μm. Green fluorescence indicates the presence of ROS, and the mean signal intensity (C) was measured using ImageJ, asterisk indicates a significant difference between NIL‐T and NIL‐S at **p < 0.001, ****p < 0.0001, n = 5–6. NBT staining indicates the superoxide content of the root tips of NIL‐T and NIL‐S under control conditions (D) and with 200 mM NaCl salt treatment for 3 days (E), the staining intensity reflects the concentration of superoxide. Bar = 1 mm. (F) Scavenging activity of the superoxide anion (O2 −) of NIL‐T (black) and NIL‐S (pink) with or without salt treatment. Asterisk indicates a significant difference between NIL‐T and NIL‐S at *p < 0.05, ****p < 0.0001, n = 3

ROS are composed of many different molecules; we initially measured the concentration of the ROS H2O2. We could detect a higher H2O2 concentration in both lines under salt treatment compared to control; but no concentration differences between NIL‐T and NIL‐S in either control or salt treatment could be measured (Figure S9), suggesting that the increase in ROS is attributed to a different molecule. Therefore, we used Nitroblue tetrazolium (NBT) to measure the production of ROS superoxide, and the results (Figure 6D and E) showed that the staining intensity was much weaker in the NIL‐T than NIL‐S after 3 days salt treatment, which indicates NIL‐T accumulated lower concentrations of superoxide than NIL‐S in the roots. The antioxidant properties of the roots of NIL‐T and NIL‐S with or without salt treatment was then analysed using guaiacol. The antioxidant enzyme activity was assayed by measuring the absorbance at 470 nm due to the enzyme dependent oxidation of guaiacol by H2O2 (Pi et al., 2016). The scavenging enzyme activity of the superoxide anion was significantly higher in NIL‐T compared to NIL‐S under control and saline conditions (Figure 6F).

DISCUSSION

Soybean cultivars harbouring GmSALT3 are better able to modulate Na+, K+, and Cl− content of shoots under salinity (Do et al., 2016; Guan et al., 2014; Liu et al., 2016; Qi et al., 2014; Qu et al., 2021), and they have a yield advantage under saline conditions compared to soybean varieties without a functional GmSALT3 (Do et al., 2016; Liu et al., 2016). The GmSALT3 ion transport protein is ER‐localised (Guan et al., 2014), suggesting a cellular function of GmSALT3 that is not directly connected to ion uptake or exclusion. We performed detailed and extended RNAseq analysis of NIL‐T and NIL‐S soybeans to decipher the cellular pathways and processes that connect GmSALT3 to soybean salinity tolerance. We could identify several genes and biological pathways involved in responses to salt stress in GmSALT3‐containing NIL‐T plants, which gives an insight into how GmSALT3 expression translates into a phenotype thorough modulating cell processes.

Illumina sequencing results showed that all the 30 constructed RNA‐seq libraries were of sufficient quality for further analysis. On average the RNA‐seq libraries had 26.8 million clean mapped reads, and reads were mapped to 53,625 annotated soybean genes with 353 additional transcripts not previously annotated—probably due to the germplasm being genetically divergent from the reference. Minimal DEGs could be detected between NIL‐T and NIL‐S at 0 h, 6 h, and 3 days under control conditions (Figure 2), this indicates that NIL‐T and NIL‐S roots had similar gene expression patterns under non‐stressed conditions. In the three consistently upregulated DEGs in NIL‐S at all time points under control conditions, there was a gene encoding for a LOX (Linoleate 13S‐lipoxygenase 3‐1) (Figure 2C). Lipoxygenases (LOXs) are involved in defence responses against biotic stresses, such as microbial pathogens and insect pests, in tomato (Hu et al., 2015), potato (Royo et al., 1999), maize (Gao et al., 2009), and rice (Zhou et al., 2009). The full length GmSALT3 during soybean natural selection and domestication was lost several times independently (Guan et al., 2014), suggesting that plants with a non‐functional GmSALT3 might have a benefit when grown under non‐saline conditions. This benefit might be related to higher LOXs expression in soybean cultivars harbouring the sensitive alleles, which could be beneficial under certain biotic stresses. In potato, depletion of one specific LOX (LOX‐H3) negatively affected their response to stress (Royo et al., 1999), but here in soybean, plants with the full length GmSALT3 had lower LOX expression than with the non‐functional GmSALT3. Loss of GmSALT3 function appears to occur in non‐saline regions so this might confer an advantage to the plant under some conditions (Guan et al., 2014). It is possible that greater LOX expression in soybean confers an advantage in certain conditions not studied here, such as biotic stress, which could be an area of further study.

Comparing DEGs, the GO terms “oxidation reduction (GO:0055114)” and “oxidoreductase activity (GO:0016491)” were consistently upregulated in NIL‐T at 6 h and 3 days of salt treatment, with a greater number of genes in this term upregulated after 3 days, which indicates that NIL‐T plants may have a greater capacity to detoxify ROS (reactive oxygen species) than NIL‐S plants. We could experimentally confirm this and show that more ROS are present in NIL‐S roots compared to NIL‐T (Figure 6A–E) and NIL‐T possessed a significantly higher scavenging enzyme activity of the superoxide anion (O2 −) (Figure 6F). Phenylalanine is the precursor for flavonoids, and flavonoids are utilised to scavenge and protect against ROS (Dastmalchi et al., 2016). The increase in gene expression related to phenylalanine synthesis seen in NIL‐T (Figure 4A and B) is consistent with the observation that GmSALT3 is connected to an increase in general ROS detoxification capacity. ROS detoxification and signalling is frequently linked to improvements in stress tolerance in plants (Dong et al., 2013; Munns & Gilliham, 2015; Munns & Tester, 2008; Roy et al., 2014). For instance, the exogenous application of scavengers such as the non‐enzymatic antioxidant ascorbic acid has previously improved plant performance under stress (Akram et al., 2017), but such treatments are impractical in the field and are likely to have additional impact on the plant and microbial communities as a carbon source. Here, the expression of GmSALT3 is correlated with the ability to detoxify ROS, and therefore represents an upstream genetic mechanism likely to enhance this capacity over germplasm not expressing GmSALT3.

In all the listed upregulated genes in NIL‐T (Table 3), a group of Cytochrome P450 enzymes‐coding genes were significantly more highly expressed in NIL‐T in response to salt‐treatment, especially Glyma.13G173500, which has the highest expression level and a fold change of 3.93. The cytochrome P450 (CYP) family is a large and essential protein superfamily in plants, and CYPs catalyse monooxygenation/hydroxylation reactions in primary and secondary metabolism pathways (Mizutani & Ohta, 2010). NaCl can induce conformational change of cytochrome P450 in animals (Oyekan et al., 1999; Yun et al., 1996), but this has not been investigated in plants. In soybean, there are 322 identified CYPs, with only few of them having been functionally characterised (Guttikonda et al., 2010). GmCYP82A3 (Glyma.13g068800) (Yan et al., 2016) and GmCYP51G1 (Glyma.07g110900) (Pi et al., 2016) were shown to be involved in plant tolerance stresses such as salinity and drought. The two P450 containing isoflavone synthases (IFS) GmIFS1 (Glyma.07G202300) and GmIFS2 (GmCYP93C1 = Glyma.13G173500) are both anchored in the ER, the same organelle where GmSALT3 is localised. GmIFS2 is a DEG identified in our study. Expression of GmIFS1 was previously shown to be induced by salt stress, and has been linked to improving salinity tolerance through the accumulation of isoflavone content (Jia et al., 2017). We found its homologue, GmIFS2 (Glyma.13G173500), had a fold change of 3.93 in NIL‐T in response to salt stress (Table 3). Combined, this suggests that the function of GmSALT3 could be important for the synthesis of isoflavones, and this may contribute to salt tolerance.

ROS production has been linked to Ca2+ influx into cells and specifically the activity of CNGCs (Ma et al., 2009). Ca2+ and ROS act as signalling molecules when plant roots are suddenly exposed to salt (Byrt et al., 2018; Stephan & Schroeder, 2014). Gene expression related to proposed Ca2+ signalling elements are different between NIL‐T and NIL‐S. As a result, our analysis suggested that the “plant–pathogen interaction” pathway could be highly activated in NIL‐S, which may lead to higher ROS production and hypersensitive responses (Figure S5a and c). Under 6 h salt treatment, several CNGCs were also significantly upregulated in NIL‐S, which may relate to activation of different signalling cascades (Figure S5d). Several CaM (calmodulin) and CML (Calmodulin‐like) genes are also downregulated in 3 days T (Figure S5b), downregulation of CaM and CML may has been proposed to restrict gaseous reactive oxygen species (ROS) production in plants (Ma & Berkowitz, 2011).

The WGCNA analysis revealed that genes related to Casparian strip (CS) development were upregulated in NIL‐T within 6 h of salt treatment. CSs play an important role in regulating nutrient uptake and salt stress tolerance, as they form the diffusion barrier for ions into the stele (Chen et al., 2011; Roppolo et al., 2011). CS localised membrane proteins (CASPs) were shown to mediate CS formation in Arabidopsis, by guiding the local lignin deposition during cell wall modification (Roppolo et al., 2011). Upon immobilisation, CASPs recruit PER64 (Peroxidase 64), RBOHF (respiratory burst oxidase Homologue F), ESB1 (enhanced suberin 1), and dirigent proteins, to establish the lignin polymerisation system (Hosmani et al., 2013; Kamiya et al., 2015; Lee et al., 2013). Our analysis identified a module of genes upregulated under salt exclusively in the NIL‐T (MKdeeppink module). This module contains eight dirigent genes, peroxidase 64 (Glyma.14g221400) as well as CASPs. The combination of genes in this model suggests that NIL‐T plants might more efficiently induce an earlier formation of the endodermal diffusion barrier as well as the formation of an exodermal diffusion barrier, in response to salt stress. We stained lignin and suberin in NIL‐T and NIL‐S roots with ClearSee and Auramine O, following the protocol developed by Ursache et al. (2018), in order to compare the CS development. Unfortunately, we could only observe autofluorescence in both NIL‐T and NIL‐S root differentiation zone that close to the root tips using confocal microscopy. So the link between GmSALT3 and CS formation needs further validation.

Lastly, CHX transporters (of which GmSALT3 is a member) were shown to be associated with the endocytic and exocytic pathways in dynamic endomembrane system (Chanroj et al., 2012; Padmanaban et al., 2007). Our results and the previously proposed vesicle trafficking function of CHXs (Chanroj et al., 2012), suggests that GmSALT3 might be important for establishing optimal conditions in the ER lumen, which could subsequently impact the secretion of compounds such as signalling molecules, ions, and soluble proteins.

To summarise, our findings suggest that the GmSALT3 containing NIL‐T has a higher salinity tolerance because NIL‐T roots are more effective in detoxifying ROS, potentially through ER‐localised flavonoid biosynthesis enzymes, which helps to prevent cellular damage. Additionally, our study has identified that Ca2+ signalling, vesicle trafficking and Casparian strip development might be modulated by GmSALT3, and are therefore areas for further study.

AUTHOR CONTRIBUTIONS

MG, YQ and RG designed experiments, with input from SW and LQ. Yue Qu and Rongxia Guan analysed the data. Lijuan Qiu and Matthew Gilliham directed the project. Oliver Berkowitz and James Whelan contributed to RNA‐sequencing. Lili Yu performed plant samples harvest for RNA‐sequencing. Rakesh David contributed to RNA‐sequencing analysis. Yue Qu, Matthew Gilliham and Rongxia Guan wrote the paper. Stefanie Wege and Lijuan Qiu assisted in editing the final version of the manuscript. All authors commented on the manuscript.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

Figure S1 GO term analysis of unique DEGs under salt treatment in NIL‐T roots after 6 h 200 mM NaCl treatment.

Figure S2 GO term analysis of unique DEGs under salt treatment in NIL‐S roots after 6 h 200 mM NaCl treatment.

Figure S3 GO term analysis of unique DEGs under salt treatment in NIL‐T roots after 3 days 200 mM NaCl treatment.

Figure S4 GO term analysis of unique DEGs under salt treatment in NIL‐S roots after 3 days 200 mM NaCl treatment.

Figure S5. DEGs between salt‐treated and control samples of NIL‐T and NIL‐S in plant‐pathogen interaction KEGG pathway at 6 h and 3 days.

Figure S6 Gene dendrogram showing the co‐expression modules defined by the WGCNA labelled by colours.

Figure S7 qPCR validation for RNA‐seq results of selected genes.

Figure S8 Glyma.13G173500 (GmIFS2) gene expression (FPKM) in all samples.

Figure S9 H2O2 concentration measurement in soybean roots.

Click here for additional data file.

Table S1 Reads mapping and quality of sequencing for 30 Glycine max root samples.

Table S2 gene annotations in each GO.

Click here for additional data file.