Plant growth

The seeds were first incubated in an oven at 45 °C for 5 days to break down seed dormancy. Then, the seeds were germinated at 37 °C over 3–4 days. The plantlets were grown in controlled greenhouse conditions at 28–30 °C and approximately 70–80% humidity. To avoid the block effect, the plants were randomly positioned using IRRISTAT software v4.0 (International Rice Research Institute (IRRI), Los Baños, Philippines) with three replicates. Individual plantlet was grown in drainable black plastic bags measuring 68 cm×16 cm in length and diameter and filled with fine river sand. The plants were grown in two different Yoshida nutrient conditions (Yoshida et al. 1971): one received a full Pi supplement (P0) of 320 μM, whereas the other received a low Pi supplement (P*) of 10 μM. Specifically, 0.05 g/L NaH₂PO₄.H₂O (320 μM P) was used in the full Pi medium, whereas 1.56 mg/L NaH₂PO₄.H₂O (10 μM P) was used to supplement in the low Pi medium. The plants were irrigated with full Pi medium (control) or with the low Pi medium (treatment) every three days for six weeks.

Sample harvesting and phenotyping experiment

The plants were harvested after six weeks. Any sand clinging to the plantlets’ root systems was carefully removed using light faucet water pressure. The plant samples were dried in an oven at 70 °C for up to one week until they reached a constant weight. Shoots and leaves, including the stem base, were weighed to obtain shoot weight (SHW), and a similar procedure was performed on the roots to obtain root weight (RTW).

The relative efficiency of Pi use (REP, %) was determined by dividing the plant dry mass (DM, g) under low Pi (10 µM) by the DM under high Pi (320 µM) according to Ozturk et al. (Ozturk et al. 2005).

The relative RTW/SHW ratio (RRSR) was obtained from the fraction of the RTW/SHW ratio under low Pi (10 µM) over the RTW/SHW ratio under high Pi (320 µM) (Ozturk et al. 2005).

The Pi content in rice shoots was analyzed using the vanadomolybdophosphoric acid colorimetric method (Rice et al. 2017). A standard curve was established with the corresponding vanadate–molybdate reagent and a range of Pi concentrations represented by hydrous KH₂PO₄. For each variety, 0.3 g of the homogeneous dry sample was subjected to dry ashing over a 6 h combustion process in a Muffle furnace (Nabertherm, New Castle, DE, USA). Wet ashing of the combusted sample by the acid hydrolysis method, which consisted of diluted HCl fuming (37%; Merck Millipore, Burlington, MA, USA) and concentrated HNO₃ (69%; Merck Millipore) was subsequently performed. Absorbance of the assay solution was read at 420 nm by a UV-1800 UV–VIS spectrophotometer (Shimadzu, Kyoto, Japan). From the obtained optical density value, the Pi concentration in each sample was determined using the standard curve with its linear regression line.

The Pi content in shoots (mg) was calculated from the multiplication of Pi concentration (mg/g) and SHW (g). To obtain the PUpE (mg Pi/g Pi) of each Pi environment, the corresponding shoot Pi content was divided by sufficient or insufficient Pi supplies (Neto et al. 2016). The relative PUpE (RPUpE) was deduced as the ratio between the PUpE under low Pi and that under high Pi.

The PPUE (g² SHW/mg Pi) was calculated as the fraction between the SHW and Pi concentration at a given Pi concentration in the rooting medium (Hammond et al. 2009). The relative PPUE (RPPUE) was derived by dividing the PPUE under low Pi by the PPUE under high Pi.

Presence and/or absence of OsPSTOL1 in rice collections

