Two sorghum restorer lines, LRNK1 (salt-tolerant (ST)) and LR2381 (salt-sensitive, SS), were used as plant materials. While they had different sensitivities to salt, the crop maturity period, growth rate, and phenotypic traits of the two sorghum genotypes (Bicolor L. Moench) were similar. LRNK1 and LR2381 seeds were collected from the sorghum seed resource bank in the Liaoning Academy of Agricultural Sciences with ID: 011427 and ID: 010362, respectively. After disinfection (surface-sterilized with 1% sodium hypochlorite for 5 min followed by rinsing three times with distilled water), the seeds were planted in a hydroponic beaker filled with sorghum hydroponic culture (Kimura solution, Supplementary Table S1 ) under artificially simulated climatic conditions: day/night temperatures of 28°C/25°C, approximately 50% relative humidity (Rh), photoperiod of 14:10 h of light (1,500 μmol·m⁻²·s⁻¹)/dark in April 2019 in an artificial climate room in Liaoning Provincial Academy of Agricultural Sciences.

Sorghum plants at the six-leaf stage were used for the salt treatment experiment. For each variety, 144 plants were randomly selected and then divided into six hydroponic beakers with 24 plants in each. These six hydroponic beakers were then randomly divided into two groups: salt treatment (T) and control (C), with three hydroponic beakers in each group. For group T, both LRNK1 and LR2381 were treated with an additional 180 mM NaCl in the basic sorghum hydroponic culture, which was based on previous salt stress simulations at 60, 120, 180, and 240 mM NaCl. The two extremes of 180 mM NaCl showed significant phenotypic differences between the salt stress and control treatments. The salt-sensitive lines survived, while the salt-sensitive lines died under 240 mM NaCl simulations. In addition, the salt stress adversity simulated by 180 mM NaCl most closely matched the actual production environment. For group C, both LRNK1 and LR2381 were treated with basic sorghum hydroponic culture (no NaCl supplement). A total of four sets of samples were prepared: LRNK1-control (STC), LRNK1-salt treatment (STT), LR2381-control (SSC), and LR2381-salt treatment (SST). After 72 h of cultivation, six plant samples were randomly selected from each hydroponic beaker and the fifth leaf at tillering on each plant was sampled (18 leaves obtained in each group). These 18 leaves were pooled and divided into six parts in each group, three of which were used for RNA-seq analysis and three for metabolomic analysis. All sampled leaves were immediately placed in liquid nitrogen for 5 min, stored at −80°C, and tested within 3 to 4 weeks. The root tissue of each group was sampled to determine the physiological indicators.

After the hydroponic salt stress sampling was completed, the seedlings of the four groups of STC, STT, SSC, and SST were transplanted into prepared pots treated with the same salinity as the hydroponics for cultivation until maturity, and the sorghum growth phenotype parameter was recorded.

The plant height and leaf area (the fifth leaf at tillering) of seedlings at the six-leaf stage with and without salt treatment were measured according to a previous report (Iqbal et al., 2010). Fresh leaves were determined by an electronic balance to obtain the fresh mass (FM) and dried at 105°C for half an hour and at 70°C for 48 h to achieve constant weight. The dry mass (DM) was recorded. Turgid mass (TM) was determined after the leaves had been immersed in distilled water for 12 h. The FM, DM, and TM of root tissues were also determined. The relative water content (RWC) of a plant tissue is expressed by RWC (%) = [(FM − DM)/(TM − DM)] * 100%. For the leaf chlorophyll content determination, the ethanol-acetone method was used as described by Zhang (1992). The same experimental protocol was repeated three times.

Total RNA was extracted using an RNAprep pure Plant Kit (Tiangen, Beijing, China) according to the manufacturer’s instructions, and the RNA quality was evaluated by gel electrophoresis. Total RNA was reverse transcribed to cDNA by a QuantScript RT Kit (Tiangen, Beijing, China). An efficient VAHTS^® mRNA-seq v2 Library Prep Kit for Illumina (Vazyme Biotech, Nanjing, China) was then used to construct the sequence libraries. All 12 mRNA libraries (six groups with three biological replicates in each group) obtained were evaluated by an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and the qRT-PCR method. Finally, all mRNA-seq libraries were sequenced on an Illumina HiSeq 4000 sequencing platform with pair-end 2 × 150 bp mode to obtain sequencing data.

