Transcriptome analysis reveals candidate genes related to phosphorus starvation tolerance in sorghum

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.


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 [1,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 [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 [4,5]. The global rock phosphate reserve is a nonrenewable resource, and approximately 70% of cultivated lands worldwide suffer from P deficiency [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 [7], exuding organic acids and phosphatases [6], and establishing symbiotic relationships with arbuscular mycorrhizal fungi [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 [9][10][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 [12][13][14][15]. Moreover, many transcription factors (TFs) could regulate the expression of PHT1;1, such as phosphate starvation response 1 (PHR1) [16], WRKY75 [17], WRKY45 [18] and MYB domain protein 62 (MYB62) [19]. The PHO1 family plays an important role in Pi translocation from root to shoot [20][21][22] and is downregulated by WRKY6 and WRKY42 [23,24]. Additionally, miR399, miR827 [25][26][27] and the zinc finger TF ZAT6 [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 [29].
Rice, maize, wheat and sorghum are important graminaceous crops. Among them, P uptake by sorghum was only lower than that by rice [30]. Furthermore, sorghum is not only a staple food for more than 500 million people but also a popular forage crop [29,31,32]. Given its high yield, extensive use and excellent adaptation to harsh environmental conditions, sorghum was planted specifically in arid and semiarid regions [33,34]. It has been suggested that the grain yield of sorghum is highly correlated with P levels [35], and approximately 0.2% of the dry weight of this crop is contributed by P [36]. Although many genes and signal pathways related to low-P stress have been identified in rice [37], wheat [38] and maize [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.

Pot experiments
To evaluate the P starvation tolerance of the 29 sorghum accessions, pot experiments were performed according to Hufnagel et al. [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. [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, 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 (http://cole-trapnell-lab.github.io/cufflinks/) [40,41]. DEGs were analyzed using Cuffdiff software [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 (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 (http://www.genome.jp/ kegg/pathway.html), which is a major public pathwayrelated 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 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. [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 [44,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. 1). Therefore, we selected accessions 12484 and 13443 as the low-P-tolerant and low-Psensitive 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 1). Moreover, the differences in RNRT, RSAR, RLR and RVR between the selected accessions 12484 and 13443 were significant (Fig. 2).
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 [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 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 [44,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. 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.

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 highthroughput sequencing technology. Approximately 512 million raw reads were generated, and 492 million clean reads were obtained after cleaning and quality checks ( Table 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

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. 4a). Moreover, 735 DEGs were upregulated and 414 DEGs were downregulated in both accessions (Fig. 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.
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. 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. 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. 5).  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) 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.

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 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 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.

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 4. Notably, almost all of the significantly enriched pathways were involved in metabolism.
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 4). Recent evidence has proven that hormones participate in the control of plants in response to P starvation [48,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. 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 Ethyleneinsensitive-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.

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 (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 5).

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-Ptolerant and -sensitive accessions in response to low-P stress (Fig. 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. 7b and c). Under low-P conditions, however, it significantly reduced the length of the primary root and the number of root tips (Fig. 7b and c).

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 [50,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-Ptolerant 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 [52][53][54]. In our study, both of the selected materials showed enhanced primary root elongation under low-P stress (Fig. 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 [55,56]. The SA signal-related gene OsAIM1 and GA signal-related gene OsSHB were also involved in primary root development [57,58], while the MYB type TF OsMYB1 was involved in lateral root development [59]. Root hairs play key roles in P uptake in plants [60]. In Arabidopsis, some TFs, including MYB types such as MEMBRANE ANCHORED MYB (maMYB) [61], bHLH types such as RHD SIX-LIKE 1 (RSL1) [62], RHD SIX-LIKE 3 (RSL3) [63], Lotus japonicus ROOT HAIR LESS-LIKE 2 (LRL2) [64] and HD-ZIP types such as HOMEODOMAIN GLABROUS 11 (HDG11) [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 [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 [67]. In our study, we found two EIL proteins, Sb02g043350 and Sb04g023730, and three bHLH-type TFs, Sb03g008290, Sb01g043570 and Sb10g005650, Fig. 6 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)    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-Ptolerant 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 [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 [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 [39], Arabidopsis [70] and rice [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 [72,73]. The same results were also reported in soybean under low-P stress [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 [74][75][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 [77], temperature [78], and nitrogen, P, or carbon nutrition [79,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 [81,82]. A kind of carotenoid derivative, strigolactone, plays roles in root responses to low-P stress in Arabidopsis [83] and tomato [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 [85,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 [43]. It has also been reported that Fe overloading modifies root growth in P-deprived plants [87][88][89]. Additionally, ROS is produced during low-P stress with the accumulation of Fe 3+ [85,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 [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 [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 [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 [93][94][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 [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 [97]. Some studies have suggested that SA is involved in the response of plants to low P by controlling their redox status [98,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 [100]. JA accumulation occurs when plants are exposed to a low-P environment, further enhancing herbivory resistance and reducing shoot growth [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 [102]. Substantial evidences have proved that TFs are also involved in phosphate homeostasis [103][104][105][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 [107]. Overexpression of OsMYB4P could activate the expression of several Pht genes and increase phosphate acquisition [108]. WRKY74 modulates the tolerance to phosphate starvation in rice [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 [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 starvationtolerant 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
Additional file 1: Table S1. Gene-specific primers used in gene expression analysis by qRT-PCR. (XLSX 13 kb) Additional file 2: 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)