Genomic DNA from the leaves of 157 rice landraces was extracted using the cetyltrimethylammonium bromide method (Murray and Thompson 1980). A polymerase chain reaction (PCR) was performed to check for the presence of OsPSTOL1 in the Vietnamese rice collection. The sequence of the forward and reverse primers used in this study were 5′-ATGCTGCTCTGTCAAAGGGCAT-3′ and 5′-CAAGCTCAAAGCCCTTTTGGTG-3′, respectively (Mukherjee et al. 2014). The 20 µL reaction mix included 2 µL Dreamtaq buffer (10X) (Dreamtaq bufer; Thermo Fisher Scientific, Waltham, MA, USA), 0.75 units of DreamTaq DNA polymerase (Dreamtaq DNA polymerase; Thermo Fisher Scientific, Waltham, MA, USA), 200 nM deoxyribonucleotide triphosphates (dNTPs; Thermo Fisher Scientific, Waltham, MA, USA), 250 nM forward primer, 250 nM reverse primer, 100 ng DNA template, and H₂O up to 20 µL. The following thermo cycle conditions were used: 95 °C for 4 min, 35 cycles at 95 °C for 30 s, 56 °C for 30 s, and 72 °C for 1 min. The reaction was terminated at 72 °C for 7 min. The PCR products were analyzed by 1% agarose gel electrophoresis at 100 V over 30 min.

Genome-wide association study

The correlation between the four investigated traits and the corresponding SNP loci was studied using the mixed linear model (MLM) implemented in TASSEL software v5.2.55. The structure matrix was identified using principal component analysis with six maximum alleles. The suggestive threshold of p<3.0E-4 was employed to determine significance (To et al. 2019).

Linkage disequilibrium (LD) heatmap and haplotype analysis

LD between SNPs was measured in terms of $r^2$, a statistical measure of correlation between two SNPs. Based on a list of significant markers from TASSEL, an LD heatmap for pair markers was generated in R. $R^2$ > 0.6 was used to consider a region as a QTL in LD blocks. Then, two main haplotypes were compared with their respective phenotypes to assess the effect of sequence variations on trait phenotype.

Screening for candidate genes