The sequencing data were filtered using FastQC (version 0.11.5) with default criteria. Low-quality reads and adaptor reads were removed from the raw data, and the clean data were assembled and mapped to the sorghum assembly genome (https://www.ncbi.nlm.nih.gov/genome/term=Sorghum+bicolor) using HISAT2 (version 2.1.0) (Kim et al., 2015). The value of FPKM (expected number of fragments per kb of exon per million of reads) of reads in each sample was calculated using Cufflinks (version 2.2.1) (Roberts et al., 2011). Principal component analysis (PCA) and Pearson’s correlation analysis were performed based on the FPKM of reads. The differentially expressed genes (DEGs) by different comparisons were identified using DESeq2 (version 1.24.0) (Anders, 2010), with the criteria of the corrected p-value (padj) <0.05 by the negative binomial distribution test and |log₂(fold change (FC))|>0. Genes with log₂FC >0 and log₂FC <0 were identified as up- and downregulated DEGs, respectively. Hierarchical clustering based on the expression profiles of DEGs was presented by pheatmap (version 1.0.10). DEGs enriched into modules correlated with the phenotypes were separately subjected to the enrichment analysis for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa et al., 2007). Significant GO biological processes (BP) and KEGG pathways were identified with the criterion of p < 0.05.

The chromatographic column (Accucore HILIC column; Thermo Fisher Scientific) was maintained at 40°C with an injection volume of 5 µl. The mobile phase (0.3 ml/min flow rate) consisted of 0.1% formic acid, 95% acetonitrile, and 10 mM ammonium acetate in solvent A, and 0.1% formic acid, 50% acetonitrile, and 10 mM ammonium acetate in solvent B. The gradient elution phase of solvent A was 98% for the first 1 min, decreased to 50% during 1–17 min, stabilized for 0.5 min (17.5 min), returned to 98% after 17.5–18 min, and stabilized at 98% for the last 2 min (18–20 min). A 5-min interval was set for the equilibrating instrument. All samples were analyzed using a Q Exactive (QE) HF-X mass spectrometer detector (Thermo Fisher Scientific) in positive and negative (polarity) ionization modes (spray voltage of 3.2 kV) with MS/MS data-dependent scans (mass range of 100–1500 m/z). The parameters of the electrospray ion source were set as 320°C drying gas temperature, 35 arb sheath gas flow rate, and 10 arb Aux gas flow rate. System stability and accuracy were validated using QC samples with an interval of five samples.

The data were recorded using the compound discoverer (CD) software, and the parameters, including retention times and mass charge ratio of compounds, were analyzed and corrected. The relative standard deviation (RSD) and signal-to-noise (S/N) ratio were analyzed. After feature screening, candidate features (RSD >30%, S/R ratio <3, signal intensity >10⁵) were retained and aligned against mzCloud (https://www.mzcloud.org/). QC samples were used for data normalization, and the relative expression of the metabolites was obtained and used for further analysis. The resulting three-dimensional data involving the peak number, sample name, and normalized peak area were fed into R package metaX for PCA and projections to latent structures-discriminant analysis (PLS-DA) (Wen et al., 2017). A Student’s t-test was used to compare peak areas of three metabolites, with the threshold of variable importance in the projection (VIP) value >1, fold change (FC) ≥2 (|log2FC|≥1), and p < 0.05. A hierarchical clustering analysis was then performed to present the expression patterns of the differential metabolites. Subsequently, the correlations between the differential metabolites were calculated by R software (v3.1.3). Differential metabolites with significant correlations (false discovery rate, FDR <0.05) were retained and used for enrichment analysis of the KEGG pathway, with significant criteria of padj < 0.05.

Differences in expression levels were investigated using a one-way analysis of variance (ANOVA). The means of the values were compared using Tukey’s multiple comparison tests. A p-value of <0.05 was considered significant. For conjoint analysis, the correlation analysis between the top 100 DEGs and the top 50 differential metabolites was performed, and the p-value was calculated. KEGG enrichment analysis was then performed on these significantly related DEGs and differential metabolites. A p-value of <0.05 was considered significant.

Physiological indicators showed that plant height, leaf area, and chlorophyll contents showed no significant difference between the two sorghum varieties under normal conditions ( Figure 1 ). However, after salt treatment, these three indicators in the salt-tolerant sorghum LRNK1 were significantly higher than in the salt-sensitive LR2381. In addition, RWC in LR2381 leaves and roots had significantly lower contents than that in LRNK1, which revealed the serious water loss in the salt-sensitive LR2381 due to salt stress. The morphological differences between salt-tolerant and salt-sensitive sorghum under salt stress are obvious in the seedling period, jointing period, and mature period ( Table 1 ; Figure 1 ).

A total of 458.10 Mb of raw reads from the 12 libraries were generated, with an average sequencing error rate of 0.02% and an average GC content of 54.74% ( Supplementary Table S2 ). The average data ratio of sequencing quality to Q20 and Q30 was 97.96% and 94.34%, respectively. After read mapping analysis, an average of 92.59% of reads were uniquely mapped to the sorghum genome sequence ( Supplementary Table S2 ). PCA results showed that the biological repeat samples in each group are clustered together ( Supplementary Figure S1A ). Sample-to-sample correlation analysis results showed a similar result with PCA, indicating good consistency and reproducibility between biological repetitions ( Supplementary Figure S1B ).

A total of 10,771 DEGs were identified by pairwise comparison of all the groups. There were 2,354, 4,162, 719, and 1,559 upregulated DEGs and 1,863, 3,797, 1,169, and 1,576 downregulated DEGs found during the comparison of STT vs. STC, SST vs. SSC, STC vs. SSC, and STT vs. SST, respectively ( Figure 2A ). The expression pattern of all DEGs showed a significant distinction between different groups ( Figure 2B ). Venn diagram analysis results showed that there were 3,326 common DEGs between STT vs. STC and SST vs. SSC ( Figure 2C ).

To reveal the role of DEGs, GO and KEGG enrichment analyses were performed. GO enrichment results revealed that 11 GO terms were significantly enriched in STT vs. STC, including carbohydrate catabolic process (GO:0016052) and transmembrane transport (GO:0055085). Most of the DEGs were upregulated in the carbohydrate catabolic process, including SbFBA (fructose-bisphosphate aldolase), SbPGK (phosphoglycerate kinase), and SbPFK (pyrophosphate-fructose 6-phosphate) family members ( Figure 3A ). In addition, some GO processes related to DNA replication were also related based on downregulated DEGs in STT vs. STC, including DNA metabolic process (GO:0006259), mismatch repair (GO:0006298), and mismatched DNA binding (GO:0030983). SbCAT (catalase), SbPUT (polyamine transporter), SbMSH (DNA mismatch repair protein), and SbBAT (amino-acid permease) were the key differentially expressed regulators in these GO terms ( Figure 3B ). For the salt-sensitive species LR2381, 16 GO terms were enriched, most of which were related to photosynthesis, including photosynthesis (GO:0015979), photosystem II (GO:0009523), and thylakoid (GO:0009579). DEGs like SbpsbP (photosystem II oxygen-evolving complex protein), Sbpsa (photosystem I assembly factor), and SbLHCs (light-harvesting complex) played key roles in these processes ( Figure 3C ). KEGG results showed that biosynthesis of amino acids (sbi01230), carbon fixation in photosynthetic organisms (sbi00710), glutathione metabolism (sbi00480), and circadian rhythm–plant (sbi04712) were significantly enriched based on the DEGs in STT vs. STC, and most DEGs involved were upregulated with salt treatment. Most of the genes in the C4-dicarboxylic acid cycle (belonging to sbi00710) were active. SbPIF (ATP-dependent DNA helicase), SbHY5 (transcription factor HY5), and SbFT (FT-interacting protein) in circadian rhythm–plants were upregulated. For DEGs in the glutathione metabolism, SbNADP⁺ was downregulated and SbGST (probable glutathione S-transferase) was upregulated. As for SST vs. SSC, the upregulated DEGs were enriched in carbon fixation in photosynthetic organisms (sbi00710), plant hormone signal transduction (sbi04075), MAPK signaling pathway–plant (sbi04016), glutathione metabolism (sbi00480), and carbon metabolism (sbi01200). The downregulated DEGs were enriched in porphyrin and chlorophyll metabolism (sbi00860), photosynthesis (sbi00195), carbon metabolism (sbi01200), starch and sucrose metabolism (sbi00500), and the TCA cycle (sbi00020). DEGs in plant hormone signal transduction showed that SbPP2C (probable protein phosphatase 2C), SbABF, SbIAA (Auxin), SbEIN3 (EIN3-binding F-box protein), and all genes related to salicylic acid (including SbJAR1, SbNPR1, and SbTGA) increased by salt treatment ( Figure 4 ). Also, SbMAP3K responded positively to salt stress. Photosynthesis-related genes, including SbPsa, SbPsb, and SbPet, were downregulated by salt stress.

The MS raw data (total ion current (TIC)) showed that there were 548 and 824 metabolites detected from all samples and quality controls under the negative and positive modes, respectively. PCA results showed that metabolisms differed between groups ( Supplementary Figures S2A, B ). XComparison analysis determined that there were 63 and 66 differential metabolites in STT vs. STC under negative and positive ion modes, respectively ( Supplementary Figures S2D, E ). As for SST vs. SSC, there were 40 and 46 differential metabolites under negative and positive ion modes, respectively ( Supplementary Figure S2C, D ).

In the negative ion mode, distinct profiles of these metabolites are shown by a heatmap in these four groups ( Figure 5A ). In salt-tolerant LRNK1, 32 and 31 metabolites were up- and downregulated, respectively, after salt treatment. As for the salt-sensitive LR2381, salt stress mediated 31 upregulated and nine downregulated metabolites ( Figure 5B ). The results of the KEGG enrichment analysis showed that plant hormone signal transduction (sbi04075) was significantly enriched, which was partially regulated by salicylic acid (Com_87_neg). Aside from the two metabolites, salt stress decreased the contents of 2-oxobutyric acid (Com_274_neg), which was involved in multiple pathways but was different in the two sorghum species. In addition, phenylalanine metabolism (map00360) and starch and sucrose metabolism (map00500) were active for energy homeostasis to cope with salt stress. Of these, five metabolites showed different profiles in STT vs. STC and SST vs. SSC ( Table 2 ; Figure 5C ), which is important.

In the positive ion mode, distinct metabolic profiling was found after salt treatment ( Figure 5D ). In LRNK1, 42 and 24 metabolites were up- and downregulated, respectively. In LR2381, salt stress mediated 37 upregulated and nine downregulated metabolites ( Figure 5E ). Metabolic profiling in salt-tolerant LRNK1, glutathione metabolism (map00480), and flavonoid biosynthesis (map00941) were significantly enriched in KEGG analysis, in which glutathione disulfide (Com_2094_pos) and gallocatechin (Com_1917_pos) act as the key regulators, respectively. As for the salt-sensitive LR2381, Betaine (Com_20_pos) in glycine, serine, and threonine metabolism and ABC transporters pathway were upregulated by salt stress. The expression of these three differential metabolites in different groups indicated that they were the key metabolites affected by salt stress ( Table 3 ; Figure 5F ).

To better understand the relationship between DEGs and differential metabolites, integration analysis of metabolomic and transcriptomic data was performed. Association results of the KEGG pathway showed that plant hormone signal transduction (sbi04075) was active in both LRNK1 and LR2381. Glutathione metabolism (sbi00480) was active in salt-tolerant LRNK1 but not enriched in LR2381 ( Figures 6A–D ). This demonstrates that plant hormone signal transduction (sbi04075) and glutathione metabolism (sbi00480) are the key pathways responding to salt stress.