1. Introduction
Phosphorus (P) is one of the essential macronutrients in plant growth and development. It is estimated that 43% (about 5.8 billion hm
2) of the world’s arable land is deficient in P, and 3/4 farmlands (about 67 million hm
2) have P shortages in China, which can result in yield reduction by 5–15% (about 25–75 billion kg) [
1]. Although soil available P deficiency can be improved by applying phosphate (Pi) fertilizer, the utilization rate of which plants apply it is no more than 20% [
2]. This is because most of P in soil exists in the form of insoluble mineral P or bound organic P, which cannot be absorbed by plants. In addition, the main source of Pi fertilizer is Pi rock, which is a non-renewable resource and is expected to be depleted soon, and the heavy use of Pi fertilizer can also cause environmental problems such as eutrophication of water [
3].
Plant adaptation to a P-deficiency environment covers a series of gene expression and morphophysiological events [
4], such as regulation of P transporters (PTs), mycorrhizal association, phosphatase secretion, organic acid exudation, and alteration in root structure [
5]. Studies have shown that OsPHR2 (Phosphate Starvation Response 2), homologous to PHR1 in
Arabidopsis, is a major transcriptional regulator of low P response in rice [
6,
7], which could activate the Pi starvation-induced genes including PHT1 (Phosphate Transporter 1) members by binding to the P1BS (PHR1 Binding Sequence; GNATATNC) motif presented in genes’ promoter region [
8,
9,
10,
11,
12,
13].
P-deficiency tolerance, however, is a complex quantitative trait controlled by many genes and is profoundly influenced by the environment [
14]. Quantitative trait locus (QTLs) analysis of P-deficiency tolerance related traits in rice showed that there were generated dozens of QTLs in different populations. These QTLs extensively distributed on chromosomes 1, 2, 3, 4, 6, 7, 9, and 12, especially on chromosomes 4, 6, 11, and 12 [
15,
16,
17]. Under different genetic backgrounds, the QTL loci of some related traits overlapped or were adjacent on the same chromosome, indicating that the traits related to low P tolerance had greater heritability.
Based on the results of QTL mapping or fine mapping, Wasaki et al. [
18] cloned an
OsPI1 gene on rice chromosome 1; Yi et al. [
19] successfully cloned and verified a transcription factor
OsPTF1 that could significantly improve P efficiency in plants; Chin et al. [
15] used the molecular marker closely linked to
Pup1 for assisted breeding, and Gamuyao et al. [
20] successfully cloned the
PSTOL1 gene. Wissuwa et al. [
21] detected 20 P utilization related locus in rice through genome-wide association analysis and identified a candidate gene on chromosome 1 through comparative variation and expression analysis. Meanwhile, some researchers (including our group) constructed interspecific hybrid population with close wild rice and obtained major QTLs for low P tolerance of wild rice, as well as created some new germplasm [
22,
23], which broadened the genetic diversity of P efficient uptake and utilization in rice and laid an important foundation for the utilization of P efficient genes in wild rice. Therefore, excavating the high efficiency P utilization gene of the crop itself will provide insights in solving the yield problem caused by P deficiency and cultivating new varieties resistant to low P.
Dongxiang wild rice (hereinafter referred to as DXWR) is a common wild rice (
O. rufipogon Griff.) found in the northernmost distribution latitude to date. It has more abundant genetic diversity than cultivated rice and contains a large number of excellent genes, including low P tolerance genes, some of which do not even exist in cultivated rice [
23,
24,
25,
26]. Therefore, DXWR is a valuable resource for the excavation and utilization of low P resistant genes. So far, some QTLs related to low P stress tolerance have been identified in DXWR [
23]. In order to understand the molecular mechanism related to low P resistance of DXWR, we detected many important differentially expressed genes associated with P-deficiency tolerance by transcriptome analysis [
26]. However, how DXWR copes with P-deficiency at the protein level is still unclear.
Label-free proteomics analysis is a method that can not only retain the authenticity of the sample to the greatest extent without relying on isotope labeling, but also can compare proteomes affected by different physiological conditions at the same time. Therefore, in this study, we used label-free proteomics analysis to detect the response of DXWR at the protein level under low P stress, and we combined it with the previous transcriptome data to further explore the low P tolerance genes in DXWR. These results would provide insights in explaining the molecular mechanism of low P resistance, as well as cloning and utilizing the P-deficiency genes from wild rice.
3. Results
3.1. Label-Free Quantitative Proteomic Analysis on DXWR
By analyzing the peptide scores obtained by MS, the results showed that about 78.22% of the peptides scored above 60, which meant that the quality of LC-MS/MS experimental data was relatively high (
Supplementary Figure S1). Through label-free proteomic data analysis, a total of 4329 protein groups (
Supplementary Table S2) and 23,598 peptides (
Supplementary Table S3) were identified from six DXWR root samples (RLP and RCK with three biological repetitions, respectively). Among these proteins, we designated 3589 proteins that were detected in at least two replicates as identified proteins (
Figure 1a–c). In addition, the clustering analysis of proteins identified in different samples showed that the similarities between the three biological repetitions of the same treatment were high and between the different treatments were very low (
Figure 1d). Based on this, it was indicated that changes in expression of these target proteins could represent a significant effect of biological treatment on samples. Using these data, we selected proteins with ≥ 1.5-fold changes as an additional standard [
35], and the volcano pot was drawn by using protein expression fold changes and
p value (
Figure 1e). The results showed that 60 protein groups were up-regulated (
Table 1) and 15 were downregulated (
Table 2). Those results indicated that P-deficiency treatment could affect the accumulation of some gene expression products in DXWR.
3.2. Functional Classification by Gene Ontology (GO)
To gain insight into the functional roles of the proteins significantly different between the RCK and RLP samples, GO annotation and enrichment was conducted and the results were listed in
Supplementary Table S4, schematically represented in three ontologies as molecular function, cellular component, and biological process, as in
Figure 2a. The enrichment of biological process involved in metabolic process and cellular process was significantly observed. The most significantly enriched molecular function were catalytic activity and binding. Significant enrichment of cellular compartments was identified, including cell part, cell, membrane, membrane part, and organelle part. From the above description, we could give a conjecture that low P stress could affect cell proliferation and enzyme synthesis as well as the ability of cell or membrane to bind to certain stimulus signals in DXWR.
3.3. Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Mapping
KEGG pathways analysis was performed on the 75 significantly different expression proteins (SDEPs, 60 up-regulated and 15 downregulated) identified in this study (
Supplementary Table S5). These proteins were involved in 31 metabolic pathways. U6 snRNA-associated Sm-like protein 8 (LSm8,
LOC_Os05g51650, up-regulated), U1 small nuclear ribonucleoprotein A (U1A,
LOC_Os05g06280, up-regulated), splicing factor of arginine/serine-rich (SR,
LOC_Os01g06290, downregulated), and pre-mRNA-splicing factor SPF27 (BCAS2,
LOC_Os01g16010, downregulated) were enriched in the spliceosome pathway which was one of the top 20 KEGG pathways predicted to be affected by low P stress (
Figure 2b). Furthermore, branched-chain-amino-acid aminotransferase (BCAT,
LOC_Os03g01600, up-regulated) and acetolactate synthase small subunit (ALS,
LOC_Os11g14950, up-regulated) were enriched in the pathway of branched-chain amino acids (BCAAs, including valine, leucine, and isoleucine) biosynthesis. Ribulose-1,5-bisphosphate carboxylase/oxygenase large subunit (rbcL, encoded by leucoplast, downregulated) and phosphoglycolate phosphatase (PGP,
LOC_Os09g08660, up-regulated) were enriched in the pathway of glyoxylate and dicarboxylate metabolism. A chitinase (
LOC_Os01g49320, up-regulated) and a glycosyl hydrolase (
LOC_Os01g47070, up-regulated) were enriched in amino sugar and nucleotide sugar metabolism. In addition, two GST (
LOC_Os10g38740 and
LOC_Os10g38360, up-regulated) proteins were enriched in the glutathione pathway. These results indicated that genes involved in these pathways might respond to low P stress in DXWR.
3.4. Protein-Protein Interaction (PPI) between the Low-P Responsive Proteins
Protein is an important component of biological organisms, which do not perform biological function independently, but through the interaction of proteins to regulate physiological and biochemical processes. Therefore, we performed a PPI analysis of the proteins identified by label-free quantitative analysis. As shown (
Figure 2c), there was interaction only between a few proteins after screening.
A2X5V0 (Uniprot ID, LOC_Os02g33850, up-regulated) with the function of elongation factor Tu family protein (EF-TU) has five mutual proteins, including Q0D840 (LOC_Os07g08840, up-regulated) with annotation of thioredoxin, A2YMF8 (LOC_Os07g36490, up-regulated) with annotation of RNA recognition, Q5TKF9 (LOC_Os05g51650, up-regulated) with annotation of U6 snRNA-associated Sm-like protein LSm8, B8AHU1 (LOC_Os02g05880, downregulated) with annotation of RNA polymerase, and B8AC69 (LOC_Os01g16010, downregulated) with annotation of BACS2 protein, which is the core component of the Prp19-related complex, along with being involved in important life activities such as splicing of precursor RNA.
Additionally, B8AC69 and A2YMF8 both have four mutual proteins. Additionally, Q0DKM4 (LOC_Os05g06280, up-regulated) with annotation of U1 small nuclear ribonucleoprotein A and A0A0P0XSK6 (LOC_Os10g07229, up-regulated) with annotation of dehydrogenase have three mutual proteins, respectively. Q6EP66 (LOC_Os09g08660, up-regulated) with annotation of phosphoglycolate phosphatase and A6MZ96 (LOC_Os01g06290, downregulated) with annotation of splicing factor have two mutual proteins, respectively. These results suggested that low P stress induced reduction of transcription-related genes, and the increase of most RNA splicing related genes along with intensification of the gene expression associated with the elongation during translation.
3.5. Analysis of Differentially Expressed Proteins Responded to P-Deficiency in DXWR
In the present study, among 75 SDEPs, there were 24 proteins with fold changes larger than 1.5 both in proteomics and transcriptome level [
26] that were associated with P-deficiency treatment in DXWR (
Table 3). In addition, we also verified its expression at the transcriptome level by qRT-PCR, which showed consistency in transcriptome data and qRT-PCR results, shown in
Figure 3. Among these proteins, 21 up-regulated and two downregulated at both transcriptome and proteome level, and only one had the opposite abundance trend (up-regulated in proteome level but downregulated in transcriptome level). There were some genes whose expression abundance increased in both protein and transcription levels, including
OsPT2 (
LOC_Os03g05640),
OsPT8 (
LOC_Os10g30790),
OsPAP10c (
LOC_Os12g44020),
OsPAP10a (
LOC_Os01g56880),
OsPHF1 (
LOC_Os07g09000), as well as a gene encoding glycerophosphoryl diester phosphodiesterase family protein (GDPD,
LOC_Os03g40670), three genes encoding glycosyl hydrolase (
LOC_Os07g23850,
LOC_Os05g15770 and
LOC_Os01g47070), and one gene encoding glucan endo-1,3-beta-glucosidase precursor (
LOC_Os07g35560). Furthermore, a gene (
LOC_Os10g35500) encoding epoxide hydrolase increased by low P stress at proteomic level, but the corresponding transcription levels decreased. The above results indicated that after low P stress treatment, there would be differences in the trend and degree of change in gene expression at the transcription and translation levels.
3.6. Conjoint Analysis of Proteomic and QTLs Related to P-Deficiency Tolerance
As shown in
Table 4, there are two P-deficiency tolerance related QTLs have been identified in DXWR [
23], and 12 QTLs existing in different positions on the chromosome related to P-deficiency stress have been found in
Oryza sativa based on the Gramene QTL database. Among genes corresponding to 75 SDEPs identified by the proteome in this study, we located nine genes among these QTL intervals, as shown in
Table 5. Among them, two genes (
LOC_Os12g44020,
OsPAP10c and
LOC_Os04g41970,
OsGLU3) have been characterized in previous studies [
36,
37], and two of the other seven uncharacterized genes (
LOC_Os12g09620 and
LOC_Os03g40670) have been detected at both transcriptome and proteome levels. Furthermore, the functional expression characteristics of the remaining five genes (
LOC_Os01g57450,
LOC_Os03g29240,
LOC_Os03g13540,
LOC_Os03g29190 and
LOC_Os06g07600) have not been reported yet.
3.7. The Expression Pattern of Genes Related to P-Deficiency Tolerance in DXWR
Based on previous studies on P-response mechanism in cultivated rice, the homologous genes of
PHR1 (
OsPHR2 and
OsPHR1),
OsPHO2,
OsPHO1, as well as its NATs, play an important role in the low P-response process [
7,
38]. Therefore, we examined the expression levels of
OsPHR2 and
OsPHO2 as well as three members of
OsPHO1 (
OsPHO1;1,
OsPHO1;2 and
OsPHO1;3) and three NATs correspondently, as shown in
Figure 4. We found that
OsPHR2 was downregulated after low P treatment in DXWR roots, which was contrary to the results of previous studies on cultivated rice [
7,
38]. We suspect that
OsPHR1 may play a major role during the P starvation signaling pathway in DXWR. For validation, we examined the expression of
OsPHR1 in DXWR, but the result was consistent with
OsPHR2 and also down-regulated. Such results may indicate that the transcription levels of
OsPHR1 and
OsPHR2 are inhibited by low P signals in DXWR. In addition,
OsPHO2 was down-regulated in both DXWR and NP. Among
OsPHO1;1,
OsPHO1;2, and
OsPHO1;3, only
OsPHO1;3 expression was up-regulated and the other two were downregulated in DXWR, but three genes all up-regulated in NP. In the corresponding NAT, only
OsPHO1;2 NAT expression change trend is consistent with
OsPHO1;2, and the other two NATs change trend is opposite to
OsPHO1;1 and
OsPHO1;3 in DXWR. However, all three NATs up-regulated in NP. These results indicate that there may be a difference in the mechanism of resistance to low P in DXWR compared to cultivated rice.
5. Conclusions
In this study, label-free proteomics analysis as well as joint analysis with transcriptome dataset were conducted to root samples to identify potential unique low P-response genes in DXWR during seedlings. 75 SDEPs were detected, 24 of which were also detected in previous transcriptome dataset verified by qRT-PCR. Furthermore, it was found that DXWR could increase PAPs’ expression, membrane location of PTs, rhizosphere area, alternative splicing and decrease ROS activity to deal with low P stress. Moreover, among the genes corresponding to 75 SDEPs, seven uncharacterized genes were located in previous P related QTL intervals, of which two genes (LOC_Os12g09620 and LOC_Os03g40670) have been detected at both transcriptome and proteome levels. In addition, the expression patterns of OsPHR1, OsPHR2, OsPHO1, and NAT-OsPHO1 in DXWR were different in cultivated rice NP, suggesting that the response mechanism of some low P tolerance in DXWR might be different from that in cultivated rice. These findings would provide insights in cloning the P-deficiency genes from wild rice, as well as elucidating the molecular mechanism of low P resistance in DXWR.