To identify the potential candidate genes responsible for Pi starvation tolerance, the positions of significantly associated SNPs, which were clustered into previously identified QTLs, were used. The Michigan State University (MSU) Rice Genome Annotation Project Database Release 7 (https://rice.plantbiology.msu.edu) (Kawahara et al. 2013) was used to browse the expressed genes situated 25 kb upward and downward from each selected significant marker. A list of candidate genes with their annotated functions was generated.

In addition, the expression profiles of candidate genes related to low Pi treatment were also verified based on transcriptome data from a public data source (https://genevestigator.com) (Hruz et al. 2008). Only the genes expressed with changes greater than a definite value of 2 were selected.

PSTOL1-containing varieties

Sixty of the 157 rice accessions (38% of the collection) possessed the OsPsTOL1 gene (Supplementary Table 2). In this study, the correlation between the presence or absence of OsPsTOL1 and Pi-related traits was analyzed and shown in Fig. 3. Three out of four investigated traits, including RRSR, RPPUE, and RPUpE, showed no significant difference between accessions with OsPsTOL1 and those without OsPsTOL1 in response to low Pi treatment. In contrast, significantly higher REP in Pi-deficient soil conditions were observed in varieties not bearing OsPsTOL1 compared to those bearing this gene (Fig. 3a). Furthermore, many accessions, such as G3, G299, G62, and G131 were able to develop well under low-Pi conditions despite the absence of OsPsTOL1 in their genome. This data suggests the existence of another QTL beside Pup1 that can overcome low Pi conditions in the growth stage of our rice collection.

Open in a separate window

Fig. 3

The correlation between with OsPSTOL1 and without OsPSTOL1-bearing varieties in a REP, b RRSR, c RPPUE and d RPUpE traits. The differences of REP, RRSR, RPPUE and RPUpE traits in the two types of rice accessions were analyszd using the Student's t-test. Asterisk (**) indicates the significant difference with p value<0.01

Genome-wide association mapping

A GWAS was performed to identify the association between SNPs and four tested traits in the full rice panel using TASSEL software v.5. The threshold of p<3.0E-4 was set as the significance level (To et al. 2019). The association results of the four traits were graphically displayed in Manhattan and Q-Q plots (Fig. 4). Overall, the Q-Q plot results showed fitted models between the experimental results and theoretical values. A total of 31 SNP-trait associations with p<3.0E-4 were detected. In detail, we identified seven common SNPs for REP, RPPUE, and RPUpE; two SNPs for only REP; two SNPs for only RRSR; 10 SNPs for only RPPUE; and 10 SNPs for only RPUpE. It should be noted that most of the significant SNPs residing on chromosomes 1 and 9 accounted for more than half of the total SNPs found in this study. The remaining SNPs were located on chromosomes 3, 4, 5, 7, 10, and 12. Detailed information regarding the positions and p-values of these significant markers can be found in Supplementary Table 3.

Open in a separate window

Fig. 4

GWAS analysis the effect of Pi starvation on four investigated traits in the rice whole panel. Manhattan plot (left) and Q-Q plot (right) for a REP, b RRSR, c RPPUE and d RPUpE traits. The x-axis indicates the positions of SNPs across 12 chromosomes. The y-axis shows the log scale of p-value for the association test at each locus. The black line represents the significant threshold

QTL identification

Out of the 31 significant SNPs, 18 significant associations, including seven QTLs for RPPUE, two QTLs for RRSR, six QTLs for RPUpE, one common QTL for REP/RPPUE/RPUpE, and two QTLs for REP, were found from our GWAS analysis (Supplementary Table 3). The detected QTLs were distributed across eight chromosomes: 1, 3, 4, 5, 7, 9, 10, and 12. The significant SNPs in each QTL were used to assemble a haplotype sequence, followed by a polymorphism combination analysis was dissected for qRPUUE9.16 (stands for qREP9.16, qRPPUE9.16, and qRPUpE9.16), which was presented in three traits: RPPUE (Fig. 5b), REP (Fig. 5c), and RPUpE (Fig. 5d). Out of the 18 identified QTLs, qRPUUE9.16 was the largest at 57 kb. Moreover, according to Fig. 5a, the LD heatmap showed seven tightly linked markers in the qRPUUE9.16 QTL. A block of these seven closely associated markers indicated a strong linkage between them. Two main haplotypes, TTTTTTT and AAAAAAA, were used to compare the effect of haplotype on traits. The results showed that the accessions containing the AAAAAAA haplotype were significantly more affected than the ones containing the TTTTTTT haplotype under Pi starvation conditions, resulting in lower RPPUE and REP. Four candidate genes located on QTL qRPUUE9.16, including OsWAK86, OsLecRLK, OsGDPD13, and OsFBX318, were found (Fig. 5e).

Open in a separate window

Fig. 5

Haplotype analysis for the QTL qRPUUE9.16 on chromosome 9. a LD heatmaps in the peak region of association analysis GWAS for qRPUUE9.16. The asterisk (*) indicates significant SNPs. “Physical Length” indicates the total length of genetic region. b, c, d Effect of allelic combination of two main haplotypes on RPPUE, REP, RPUpE traits. N indicates the number of accessions for haplotype “AAAAAAA” and “TTTTTTT”. e Candidate genes and their positions on qRPUUE9.16. The level of significant difference between two haplotypes was analysed using the Student's t-test. Asterisk (**) indicates the significant difference with p-value<0.01

Co-locations of the identified QTLs in previous studies

The co-locations of our observed QTLs with previously published QTLs related to Pi starvation were identified. Interestingly, our results showed some similarity with the previous studies, which confirms the reliability of our findings. Specifically, the qRPPUE12.18 for RPPUE located on chromosome 12 shared a similar region with a QTL, found by Wissuwa et al. (1998) for PUP and PUE which was flanked by the two markers G2140 and C443 and qRPPUE1.2 overlapped with a QTL, identified by Wang et al. (2014) for Pi translocation efficiency that was in the proximity of the BIN31-BIN32 markers on chromosome 1. qRPUpE1.5 shared a similar region with a QTL, also found by Wang et al. (2014) for the total above-ground Pi uptake which was flanked by the two BIN46-BIN47 markers on chromosome 1. qRRSR1.4 was close to QTLs, which were detected by Wang et al. (2014) flanked by marker RM283 (Xiang et al. 2015) and markers BIN33-BIN34 (Wang et al. 2014) for relative root dry weight, panicle number per plant, and Pi translocation efficiency.

It is worth noting that the common QTL qRPUUE9.16 was found for three traits, including REP, RPPUE, and RPUpE. Moreover, this QTL was found to be associated with the dry shoot weight, dry root weight, and total dry weight traits (Mai et al., 2020, under revision). These finding indicate that Pi uptake, Pi use efficiency, and REP were strongly correlated to the dry weight trait of rice plants.

Candidate genes for low Pi tolerance in the Vietnamese rice collection

Many mechanisms are activated in response to Pi starvation conditions in plants. Numerous studies have showed that protein kinases and phosphatases play crucial roles in signal transduction in response to stressors at the cellular level in plants (Sopory and Munshi 1998; Luan 1998; Wang et al. 2007). It is interesting to note that among the low Pi-responsive gene list, several genes encoding protein kinases (LOC_Os01g57420.1, LOC_Os01g57480.1, LOC_Os07g42940.1, LOC_Os09g12290.1, LOC_Os09g12300.1, and LOC_Os09g17010.1) and one gene encoding protein phosphatase (PP; LOC_Os10g39540.1) were found.

One very important gene, Os07g0622000 or OsSAPK2, which encodes for stress-activated protein kinase and belongs to the sucrose non-fermenting-1-related protein kinase 2 family (SnRK2s), was identified. This kinase plays a crucial role in energy sensing and in adaptive responses to stress in plants (Colina et al. 2019). OsSAPK2 was found to be located on chromosome 7 near the Dj07_25716591R marker and belongs to qRPPUE7.14. Recently, SnRK2s were reported to be negatively regulated by PP2C via dephosphorylation in the PP2C-SnRK2 complex (Hirayama and Umezawa 2010). Surprisingly, OsSAPK2 has also been shown to be involved in nutrient starvation. SnRK2.8, a homolog of OsSAPK2 in A. thaliana, was downregulated in the roots by Pi deprivation and consequently led to decrease in the growth of A. thaliana under Pi deprivation conditions (Umezawa et al. 2004; Shin et al. 2007). In rice, the expression level of OsSAPK2 was also markedly reduced from 6 to 24 h following the onset of Pi starvation conditions, which resulted from the reduction of the P content in rice grains (Lou et al. 2020). In addition, a number of studies have been carried out to demonstrate the function of this protein in plant defense mechanisms against biotic and abiotic stressors. For example, the expression level of OsSAPK2 was increased in response to Xanthomonas oryzae pv. Oryzae strain PXO99^A (Hu et al. 2015) and Xanthomonas oryzae pv. oryzicola strain RS105. OsSAPK2 was also involved in the abscisic acid signal transduction pathway (Fujita et al. 2011), during drought, high salinity, and polyethylene glycol treatments in rice (Umezawa et al. 2004; Lou et al. 2017; Jadamba et al. 2020). Based on these findings, OsSAPK2 appears to be a promising gene for regulating stress response in rice.

Phosphatases (PPs) are hydrolytic enzymes that cleave a Pi group from its unavailable organic phosphates into assimilable Pi for plant uptake. Since 50–80% of the total Pi in agricultural and natural soils is presented in organic forms (Wang et al. 2009), PPs play an essential role in releasing Pi from these organic sources (Brinch-Pedersen et al. 2002). In plants, the serine/threonine PPs are divided into two types, PP1 and PP2, and PP2 is further divided into three groups, 2A, 2B, and 2C, based on their dependence on divalent cations (Ingebritsen and Cohen, 1983). A number of studies have demonstrated the function of these PPs in signal transduction, metabolism, and hormone regulation and stress response (Luan 1998; Schweighofer et al. 2004; Farkas et al. 2007). The current GWAS identified a PP2C gene, Os10g0541200 (LOC_Os10g39540), located on chromosome 10 near the Dj10_21135126R marker and co-located with qREP10.17. Various studies have been carried out to demonstrate the involvement of this gene in rice plants dealing with salt stress and abscisic acid treatment (Wang et al. 2016) as well as in regulating seed dormancy (Wu et al. 2016). Up till recently, there has not been any research showing the involvement of Os10g0541200 in the Pi starvation response; however, a recent study demonstrated the function of this gene in response to low ammonium nitrate treatment (Yang et al. 2017). Furthermore, a gene encoding protein PP2A in Zea mays L. was shown to be upregulated over 60 times in the roots following a 7 h Pi starvation (Wang et al. 2017). This result showed the role of PP2A in response to Pi starvation. Therefore, Os10g0541200 may be a potential candidate gene for improving plant tolerance to Pi deficiency.

Furthermore, the RNA-Seq transcriptome data was used to search for Pi-responsive genes among our list of candidate genes. In fact, the activation of Pi starvation-induced genes was regulated by a number of TFs, of which PHR1 is the key factor in vascular plants (Devaiah et al. 2007). In our study, although PHR1 was not found, the transcriptional data showed that two other TFs, including bHLH (LOC_Os03g59670.1) and NAC (LOC_Os12g29330.1), were found to be significantly upregulated in the roots and shoots of rice plants exposed to Pi starvation conditions (https://genevestigator.com). In addition, LOC_Os12g29330.1 was observed to be upregulated to response to osmotic stress (Dong et al. 2019) and Magnaporthe oryzae (Zhang et al. 2016). Among the list of candidate genes generated from the transcriptional data that were regulated by low Pi treatment, LOC_Os05g31670.1 named OsPM1 (O. sativa PLASMA MEMBRANE PROTEIN 1) was not only reported to be highly upregulated 214.41 times more in the roots of Nipponbare rice exposed to treated medium without iron compared to medium without both iron and P but also strongly responsive to drought in various tissues obtained from young seedlings to the leaves, roots, stems, and panicles (Yao et al. 2018) as well as to salt stress in the roots (Kong et al. 2019). The differential expression pattern from https://genevestigator.com also revealed another interesting Pi-responsive gene named OsCML15, a calcium signaling-related gene, which was downregulated 15.73 times more in the roots of a Pi-tolerant variety (Dular) compared to a Pi-sensitive one (Pusa Basmati 1) under low Pi supply. Notably, this gene also responded to salt stress by extensively increasing the expression over 100 times more under high salt treatment in rice (O. rufipogon Griff.) (Wang et al. 2017). OsCML15 was also stimulated by brown planthopper (Lv et al. 2014) and Xanthomonas strain PXO99^A (Pérez-Quintero et al. 2013) as well as heat stress (Jung and An 2012). These interesting results corroborate our findings and suggest a new research question for validating the function of these genes in our Vietnamese rice collection.

In addition, glycerophosphoryl diester phosphodiesterase proteins (GDPDs) function in remobilizing glycerophosphate in cells, thus maintaining phosphate homeostasis (Corda et al. 2014). In our study, Os09g0339800, which is located in qRPUUE9.16 close to the Dj09_10404314F marker and named GDPD13, was found. Impressively, a number of studies have demonstrated highly induced expression levels of GDPDs under low Pi supply in Arabidopsis (Misson et al. 2005; Gaude et al. 2008), chickpea (Cicer arietinum L.) (Mehra and Giri 2016), maize (Z. mays L.) (Calderon-Vazquez et al. 2008), and white lupin (Lupinus albus L.) (Cheng et al. 2011) as well as in rice (Mehra and Giri 2016; Mehra et al. 2019). For instance, the expression levels of OsGDPD1 and OsGDPD5 were 42 and 67 times higher in low Pi conditions, respectively, compared to that upon full Pi supply after 7 days in low Pi medium. The expression levels increased further when rice plants were grown in low Pi medium for 15 days. Approximately 65 and 147 times higher expression levels of the OsGDPD1 and OsGDPD5 genes, respectively, were observed when plants were treated with a low Pi supply (Mehra and Giri 2016). GDPD13 was also found in the QTL for root weight, shoot weight, and total weight in our previous study (Mai et al., 2020, under revision). These expression profiles highlight the importance of GDPDs in plant adaptation to Pi starvation conditions and suggest further studies on GDPD13 to understand the function of this gene in acclimatization to Pi starvation.

We would like to thank Prof. Michel Lebrun for the financial support to multiply the rice collection, Dr. Mai Duc Chung, Dr. Phung Thi Phuong Nhung, Ms. Nguyen Van Anh, Ms. Ngo Le Na for their extensive help with experimental procedures and sample harvesting.