Abstract Background Phosphorus (P) deficiency in soil is a worldwide issue and a major constraint on the production of sorghum, which is an important staple food, forage and energy crop. The depletion of P reserves and the increasing price of P fertilizer make fertilizer application impractical, especially in developing countries. Therefore, identifying sorghum accessions with low-P tolerance and understanding the underlying molecular basis for this tolerance will facilitate the breeding of P-efficient plants, thereby resolving the P crisis in sorghum farming. However, knowledge in these areas is very limited. Results The 29 sorghum accessions used in this study demonstrated great variability in their tolerance to low-P stress. The internal P content in the shoot was correlated with P tolerance. A low-P-tolerant accession and a low-P-sensitive accession were chosen for RNA-seq analysis to identify potential underlying molecular mechanisms. A total of 2089 candidate genes related to P starvation tolerance were revealed and found to be enriched in 11 pathways. Gene Ontology (GO) enrichment analyses showed that the candidate genes were associated with oxidoreductase activity. In addition, further study showed that malate affected the length of the primary root and the number of tips in sorghum suffering from low-P stress. Conclusions Our results show that acquisition of P from soil contributes to low-P tolerance in different sorghum accessions; however, the underlying molecular mechanism is complicated. Plant hormone (including auxin, ethylene, jasmonic acid, salicylic acid and abscisic acid) signal transduction related genes and many transcriptional factors were found to be involved in low-P tolerance in sorghum. The identified accessions will be useful for breeding new sorghum varieties with enhanced P starvation tolerance. Electronic supplementary material The online version of this article (10.1186/s12870-019-1914-8) contains supplementary material, which is available to authorized users. Keywords: Sorghum, Transcriptome analysis, Phosphorus starvation, Candidate genes, Malate, Root development Background Phosphorus (P) is an essential macronutrient for plant growth and development and is considered to regulate energy metabolism, enzymatic reactions, and signal transduction processes [[37]1, [38]2]. Inorganic orthophosphate (Pi) is the primary source of P for plants. However, Pi can be readily fixed with aluminum and iron in acidic soil and with calcium in alkaline soil [[39]3]. Since Pi content in soil is too low to satisfy the requirements for plant growth, excessive quantities of P fertilizer must be applied, which possibly cause environmental pollution [[40]4, [41]5]. The global rock phosphate reserve is a nonrenewable resource, and approximately 70% of cultivated lands worldwide suffer from P deficiency [[42]6]. Therefore, the development of plants that are better adapted to low-P environments is a sustainable and economical approach in agricultural production. To adapt to persistent P deficiency, plants have developed several strategies, including changing root morphology [[43]7], exuding organic acids and phosphatases [[44]6], and establishing symbiotic relationships with arbuscular mycorrhizal fungi [[45]8]. These strategies in plants are dependent on changes in gene expression. Numerous genes related to low-P stress have been identified and recently reported in plants [[46]9–[47]11]. The phosphate transporter1 (PHT1) gene family participates in the uptake of Pi from the soil. As a member of the PHT1 family, PHT1;1 plays an important role in Pi uptake from the soil [[48]12–[49]15]. Moreover, many transcription factors (TFs) could regulate the expression of PHT1;1, such as phosphate starvation response 1 (PHR1) [[50]16], WRKY75 [[51]17], WRKY45 [[52]18] and MYB domain protein 62 (MYB62) [[53]19]. The PHO1 family plays an important role in Pi translocation from root to shoot [[54]20–[55]22] and is downregulated by WRKY6 and WRKY42 [[56]23, [57]24]. Additionally, miR399, miR827 [[58]25–[59]27] and the zinc finger TF ZAT6 [[60]28] modulate phosphate homeostasis. Most of these factors were identified in model plants, such as Arabidopsis and rice. In sorghum, however, only the multiple homologs of the phosphorus starvation tolerance 1 (PSTOL1) gene have been identified, which may be associated with changes in root morphology and root system architecture under low-P conditions [[61]29]. Rice, maize, wheat and sorghum are important graminaceous crops. Among them, P uptake by sorghum was only lower than that by rice [[62]30]. Furthermore, sorghum is not only a staple food for more than 500 million people but also a popular forage crop [[63]29, [64]31, [65]32]. Given its high yield, extensive use and excellent adaptation to harsh environmental conditions, sorghum was planted specifically in arid and semiarid regions [[66]33, [67]34]. It has been suggested that the grain yield of sorghum is highly correlated with P levels [[68]35], and approximately 0.2% of the dry weight of this crop is contributed by P [[69]36]. Although many genes and signal pathways related to low-P stress have been identified in rice [[70]37], wheat [[71]38] and maize [[72]39] through transcript profiling, the genes and pathways related to low-P stress in sorghum remain unclear. In this study, the P starvation tolerance of 29 sorghum accessions was evaluated. Accessions 12484 and 13443 were identified as low-P-tolerant and low-P-sensitive accessions, respectively. The gene expression in the roots of the two accessions, which both underwent P treatment for 8 days, was analyzed by mRNA sequencing (RNA-seq). Differentially expressed genes (DEGs) related to P starvation tolerance were identified through transcriptome profiles. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed. Moreover, the effect of malate application on root development in sorghum was investigated. The findings reported in this work increase our understanding of the molecular mechanisms of P starvation in sorghum. Methods Pot experiments To evaluate the P starvation tolerance of the 29 sorghum accessions, pot experiments were performed according to Hufnagel et al. [[73]29]. Generally, the sorghum seeds were sown in plastic pots in a greenhouse in the experimental station of Nanjing Agricultural University in 2015 and 2016. The Olsen-P of soil used in the experiments was 1.25 μg/g. After 7 days, two accessions were allowed to continue growing under sufficient-P (30 μg/g) and low-P (3 μg/g) conditions in each pot for 40 days. The roots of sorghum were harvested and cleaned carefully after breaking the pots. Subsequently, electronic images of root morphology were captured using an EPSON scanner and analyzed by WinRHIZO software. Roots and shoots were heated at 105 °C for 30 min and then oven-dried at 65 °C for 72 h. The samples were weighed and milled. A portion of the milled samples was digested with HNO[3]. The total P content of shoots was determined by ICP-OES (PerkinElmer, Optima 8000, USA). The relative length of roots (RLR) was expressed by the ratio of the root length of plants treated with low P to that treated with sufficient P. The relative surface area of root (RSAR), relative average diameter of root (RADR), relative volume of root (RVR), relative number of root tips (RNRT), relative dry weight of root (RDWR), relative dry weight of shoot (RDWS) and relative P content of shoot (RPCS) were calculated using similar methods. Hydroponic experiments To determine the duration of the low-P treatment and facilitate the harvest of roots for RNA-seq, hydroponic experiments were performed according to Hufnagel et al. [[74]29]. Sorghum seeds were sterilized with 75% ethanol for 5 min, dried, and then germinated in moistened filter paper. After 3 days, seedlings with similar growth vigor were transplanted into distilled water. When the seedlings reached the three-leaf stage, they were subjected to sufficient-P conditions (Hoagland (5 mmol/L Ca(NO[3])[2]·4H[2]O, 5 mmol/L KNO[3], 2 mmol/L MgSO[4]·7H[2]O, 0.025 mmol/L Fe-EDTA, 1 μmol/LH[3]BO[3], 1 μmol/L MnSO[4]·H[2]O, 1 μmol/L ZnSO[4]·7H[2]O, 0.5 μmol/L CuSO[4]·5H[2]O, and 0.005 μmol/L (NH[4])[6]Mo[7]O[24]·4H[2]O) with 1.0 mmol/L KH[2]PO[4]) and low-P conditions (Hoagland with 1 μmol/L KH[2]PO[4] and 1.0 mmol/L KCl) for 2 weeks. The nutrient solution with a pH of 5.8 ± 0.1 was replaced every 3 days. The seedlings were cultured in a growth chamber under a 12 h/12 h (day/night) photoperiod and a temperature cycle of 28 °C/22 °C (day/night). Each treatment was replicated three times. Roots and shoots of plants subjected to treatment with different P concentrations (sufficient-P and low-P conditions) were harvested every 2 days and used for measuring the ratio of root dry weight to shoot dry weight (R/S ratio) and gene expression analysis. RNA isolation and transcriptome sequencing Total RNA was isolated using TRIZOL (Invitrogen, USA) from the roots of accessions 12484 and 13443 harvested after 8 days of P treatment. The quality and integrity of the RNA were checked by an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The total RNA concentration was assessed using a QUBIT RNA ASSAY KIT (Invitrogen, USA). mRNA was enriched by the NEBNext® Poly(A) mRNA Magnetic Isolation Module (Invitrogen, USA), and the mRNA molecules were fragmented and subsequently used in first- and second-strand cDNA syntheses. The cDNA was subsequently subjected to terminal repair and poly (A) and unique adapter ligation. Prior to sequencing, the cDNA fragments were amplified and purified. The purified amplification products were sequenced on an Illumina Hiseq 2500 platform (CapitalBio Technology, Beijing, China). RNA-seq data analysis The original sequencing data were defined as raw reads. The clean reads were generated from the raw reads after removing the low-quality reads, mismatches, and adaptor sequences. The sequences of the clean reads were aligned with the sorghum transcript sequence (Ensembl, [75]http://plants.ensembl.org/index.html). The perfectly matching sequences were used for further analyses. The gene expression levels were normalized by fragments per kilo base of transcript per million fragments mapped (FPKM) and analyzed using Cufflinks software ([76]http://cole-trapnell-lab.github.io/cufflinks/) [[77]40, [78]41]. DEGs were analyzed using Cuffdiff software [[79]42]. The DEGs were identified based on the following criteria: the absolute value of the Log fold change had to be ≥1, the P value adjusted using the false discovery rate (FDR) method had to be ≤0.05, and three biological replicates had to be used. DEGs were subjected to GO analyses, as follows: 1) all DEGs were mapped to GO terms in the database ([80]http://www.geneontology.org/); 2) the gene numbers of each term were calculated; 3) the hypergeometric test was used to analyze the significantly enriched GO terms in the DEGs; 4) the input frequency represented the ratio of the number of DEGs with annotation to the total number of DEGs, and the background frequency represented the ratio of the number of genes with annotation to the total number of genes. To analyze the pathways that were significantly associated with DEGs, we used the same method to blast the DEGs with the KEGG database ([81]http://www.genome.jp/kegg/pathway.html), which is a major public pathway-related database. Then, the gene numbers of each pathway were calculated, and the hypergeometric test was used to analyze the significantly enriched pathways. Quantitative real-time polymerase chain reaction (qRT-PCR) Total RNA was extracted from each sample by using the Plant RNA Extract Kit (TIANGEN, Beijing, China), and approximately 0.5 μg of total RNA was used for cDNA synthesis using HiScript® II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China). After diluting the cDNA reaction mixture five times, 2 μL of the reaction mixture was used as template in a 20-μL reaction system. In addition, the reaction system contained 0.8 μL of 10 μmol/L gene-specific primers (Additional file [82]1: Table S1) and 10 μL of AceQ® qPCR SYBR® Green Master Mix (Low ROX Premixed) (Vazyme, Nanjing, China). qRT-PCR was performed on an ABI 7500 real-time PCR system (Applied Biosystems, Forster City, CA, USA), and the data were analyzed using the ABI 7500 Sequence Detection System software v.1.4. A sorghum constitutive expression gene, 18S rRNA, was used as the reference gene for normalization. Three technical replicates were carried out for each sample. Malate application experiments Accession 13443 was used as the experimental material. Hydroponic experiments were performed according to Mora-Macías et al. [[83]43]. After seeds germinated in distilled water for 1 week, the seedlings displaying similar growth vigor were picked and subjected to sufficient P (1.0 mmol/L KH[2]PO[4]) with and without 1.0 mmol/L malate or low P (1.0 μmol/L KH[2]PO[4]) with and without malate for 2 weeks. Each treatment was replicated three times. The lengths of the primary root and root tips were measured. Results Identification of the low-P-tolerant accession (12484) and the low-P-sensitive accession (13443) It is well documented that the biomass of the shoot is greatly reduced by P-deficient treatment [[84]44, [85]45], and shoot biomass is a good indicator of plant response to P deficiency tolerance. We expected that the shoot biomass of the low-P-tolerant accession would be less negatively affected by low-P treatment than by sufficient-P treatment; thus, the low-P-tolerant accession would display a high RDWS, and vice versa for the low-P-sensitive accession. A total of 29 sorghum accessions were grown in both low-P and sufficient-P conditions. These accessions displayed great variability in RDWS with values ranging from 0.2 to 0.9. The RDWS results were generally consistent between two independent experiments carried out in two consecutive years (2015 and 2016). Among them, accession 12484 displayed the highest RDWS, while accession 13443 showed the lowest values in both years (Fig. [86]1). Therefore, we selected accessions 12484 and 13443 as the low-P-tolerant and low-P-sensitive accessions, respectively, for further analysis. Meanwhile, the root morphology of each accession was tested. The RPCS showed a significant positive correlation with RDWS, and there were significant positive correlations with RDWR, RNRT, RSAR, RLR, and RVR (Table [87]1). Moreover, the differences in RNRT, RSAR, RLR and RVR between the selected accessions 12484 and 13443 were significant (Fig. [88]2). Fig. 1. [89]Fig. 1 [90]Open in a new tab Physiological changes of sorghum accessions in response to low-P stress. a and b Relative dry weight of shoot (a) and relative P content of shoot (b) of different accessions from 2015 and 2016 experiments. Asterisks denote statistically significant differences between the results from two independent experiments (P < 0.05). Bars represent the standard error of the mean (n = 3) Table 1. Correlation analysis between root morphology and P starvation tolerance. Asterisks note statistically significant differences between control and treatment groups (P < 0.05), double asterisks note statistically significant differences between control and treatment groups (P < 0.01). Bars represent the standard error of the mean (n = 3) 2015/2016 Relative length of root Relative surface area of root Relative average diam of root Relative volume of root Relative number of root tips Relative dry weight of root Relative dry weight of shoot Relative P content of shoot Relative length of root 1 0.893^**/0.948^** 0.101/0.137 0.719^**/0.812^** 0.918^**/0.840^** 0.969^**/0.422^** 0.539^**/0.475^** 0.572^**/0.507^** Relative surface area of root 0.893^**/0.948^** 1 0.486^**/0.435^** 0.953^**/0.954^** 0.695^**/0.694^** 0.977^**/0.529^** 0.591^**/0.552^** 0.521^**/0.578^** Relative average diam of root 0.101/0.137 0.486^**/0.435^** 1 0.678^**/0.665^** −0.129/−0.208 0.315/0.444^* 0.429^*/0.360 0.138/0.342 Relative volume of root 0.719^**/0.812^** 0.953^**/0.954^** 0.678^**/0.665^** 1 0.468^**/0.488^** 0.868^**/0.587^** 0.543^**/0.567^** 0.426^**/0.596^** Relative number of root tips 0.918^**/0.840^** 0.695^**/0.694^** −0.129/−0.208 0.468^**/0.488^** 1 0.821^**/0.405^* 0.485^**/0.495^** 0.542^**/0.524^** Relative dry weight of root 0.969^**/0.422^** 0.977^**/0.529^** 0.315/0.444^* 0.868^**/0.587^** 0.821^**/0.405^* 1 0.582^**/0.647^** 0.560^**/0.687^** Relative dry weight of shoot 0.539^**/0.475^** 0.591^**/0.552^** 0.429^*/0.360 0.543^**/0.567^** 0.485^**/0.495^** 0.582^**/0.647^** 1 0.797^**/0.876^** Relative P content of shoot 0.572^**/0.507^** 0.521^**/0.578^** 0.138/0.342 0.426^**/0.596^** 0.542^**/0.524^** 0.560^**/0.687^** 0.797^**/0.876^** 1 [91]Open in a new tab * menas P < 0.05; ** means P < 0.01 Fig. 2. Fig. 2 [92]Open in a new tab Effect of low-P stress on root morphologies of sorghum accessions 12484 and 13443 grown in soil. Asterisks denote statistically significant differences between the two accessions (P < 0.05) by Student’s t-test. Bars represent the standard error of the mean (n = 3) Increased P content in the shoot accounts for enhanced tolerance to low P Plants can achieve tolerance to low P availability by optimizing the P utilization and/or by improving P acquisition from the soil [[93]46]. To understand whether the variations in low-P tolerance of the 29 sorghum accessions were due to their improved ability to acquire P, we determined the P content in the shoot tissues of these accessions grown under both low-P and sufficient-P conditions. Figure [94]1b shows that low-P-tolerant accessions displayed high P content, while sensitive accessions displayed low P content. In general, there was a good correlation between RPCS and RDWS, with correlation coefficients equal to 0.797 and 0.876 in 2015 and 2016, respectively, indicating variation in P acquisition ability accounts for the variation in low-P tolerance in these selected sorghum accessions. Determining optimal sampling time for RNA-seq To analyze the transcript profiles of low-P-tolerant and -sensitive accessions in responding to low P, we first determined the optimal sampling time by analyzing the R/S of sorghum grown hydroponically under low P or sufficient P supply. An increasing R/S ratio is a hallmark indicator of plant responses to P-deficient conditions [[95]44, [96]47]. The R/S ratio of accession 12484 seedlings grown under low-P conditions was not significantly different than that grown under sufficient-P conditions after six successive days of treatment (Fig. [97]3). In contrast, from the 8th to the 14th day after low-P treatment, the R/S ratio was significantly higher than that of plants grown under sufficient-P conditions. Moreover, a similar trend was found in accession 13443. Therefore, the roots of seedlings grown under low P and sufficient P for 8 days were sampled for RNA-seq analysis. Fig. 3. [98]Fig. 3 [99]Open in a new tab Physiological responses of sorghum accessions 12484 and 13443 to low-P stress in the hydroponic system. a and b Root/shoot dry weight ratio of accession 12484 (a) and accession 13443 (b). c Photographs of representative plants. Asterisks denote significant differences between control and treatment groups by Student’s t-test (* means P < 0.05; ** means P < 0.01). Bars represent the standard error of the mean (n = 3) RNA-seq analysis and de novo assembly The libraries from the roots of accessions 12484 and 13443 grown under low-P and sufficient-P conditions for 8 days were sequenced using Illumina high-throughput sequencing technology. Approximately 512 million raw reads were generated, and 492 million clean reads were obtained after cleaning and quality checks (Table [100]2). The clean rates and the Q20 and Q30 values of the samples reached up to 96, 97, and 93%, respectively, indicating the high quality of the sequencing results. Approximately 87% of the total clean reads were mapped to the sorghum reference sequence Sorghum_bicolor.Sorbi1.33 (Ensembl; [101]http://plants.ensembl.org/index.html), > 97.5% of which were uniquely mapped reads. Table 2. Summary of RNA-seq data and de novo assembly accession Sample Raw reads number Clean reads number Clean rate(%) Q20(%) Q30(%) Mapped reads rate(%) Uniquely mapped reads rate(%) Gene Number Transcript Number 13443 sufficient P_1 44920566 43070114 95.88 97.62 93.31 89.36 97.45 25489 33223 sufficient P_2 42063644 40616024 96.56 97.76 93.63 85.21 97.43 26110 33633 sufficient P_3 30906152 29804314 96.43 97.73 93.55 88.59 97.54 24579 31528 low P_2 37808338 36157312 95.63 97.6 93.28 88.36 97.62 24765 31741 low P_3 48044802 46225956 96.21 97.71 93.54 87.75 97.66 25942 33394 low P_5 51611204 49267030 95.46 97.54 93.16 89.32 97.41 25331 33283 12484 sufficient P_1 41036008 39463806 96.17 97.69 93.48 86.76 97.46 25018 32024 sufficient P_2 44683368 42943164 96.11 97.71 93.56 84.96 97.37 25196 32012 sufficient P_3 43627358 41878828 95.99 97.67 93.46 84.02 97.46 25308 32157 low P_1 41144016 39708924 96.51 97.76 93.63 88.35 97.55 24587 31120 low P_2 46931232 45299694 96.52 97.75 93.61 86.93 97.52 24935 31395 low P_5 39105304 37577298 96.09 97.65 93.39 87.27 97.61 24670 31545 [102]Open in a new tab Screening of the candidate genes related to P starvation tolerance in sorghum Comparing the transcripts of sorghum grown under low-P conditions with those under sufficient-P conditions, a total of 2627 DEGs in accession 13443 and 5240 DEGs in accession 12484 were identified. Among these DEGs, 1815 were upregulated and 812 were downregulated in accession 13443, while 2211 and 3029 DEGs were upregulated and downregulated in accession 12484, respectively (Fig. [103]4a). Moreover, 735 DEGs were upregulated and 414 DEGs were downregulated in both accessions (Fig. [104]4b). Interestingly, the numbers of DEGs, regardless of their upregulation or downregulation, in accession 13443 were lower than those in accession 12484. Here, the 5275 nonredundant DEGs and the 1296 overlapping DEGs, for a total of 6571 DEGs, were identified as DEGs in response to low-P stress in sorghum. Fig. 4. [105]Fig. 4 [106]Open in a new tab Comparison of the number of DEGs. 1: low-P-sensitive accession 13443; 2: low-P-tolerant accession 12484; L: low P; C: sufficient P. a The number of differentially expressed genes in each part. b Venn diagram illustrating the genes of the two accessions in response to low-P stress. c Response of DEGs to low-P stress and different accessions under low-P conditions A total of 4762 nonredundant DEGs, of which 2237 were upregulated and 2525 were downregulated, were identified in the low-P-tolerant accession 12484 vs. the low-P-sensitive accession 13443 under low-P conditions (Fig. [107]4a). During our study, these 4762 nonredundant DEGs were identified as DEGs in different accessions under low-P conditions. Finally, comparing these two groups, a total of 2089 common DEGs were found (Fig. [108]4c), meaning these genes were not only differentially expressed in response to low-P stress but also differentially expressed in different accessions under low-P stress. These DEGs were identified as the candidate genes related to P starvation tolerance in sorghum. qRT-PCR verification To verify the RNA-seq results, we determined the abundance of 20 randomly selected DEGs by qRT-PCR assay. The qRT-PCR assay showed that Sb05g022855, Sb05g026550, Sb02g007580, Sb01g023270, Sb09g024950, Sb01g002580, Sb08g007700, Sb06g027670, Sb04g007280, Sb09g026280 and Sb10g022080 were upregulated, whereas Sb06g021250, Sb03g026840, Sb07g022650, Sb01g016730, Sb03g009010, Sb01g034670, Sb06g031920 and Sb09g022260 were downregulated in both accessions (Fig. [109]5). Additionally, Sb01g042040 was downregulated in accession 12848 but was upregulated in accession 13443. These results were consistent with the data obtained from the RNA-seq analysis, indicating the reliability of the RNA-seq results. Fig. 5. Fig. 5 [110]Open in a new tab Validation of transcript abundance obtained from RNA-seq using qRT-PCR. Twenty randomly chosen genes were used for validation. a and b Relative expression by RNA-seq and qRT-PCR of accession 13443 (a) and accession 12484 (b). Bars represent the standard error of the mean (n = 3) Annotation and function classification of DEGs To evaluate the potential functions of those DEGs, GO analyses were performed on the DEGs identified during each comparison. Among the DEGs identified in response to low-P stress, 1578 of the 2627 DEGs in accession 13443 were assigned to at least one GO term, whereas 3147 of the 5240 DEGs in accession 12484 were assigned to at least one GO term. Among the DEGs in different accessions under low-P conditions, 2317 of the 4762 DEGs were assigned to different GO terms. Among the candidate DEGs, 1006 of the 2089 DEGs were assigned to at least one GO term. To deeply investigate the functions of DEGs, GO enrichment of each group was analyzed, and the top 30 GO enrichment terms were listed (Additional file [111]2: Figure S1). The results showed the following: 1) For the DEGs in the different accessions under low-P conditions, many of the top 30 significantly enriched GO terms were in the molecular function category, suggesting these two accessions had significant differences in molecular functions under low-P conditions. 2) Candidate genes were significantly enriched in the molecular function and biological process categories, implying that molecular functions and biological processes play important roles in P starvation tolerance of sorghum. Within the molecular function category, GO terms related to oxidoreductase activity (GO:0016491, GO:0016705, GO:0016684, GO:0051213, GO:0004601, GO:0004497) were significantly enriched. The candidate genes enriched in GO terms for which the Input frequency/Background frequency ratio was higher than 1 were examined and are listed in Table [112]3. The results showed that a total of 20 GO terms exhibited a higher Input frequency than Background frequency, suggesting that the expression changes of these genes might affect the functions of these GO terms. Among them, GO terms related to antioxidant activity, transporter activity, nucleic acid binding TF activity, catalytic activity, extracellular region, membrane part and response to stimulus also showed a higher Input frequency than Background frequency in the DEGs in response to low-P stress and the DEGs in different accessions under low-P conditions. Table 3. List of the focused significantly enrichedment GO terms Namespace ID Term Input frequency/Background frequency 13443(L/C) 12484(L/C) 12484 L vs 13443 L candidate genes molecular function GO:0016209 antioxidant activity 1.9422 2.0953 1.8839 2.6634 molecular function GO:0005215 transporter activity 1.2064 1.3471 1.2791 1.3071 molecular function GO:0001071 nucleic acid binding transcription factor activity 1.3734 1.3157 1.2425 1.2459 molecular function GO:0003824 catalytic activity 1.1786 1.1433 1.1391 1.2145 cellular component GO:0005576 extracellular region 1.1444 1.8283 1.7582 2.4556 cellular component GO:0044425 membrane part 1.3049 1.2736 1.0791 1.0349 biological process GO:0050896 response to stimulus 1.0837 1.1235 1.0924 1.0360 molecular function GO:0060089 molecular transducer activity 0.8920 1.0287 1.0935 1.2583 molecular function GO:0030234 enzyme regulator activity 0.8054 1.0453 0.8389 1.1622 molecular function GO:0009055 electron carrier activity 1.0321 0.9704 0.9665 1.1497 molecular function GO:0098772 molecular function regulator 0.8312 1.0748 0.8044 1.1363 molecular function GO:0045735 nutrient reservoir activity 0.4629 1.7021 0.9457 1.1098 cellular component GO:0016020 membrane 1.1236 1.1582 0.9954 1.0698 biological process GO:0001906 cell killing – – – 2.9594 biological process GO:0048511 rhythmic process 0.5311 1.3314 1.6275 2.4213 biological process GO:0040011 locomotion 1.1418 1.0370 – 1.6646 biological process GO:0051704 multi-organism process 0.9849 1.0975 0.9937 1.5259 biological process GO:0044699 single-organism process 1.0054 1.0551 0.9164 1.1612 biological process GO:0008152 metabolic process 0.9470 0.9064 0.8909 1.1021 biological process GO:0002376 immune system process 0.7807 1.1255 1.1963 1.0871 [113]Open in a new tab “—” means GO term not significantly enreached Significant enrichment pathways Genes usually play roles in certain biological processes by interacting with each other. Therefore, pathway analyses are helpful for further understanding the biological functions of genes. Pathway enrichment analysis of each group’s DEGs based on the KEGG database was performed, and the significantly enriched pathways are listed in Table [114]4. Notably, almost all of the significantly enriched pathways were involved in metabolism. Table 4. List of significantly enrichedment pathways Term Pathways Identifiers P-Value First class Second class 13443 (L/C) 12484 (L/C) 12484 L vs 13443 L candidate DEGs Phenylpropanoid biosynthesis map00940 4.2902E-06 6.0067E-06 0.0001 6.3352E-09 Metabolism Biosynthesis of other secondary metabolites Phenylalanine metabolism map00360 2.5179E-05 4.1713E-06 0.0010 1.2791E-06 Metabolism Amino acid metabolism Biosynthesis of secondary metabolites map01110 0.0002 0.0001 0.0001 2.2959E-06 Metabolism Global and overview maps Metabolic pathways map01100 – 0.0023 0.0084 0.0068 Metabolism Global and overview maps Flavone and flavonol biosynthesis map00944 0.0278 – 0.0244 0.0116 Metabolism Biosynthesis of other secondary metabolites Carotenoid biosynthesis map00906 – 0.0038 0.0188 0.0189 Metabolism Metabolism of terpenoids and polyketides Starch and sucrose metabolism map00500 – 0.0010 – 0.0217 Metabolism Carbohydrate metabolism Monoterpenoid biosynthesis map00902 – 0.0218 – 0.0223 Metabolism Metabolism of terpenoids and polyketides Galactose metabolism map00052 – 0.0168 0.0056 0.0368 Metabolism Carbohydrate metabolism Flavonoid biosynthesis map00941 – 0.0398 – 0.0403 Metabolism Biosynthesis of other secondary metabolites Nitrogen metabolism map00910 – 0.0359 – 0.0483 Metabolism Energy metabolism Glutathione metabolism map00480 – 0.0016 – – Metabolism Metabolism of other amino acids Amino sugar and nucleotide sugar metabolism map00520 – 0.0045 – – Metabolism Carbohydrate metabolism Plant hormone signal transduction map04075 – 0.0082 – – Environmental Information Processing Signal transduction alpha-Linolenic acid metabolism map00592 – 0.0250 – – Metabolism Lipid metabolism Zeatin biosynthesis map00908 – 0.0302 – – Metabolism Metabolism of terpenoids and polyketides Glycerolipid metabolism map00561 0.0004 – – – Metabolism Lipid metabolism Cutin, suberine and wax biosynthesis map00073 0.0066 – – – Metabolism Lipid metabolism Phenylalanine, tyrosine and tryptophan biosynthesis map00400 0.0213 – – – Metabolism Amino acid metabolism Taurine and hypotaurine metabolism map00430 0.0277 – – – Metabolism Metabolism of other amino acids Degradation of aromatic compounds map01220 0.0278 – – – Metabolism Global and overview maps Stilbenoid, diarylheptanoid and gingerol biosynthesis map00945 0.0362 – – – Metabolism Biosynthesis of other secondary metabolites Glycerophospholipid metabolism map00564 0.0481 – – – Metabolism Lipid metabolism Cysteine and methionine metabolism map00270 – – 0.0187 – Metabolism Amino acid metabolism Thiamine metabolism map00730 – – 0.0260 – Metabolism Metabolism of cofactors and vitamins [115]Open in a new tab “—” means that the P-Value is greater than 0.05 For those pathways, 1) 3 pathways (phenylpropanoid biosynthesis, phenylalanine metabolism, and biosynthesis of secondary metabolites) were significantly enriched by low-P treatment in both accessions; 2) pathways specifically enriched in accession 12484 by low-P stress include sugar metabolism (starch and sucrose, galactose, amino sugar and nucleotide sugar) and carotenoid biosynthesis, while those in accession 13443 include amino acid metabolism (phenylalanine, tyrosine, tryptophan, taurine and hypotaurine) and lipid metabolism (glycerolipid and glycerophospholipid); 3) for the candidate DEGs, 11 pathways were significantly enriched, and most of them were only significantly enriched in accession 12484 in response to low-P stress. Candidate genes involved in plant hormone signal transduction Almost all of the significantly enriched pathways were involved in metabolism. Interestingly, only one term (plant hormone signal transduction) involved in signal transduction was from the class of environmental information processing (Table [116]4). Recent evidence has proven that hormones participate in the control of plants in response to P starvation [[117]48, [118]49]. Therefore, the DEGs related to plant hormone signal transduction were chosen for further analysis. A total of 21 candidate DEGs were found to be associated with plant hormone (including auxin (AUX), abscisic acid (ABA), jasmonic acid (JA), ethylene (ET), and salicylic acid (SA)) signal transduction (Fig. [119]6). Among them, five DEGs were AUX signaling pathway-related genes, including two small auxin upregulated RNA (SAUR) genes and three Aux/IAA genes. Moreover, four candidate genes were ABA signaling pathway-related genes, including one protein phosphatase 2C (PP2C) and three PYR/PYL genes. Furthermore, four candidate genes were involved in the JA signaling pathway, including three jasmonatezim-domain (JAZ) genes and one coronatine insensitive 1 (COI1) gene. In addition, there were seven candidate genes related to the ET signaling pathway, including one EBF1/2 gene, two EIN3 genes, two Ethylene-insensitive-3-like (EIL) genes and two ETR genes, and three candidate genes related to the SA signaling pathway, including one PR-1 gene and two TGA genes. Fig. 6. Fig. 6 [120]Open in a new tab Heatmap of the candidate genes involved in plant hormone signal transduction. The log2 fold change of the candidate genes involved in plant hormone signal transduction under low-P conditions compared with that under sufficient-P conditions in each section is represented by a color scale consisting of red (upregulated), white (not regulated) and green (downregulated) Transcriptional factors (TFs) identified from candidate genes A total of 127 nonredundant TFs were identified from the candidate DEGs using the plant TF database PlantTFDB version 4.0 ([121]http://planttfdb.cbi.pku.edu.cn/). In general, 109 were responding to low P in accession 12848, while 50 were responding to low P in accession 13443. Among them, bHLH, WRKY and MYB were the top three most abundant TFs, numbering 17, 13, and 12, respectively (Table [122]5). Table 5. List of genes belonging to TFs TF-Family Gene Expression (in response to low P) TF-Family Gene Expression (in response to low P) 13443 12484 13443 12484 AP2 Sb01g003400 up up HD-ZIP Sb01g022420 – up B3 Sb09g004870 – up HD-ZIP Sb06g025750 – up bHLH Sb03g008290 – up HD-ZIP Sb04g023410 – down bHLH Sb10g005650 – up HD-ZIP Sb06g024000 – down bHLH Sb02g027210 – up HSF Sb01g021490 up up bHLH Sb01g045710 – up HSF Sb01g008380 up – bHLH Sb01g043570 – up HSF Sb03g033750 – down bHLH Sb03g036100 up up LSD Sb07g004050 – up bHLH Sb07g004190 up up LBD Sb09g030780 up – bHLH Sb2250s002010 up – LBD Sb01g014800 – down bHLH Sb09g006220 up – MIKC_MADS Sb07g001250 – up bHLH Sb04g017390 down down MIKC_MADS Sb06g017660 – up bHLH Sb04g027280 – down MYB Sb03g004600 – up bHLH Sb01g041960 – down MYB Sb03g032260 – up bHLH Sb03g042860 – down MYB Sb02g040480 – up bHLH Sb05g023730 – down MYB Sb08g016620 up – bHLH Sb06g028750 – down MYB Sb02g030900 up – bHLH Sb06g020810 – down MYB Sb08g018840 down up bHLH Sb03g005250 – down MYB Sb03g003120 down down bZIP Sb08g020600 – up MYB Sb09g001590 – down bZIP Sb04g008840 – up MYB Sb03g012310 – down bZIP Sb04g007060 – up MYB Sb02g024260 – down bZIP Sb07g025490 – up MYB Sb02g043420 – down bZIP Sb01g037520 up up MYB Sb09g002680 – down bZIP Sb04g000300 up – MYB_related Sb09g029560 – up bZIP Sb09g021840 down down MYB_related Sb03g006440 – down bZIP Sb09g024290 – down MYB_related Sb04g031590 – down bZIP Sb03g040530 – down MYB_related Sb03g011280 – down bZIP Sb02g027410 – down NAC Sb04g036340 – up C2H2 Sb04g021440 – up NAC Sb03g035820 up – C2H2 Sb01g031900 up up NAC Sb02g028870 up down C2H2 Sb01g031920 up up NAC Sb04g023990 – down C2H2 Sb03g025790 – down NAC Sb01g048130 – down C2H2 Sb07g028010 – down NAC Sb01g006410 – down C2H2 Sb02g028220 – down NAC Sb07g005610 – down C2H2 Sb01g043960 – down NAC Sb06g017190 – down C2H2 Sb10g004570 – down NAC Sb08g022560 – down C3H Sb03g003110 down down NAC Sb02g043210 down down C3H Sb09g006050 down – NF-YB Sb02g038870 down down CO-like Sb04g029480 – down NF-YC Sb06g033380 up – CO-like Sb06g021480 down down Nin-like Sb04g002940 – down EIL Sb02g043350 – down RAV Sb03g031860 – down EIL Sb04g023730 down down TALE Sb04g008030 – down ERF Sb01g040280 – up TALE Sb05g003750 – down ERF Sb06g025890 up up TALE Sb02g002200 – down ERF Sb06g025900 up up TCP Sb06g023130 up – ERF Sb06g024540 up up TCP Sb10g008030 down down ERF Sb07g023803 up – TCP Sb02g030260 – down ERF Sb03g012890 up down Trihelix Sb03g013050 up – ERF Sb10g004580 down down WRKY Sb03g038170 – up ERF Sb03g042060 – down WRKY Sb10g025600 up up ERF Sb09g020690 – down WRKY Sb10g025590 up up ERF Sb01g044410 – down WRKY Sb09g029850 up – FAR1 Sb01g040730 – down WRKY Sb05g001170 up – G2-like Sb07g021290 up – WRKY Sb03g028530 up down G2-like Sb02g001600 up – WRKY Sb03g029920 up down G2-like Sb06g031970 down down WRKY Sb03g047350 down – G2-like Sb01g036680 – down WRKY Sb02g043030 down down G2-like Sb04g003140 – down WRKY Sb04g009800 – down GeBP Sb03g009480 down down WRKY Sb08g005080 – down GeBP Sb03g005180 down down WRKY Sb02g024760 – down GRAS Sb02g034550 – down WRKY Sb01g036180 – down GRAS Sb06g017860 – down ZF-HD Sb05g001690 up down GRF Sb04g034800 – up [123]Open in a new tab Interestingly, among the TFs, two WRKY (Sb03g028530 and Sb03g029920), one NAC (Sb02g028870), one ERF (Sb03g012890) and one ZF-HD (Sb05g00169) were downregulated in accession 12848 but were upregulated in accession 13443 by low-P stress, while the opposite was observed for one MYB (Sb08g018840). Root development of sorghum in response to low P with and without malate Under low-P stress, plants often secrete organic acid in roots. Among the candidate genes, three were involved in malate metabolism, encoding malate dehydrogenase (Sb07g023910), malate synthase (Sb06g020720) and malic enzyme (Sb09g005810) respectively. And all of them showed different expression patterns in low-P-tolerant and -sensitive accessions in response to low-P stress (Fig. [124]7a). Our results showed that the primary root length and the numbers of root tips were significantly correlative with P starvation tolerance in sorghum. Therefore, we wondered whether malate affects the primary root length and the number of root tips in sorghum. The results showed that under sufficient-P conditions, application of malate did not affect the length of the primary root and the number of tips significantly (Fig. [125]7b and c). Under low-P conditions, however, it significantly reduced the length of the primary root and the number of root tips (Fig. [126]7b and c). Fig. 7. [127]Fig. 7 [128]Open in a new tab Effect of malate application on root development. a Expression profile of the malate metabolism-related genes XLOC_030929 (Sb07g023910), XLOC_027815 (Sb06g020720) and XLOC_035509 (Sb09g005810). b Effect of applying malate on the length of the primary root. c Effect of applying malate on the number of root tips. The different lowercase letters denote significant differences among different treatments. Bars represent the standard error of the mean (n = 3) Discussion Different sorghum accessions displayed dramatic variability in their tolerance of P starvation, which is associated with shoot internal P content, indicating that P uptake contributes to P starvation tolerance variability. Root morphology is a key trait for optimizing the efficiency of P acquisition in plants [[129]50, [130]51]. Our results showed that P starvation tolerance in sorghum was significantly correlated with root morphology, implying roots played important roles in P starvation tolerance in sorghum. Therefore, the transcription profiles in sorghum roots were studied. In our study, the majority of the accessions displayed consistent results between experiments carried out in 2 years. However, several accessions, such as 357 and 2349, showed significant differences between the 2 years. This led to the ranking of the tolerance of different accessions to low P. We chose two accessions, 12484 and 13443, as low-P-tolerant and -sensitive varieties, respectively, because these two accessions were constantly at the two extremes of the tolerance ranking. The ranking variations of tolerance to low P were concurrent with those of internal P content. The variations in internal P content of the same accession in different years indicate that there are other environmental factors beyond genotype contributing to P acquisition. The two accessions showed great differences in root morphology under low-P stress It was suggested that most plant responses to P deficiency included the remodeling of root morphology, which involves suppressing primary root growth and enhancing the production of root hairs and lateral roots [[131]52–[132]54]. In our study, both of the selected materials showed enhanced primary root elongation under low-P stress (Fig. [133]3), while more lateral roots were promoted in tolerant accessions than sensitive accessions, implying that the genes regulating the elongation of the primary root and the enhancement of lateral roots play important roles in the response to low-P stress. It is well documented that root morphology is regulated by many signaling pathways. For instance, the AUX signal-related genes OsARF12 and OsTOP1 were involved in primary root development in rice [[134]55, [135]56]. The SA signal-related gene OsAIM1 and GA signal-related gene OsSHB were also involved in primary root development [[136]57, [137]58], while the MYB type TF OsMYB1 was involved in lateral root development [[138]59]. Root hairs play key roles in P uptake in plants [[139]60]. In Arabidopsis, some TFs, including MYB types such as MEMBRANE ANCHORED MYB (maMYB) [[140]61], bHLH types such as RHD SIX-LIKE 1 (RSL1) [[141]62], RHD SIX-LIKE 3 (RSL3) [[142]63], Lotus japonicus ROOT HAIR LESS-LIKE 2 (LRL2) [[143]64] and HD-ZIP types such as HOMEODOMAIN GLABROUS 11 (HDG11) [[144]65] were found to be involved in root hair growth. Moreover, P starvation-induced genes involved in root hair development have been identified. For instance, OsAUX1 encoding an auxin influx transporter could promote root hair elongation under low-P stress in rice [[145]66]. In Arabidopsis, the ET signaling gene EIN3 and its closest homolog EIL1 were demonstrated to be involved in P starvation-induced root hair development [[146]67]. In our study, we found two EIL proteins, Sb02g043350 and Sb04g023730, and three bHLH-type TFs, Sb03g008290, Sb01g043570 and Sb10g005650, similar to RSL1, RSL3 and LRL2, respectively, that were candidate genes, suggesting that these genes might regulate P tolerance by mediating root hair growth in sorghum. Thus, our data help reveal the related genes or mechanism controlling root morphology in sorghum under low-P stress. More active responses to P starvation in the low-P-tolerant accession Under different P treatments, the number of DEGs in the low-P-tolerant accession 12484 was nearly two-fold higher than that in the low-P-sensitive accession 13443, implying that the low-P-tolerant accession responded to low P more positively and dramatically. Moreover, the numbers of unique DEGs identified in accession 13443 or accession 12484 were greater than the number of overlapping DEGs, indicating that there are differences between the two accessions in response to low P. Similarly, many more DEGs in response to low P were found in low-P-tolerant accessions than in low-P-sensitive accessions in soybean [[147]68]. In maize, however, the number of DEGs in low-P-tolerant materials was lower than that in low-P-sensitive materials under low-P stress [[148]69]. Thus, it is not valid to evaluate whether an accession is tolerant of low P or not solely based on the number of DEGs under low-P conditions. Complex mechanism of P starvation tolerance in sorghum Among the candidate genes, the significantly enriched GO terms referred to molecular function, cellular component and biological process, implying that achieving P starvation tolerance required many processes. The top 30 significantly enriched GO terms of DEGs under low-P conditions suggested that the two accessions had differences in molecular functions under low-P conditions. Further analysis showed that the input frequency of antioxidant activity, catalytic activity, transporter activity and nucleic acid binding TF activity, belonging to molecular function, were higher than the background frequency, suggesting the key roles of genes belonging to these terms in low-P tolerance in sorghum. Candidate genes related to P starvation tolerance were enriched in many pathways. Notably, genes participating in phenylpropanoid biosynthesis, secondary metabolism, starch and sucrose metabolism and nitrogen metabolism also actively responded to low-P stress in maize [[149]39], Arabidopsis [[150]70] and rice [[151]71]. Phenylpropanoid biosynthesis, phenylalanine metabolism, and biosynthesis of secondary metabolites were the top 3 enriched pathways, suggesting these pathways responded actively when sorghum suffered from low-P stress. Phenylpropanoids contribute to all aspects of plant responses toward biotic and abiotic stimuli, and phenylpropanoids are not only indicators of plant stress responses under varied light, mineral treatment, and pest stresses but are also key mediators in plants [[152]72, [153]73]. The same results were also reported in soybean under low-P stress [[154]53]. Moreover, the candidate genes were significantly enriched in flavone and flavonol biosynthesis, carotenoid biosynthesis and flavonoid biosynthesis pathways. Flavonol is one of the main classes of phenylpropanoid pathway derivatives [[155]74–[156]76]. Flavonoids and carotenoids are widely known for their influence on the colors of plant tissues, which further contributes to plant fitness and food quality. The flavonoid pathway has been suggested to play important roles in protecting plants from oxidative stresses induced by drought [[157]77], temperature [[158]78], and nitrogen, P, or carbon nutrition [[159]79, [160]80]. Carotenoids are crucial for driving biological processes in plants, including the assembly of photosystems and light harvesting antenna complexes for photosynthesis and the regulation of growth and development [[161]81, [162]82]. A kind of carotenoid derivative, strigolactone, plays roles in root responses to low-P stress in Arabidopsis [[163]83] and tomato [[164]84]. In general, the results of our study show that many GO terms and pathways are enriched in sorghum suffering from low-P stress, implying that a complex mechanism of P starvation tolerance exists in sorghum. Malate play key roles in root development in sorghum under low-P stress Notably, we found that applying malate reduced the length of the primary root and the number of root tips under P stress, implying malate was involved in P starvation tolerance in sorghum by affecting root development. Recently, reactive oxygen species (ROS) were also found to trigger callose deposition, which further adjusts apical meristem activity in roots under low-P stress [[165]85, [166]86]. For the candidate DEGs, many GO terms related to oxidoreductase activity were found, which suggested that ROS generation and elimination were active in sorghum suffering from low-P stress. In our study, applying malate only affected root development under low-P conditions, suggesting that low-P stress might generate some factors that perhaps promote malate to affect root development. Moreover, in Arabidopsis, malate was found to inhibit root growth by adjusting the accumulation of Fe in plants under low-P stress, which is dependent on ALMT [[167]43]. It has also been reported that Fe overloading modifies root growth in P-deprived plants [[168]87–[169]89]. Additionally, ROS is produced during low-P stress with the accumulation of Fe^3+ [[170]85, [171]90]. Our GO enrichment results showed that many P candidate genes were involved in Fe binding. From these results, we speculate that malate plays key roles in root development in sorghum under low-P stress, perhaps through Fe and ROS. Plant hormone signal transduction and TFs related to P starvation tolerance in sorghum Plant hormones are involved in many biological processes in enhanced resistance to environmental stresses, diseases and pathogen infections [[172]91]. In the present study, many DEGs were related to plant hormone (including AUX, ABA, JA, ET, and SA) signal transduction. AUX has been suggested to stimulate root growth and lateral root proliferation upon P starvation [[173]7]. Notably, the 3 AUX/IAA genes (Sb03g001490, Sb03g035500, Sb06g023800) related to AUX signal transduction were downregulated in low-P-tolerant material by low-P stress but were not significantly changed in low-P-sensitive material. In Arabidopsis, the phloem-mobile transcript of IAA could target the root tip and then regulate root morphology [[174]92]. We speculate that the 3 AUX/IAA genes play similar roles in root morphology regulation, but the specific functions of these genes under low-P stress in sorghum should be further studied. ET was also involved in low-P stress. P starvation could alter ET biosynthesis in plants, while ET could regulate root growth under low-P stress [[175]93–[176]95]. In our study, several ET signal-related genes were considered to be related to low-P tolerance in sorghum. In Arabidopsis, the eto1 mutant overproducing ET exhibited reduced primary root growth and increased production of root hairs [[177]96]. Whether a similar regulatory mechanism also exists in sorghum should be further investigated. Although the exact roles of ABA and SA in plants suffering from low P are not clear, these hormones were found to respond to low P. In Arabidopsis, ABA biosynthesis mutants exposed to low-P stress displayed reduced expression of P starvation-responsive genes and accumulation of anthocyanins [[178]97]. Some studies have suggested that SA is involved in the response of plants to low P by controlling their redox status [[179]98, [180]99]. In the present study, two candidate bZIP genes (Sb03g040530, Sb09g021840) encoding TGA were involved in the SA signaling pathway. In Arabidopsis, type-II TGA TFs were demonstrated to be key activators of JA/ET-induced immune reactions [[181]100]. JA accumulation occurs when plants are exposed to a low-P environment, further enhancing herbivory resistance and reducing shoot growth [[182]101]. However, the role that SA plays in low-P stress is unclear. Thus, we speculate that SA can regulate P starvation through crosstalk with the JA/ET signaling pathway or other pathways, but the underlying mechanism should be further studied. TFs were suggested to play a key role in the regulation of gene expression at the transcriptional level by binding to DNA regulatory elements [[183]102]. Substantial evidences have proved that TFs are also involved in phosphate homeostasis [[184]103–[185]106]. In our study, besides ET- and SA-related TFs, many other TFs including bHLHs, WRKYs, MYBs, bZIPs, NACs and C2H2 zinc finger proteins, were identified. In rice, OsMYB2P-1 is involved in the regulation of P starvation responses and root architecture by suppressing or activating downstream genes [[186]107]. Overexpression of OsMYB4P could activate the expression of several Pht genes and increase phosphate acquisition [[187]108]. WRKY74 modulates the tolerance to phosphate starvation in rice [[188]109]. In wheat, TaZAT8, a C2H2-ZFP-type TF gene, plays critical roles in mediating wheat tolerance to P deprivation by regulating P acquisition, ROS homeostasis and root system establishment [[189]110]. Our results indicate that TFs may also play important roles in sorghum tolerance to low P, but their functions should be further studied. Conclusions Transcriptome analysis showed that a P starvation-tolerant accession exhibited more active responses to low-P stress. A total of 2089 genes were identified and enriched in many GO terms and pathways, suggesting that P starvation tolerance of sorghum is a complex mechanism. Malate significantly reduced the length of the primary root and numbers of root tips in sorghum suffering from low-P stress. Plant hormone (including AUX, ET, JA, SA and ABA) signal transduction-related genes and many TFs were found to be involved in low-P tolerance in sorghum. The findings reported herein increase our understanding of the molecular characteristics of sorghum tolerant of low-P stress. The identified accessions will be useful for inbreeding new sorghum varieties with a high P starvation tolerance. Additional files [190]Additional file 1:^ (14KB, xlsx) Table S1. Gene-specific primers used in gene expression analysis by qRT-PCR. (XLSX 13 kb) [191]Additional file 2:^ (1.3MB, pdf) Figure S1. GO annotations and GO enrichment of DEGs. a and e DEGs in accession 12484 in response to low-P stress. b and f DEGs in accession 13443 in response to low-P stress. c and g DEGs in different accessions under low-P conditions. d and h Active DEGs. DEGs in accession 12484 in response to low-P stress. (PDF 1359 kb) Acknowledgements