Identification of Phosphorus Stress Related Proteins in the Seedlings of Dongxiang Wild Rice (Oryza Rufipogon Griff.) Using Label-Free Quantitative Proteomic Analysis

Phosphorus (P) deficiency tolerance in rice is a complex character controlled by polygenes. Through proteomics analysis, we could find more low P tolerance related proteins in unique P-deficiency tolerance germplasm Dongxiang wild rice (Oryza Rufipogon, DXWR), which will provide the basis for the research of its regulation mechanism. In this study, a proteomic approach as well as joint analysis with transcriptome data were conducted to identify potential unique low P response genes in DXWR during seedlings. The results showed that 3589 significant differential accumulation proteins were identified between the low P and the normal P treated root samples of DXWR. The degree of change was more than 1.5 times, including 60 up-regulated and 15 downregulated proteins, 24 of which also detected expression changes of more than 1.5-fold in the transcriptome data. Through quantitative trait locus (QTLs) matching analysis, seven genes corresponding to the significantly different expression proteins identified in this study were found to be uncharacterized and distributed in the QTLs interval related to low P tolerance, two of which (LOC_Os12g09620 and LOC_Os03g40670) were detected at both transcriptome and proteome levels. Based on the comprehensive analysis, it was found that DXWR could increase the expression of purple acid phosphatases (PAPs), membrane location of P transporters (PTs), rhizosphere area, and alternative splicing, and it could decrease reactive oxygen species (ROS) activity to deal with low P stress. This study would provide some useful insights in cloning the P-deficiency tolerance genes from wild rice, as well as elucidating the molecular mechanism of low P resistance in DXWR.


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.

Plant Materials and Phosphorus Deficiency Treatment
In the present study, DXWR from Jiangxi academy of agricultural sciences was carried out as experimental material. The DXWR seeds were surface sterilized using mixed solutions of NaClO (10%) for 15min and soaked in petri dishes containing 20 mL deionized water at room temperature for 3 days. Then, the seeds were selected with the same growth trend and planted in the plastic pot containing quartz sand in a climate control chamber at day/night 14 h/10 h (30 • C/26 • C) [26]. A 1/2 Yoshida culture medium (pH 5.8) was added once for growth when germinated seeds had coleoptiles 10 mm approximately in length [27]. At the emergence of the third leaf (about 15 days), plants were transferred into either a culture medium with low P concentration (0.016 mM NaH 2 PO 4 ) or normal P concentration (0.32 mM NaH 2 PO 4 ), which corresponded to the −P treatments (RLP) and +P treatments (RCK), respectively. The culture medium was replaced every three days. There were 20 seedlings per treatment with three biological replications. Roots were harvested at 9 days after the experimental treatments started, and the samples were frozen immediately using liquid nitrogen and stored at −80 • C for further analyzing. Root samples of cultivated rice Nipponbare (NP) were obtained in the same way at the same time and used to determine the expression of OsPHR2, OsPHR1, OsPHO2 (Phosphate2), and OsPHO1 (Phosphate Transporter 1), as well as its natural reverse transcripts (Cis-Natural Antisense Transcripts, NATs).

Protein Extraction and Enzymatic Hydrolysis
Root proteins of six DXWR samples were extracted using trichloroacetic acid (TCA)acetone precipitation method [28]. Briefly, the sample was ground into powder in liquid nitrogen and then suspended in extraction buffer (8 M urea, 1% DTT, 0.1 M Tris-HCl, pH 8.8, 1% complete protease inhibitors (Roche, Mannheim, Germany)). Repeated vortex of the sample and then removal of insoluble precipitation was performed by centrifugation at 14,000× g for 40 min. The supernatant was precipitated overnight with 20% (v/v) TCA, washed three times with cold acetone, and solubilized in extraction buffer. All operations were performed at low temperatures. The protein concentration was measured using the BCA Protein Assay Kit (Bio-Rad, Hercules, CA, USA). The final concentration of urea in protein solution was adjusted to 2 M with 40 mM NH 4 HCO 3 solution. An amount of 4 µg trypsin (Promega, Madison, WI, USA) was added to each sample containing 200 µg protein and incubated overnight at 37 • C following the instructions of the manufacturer.

Liquid Chromatography and Tandem Mass Spectrometry Proteomics Analysis (LC-MS/MS) of the RLK and RCK Samples
LC-MS/MS analysis was performed on a Q-Exactive mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA) that was coupled to an Easy nLC Biosystem (Thermo Fisher Scientific, Waltham, MA, USA) with the help provided by Shanghai Applied Protein Technology (Shanghai, China). Balance chromatographic column with buffer A (0.1% formic acid, 3% acetonitrile and 97% H 2 O). Each sample was automatically injected into the prepacked column (2 cm × 100 µm 3 µm-C 18 ), then flowed into an analytical column (10 cm × 75 µm 3 µm-C 18 ) at a speed of 250 nL/min controlled by intelliflow technology. After that, the sample was separated with a linear gradient of buffer B from 6% (80% acetonitrile, 0.08% formic acid) to 95% over 116 min and then followed by an equilibration of the column at 6 % buffer B for 4 min. MS/MS spectra were searched using MaxQuant software (version 1.5.3.17) against the UniProt proteome database (Uni-port_Oryza sativa_168264_20171201.fasta), and the label-free quantitation algorithm was performed for quantitative analysis. The maximum missed cleavages used for the database search were set to 2. The mass tolerance was set to 20 ppm on full scans. For label-free quantitative methods, retention time matching between runs was performed within a time window of 2 min. The peptide false discovery rate (FDR) and protein FDR did not exceed 0.01. This rigorous analysis tool named Andromeda was used for analysis-obtained excellent peptide score distribution to judge the quality of MS experimental data [29]. Quantifiable proteins were defined as those identified at least twice in the three biological replicates. Proteins with an adjusted p value < 0.05 were assigned as differentially expressed between the RLK and RCK samples [30].

Bioinformatics Analysis
With the help provided by Shanghai Applied Protein Technology (Shanghai, China), bioinformatics analysis was performed on the obtained proteome data. First, we searched the EBI database for conserved motifs that matched the target protein through Inter-ProScan [31] and annotated the motif-related functional information to the target protein sequence to achieve gene ontology (GO) functional annotation [32,33]. Then, proteins were matched to the Kyoto encyclopedia of genes and genomes (KEGG) database to obtain the pathway they might participate in. Last but not least, proteome data set was used to construct protein-protein interaction (PPI) network by using STRING (https://cn.string-db.org/), and the parameter was set to moderate confidence (0.400).

Gene Expression Analysis by Quantitative Real-Time PCR (qRT-PCR)
For qRT-PCR analysis, total RNA in roots samples were isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. The first cDNA was synthesized with 2 µg of total RNA, using ReverTra Ace ® qPCR RT Master Mix with gDNA Remover (TOYOBO, Osaka, Japan). Briefly, RNA was pre-denatured at 65 • C, and gDNA was removed by adding DN Master Mix with gDNA Remover at 37 • C for 5 min. The first strand cDNA was synthesized by using RT Master Mix through the following three-step reaction: 37 • C for 15 min, 50 • C for 5 min, and 98 • C for 5 min. Then, synthesized cDNAs were used as templates for qRT-PCR with THUNDERBIRD™ SYBR ® qPCR Mix (TOYOBO, Osaka, Japan) and the LightCycler ® 96 instrument (Roche, Mannheim, Germany) according to the manufacturer's manuals. The conditions of the qRT-PCR reaction are set as follows: 95 • C for 30 s with 1 cycle for pre-denaturation; two-step reaction with 40 cycles of 10 s at 95 • C and 30 s at 60 • C for amplification; three-step reaction for melt curve stage with 95 • C for 10 s, 65 • C for 60 s, and 97 • C for 1 s; 37 • C for 30 s for cooling. OsActin1 (LOC_Os03g50885) was used as an internal control [34]. The 2 −∆∆Ct method was used for relative quantification. The statistical significance was evaluated by t-test analysis. The primers used are listed in Supplementary Table S1.

Conjoint Analysis of Proteomic and Transcriptomics Related to Low P Stress
We selected the DXWR root transcriptome data obtained from samples treated at the same time with this study for comparison [26], and we screened the transcriptome data with the gene expression change fold ≥ 1.5 times and p value ≤ 0.05 after low P treatment for analysis.

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.  Table 2). Those results indicated that P-deficiency treatment could affect the accu tion of some gene expression products in DXWR.   The protein expression fold change (X-axis) was plotted against the p value obtained from t test log10-value (Y-axis). Small circle represents protein. The red circle represents a protein with a change fold greater than 1.5.

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.

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

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.

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-betaglucosidase 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. Table 3. Significant differential expression proteins (p ≤ 0.05) with fold changes both in transcriptome and proteomics level larger than 1.5.

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.

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 downregulated 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.  Figure 3. Quantitative real-time PCR analysis of 24 significant differential expression proteins with fold changes both transcriptome and proteomics levels larger than 1.5 in DXWR. Bars mean SD. Expression change fold refers to the chan of the treatment group compared with the control group.

Conjoint Analysis of Proteomic and QTLs Related to P-deficiency Tolerance
As shown in Table 4, there are two P-deficiency tolerance related QTLs hav identified in DXWR [23], and 12 QTLs existing in different positions on the chrom related to P-deficiency stress have been found in Oryza sativa based on the Gramen database. Among genes corresponding to 75 SDEPs identified by the proteome study, we located nine genes among these QTL intervals, as shown in Table 5. A them, two genes (LOC_Os12g44020, OsPAP10c and LOC_Os04g41970, OsGLU3) hav characterized in previous studies [36,37], and two of the other seven uncharac genes (LOC_Os12g09620 and LOC_Os03g40670) have been detected at both transcri  These results indicate that there may be a difference in the mechanism of resistance to low P in DXWR compared to cultivated rice.

Low P Stress Leads to Differential Expression of P Absorption Efficiency Related Genes in DXWR
A previous study [39] has shown that OsPHF1 regulates the plasma membrane localization of low-affinity Pi transporter OsPT2 and the high-affinity Pi transporter OsPT8, of which ortholog in Arabidopsis reported to be only an important factor for the localization of high-affinity Pi transporters to the plasma membrane [40]. Subcellular location experiments show that mutation of OsPHF1 lead to the retention of OsPT2 and OsPT8 in the endoplasmic reticulum and reduce the accumulation of Pi in shoots due to overexpression of OsPHR2 [39,41]. To the contrary, overexpression of OsPHF1 results in excessive Pi accumulation in leaf and root. Furthermore, OsPT2 is the only low affinity transporter in the PHT1 family induced by Pi deprivation under the transcriptional control of OsPHR2, whereas OsPT8 is constitutively expressed high-affinity Pi transporters in rice whose expression is not affected by external Pi levels. In our study, the expression levels of OsPT2 (LOC_Os03g05640), OsPT8 (LOC_Os10g30790) and OsPHF1 (LOC_Os07g09000) were upregulated. The up-regulated expression of OsPHF1 may increase the plasma membrane localization of OsPT2 and OsPT8.
On the other hand, purple acid phosphatase (PAP) is a family of metals phosphoesterases involved in a variety of physiological functions, especially in low Pi adaptations in plants [42]. PAPs have non-specific acidic phosphatase activity, which can catalyze the hydrolysis of various organic P into Pi under acidic pH conditions, thus providing more Pi for plants [36,43]. PAP plays a critical role in the plant's ability to utilize organic P in growth medium. There were two PAP genes, OsPAP10c (LOC_Os12g44020) and OsPAP10a (LOC_Os01g56880) induced by low P stress in DXWR. It is likely that they play a crucial role in the ability of plants to use organic P.
In addition to the above-mentioned genes that may increase the efficiency of P absorption of DXWR, there are other functional genes that are differentially expressed. Studies have shown that PGP inactivation attenuated triosephosphate isomerase activity, thereby increasing triglyceride levels at the expense of the cellular phosphatidylcholine content, and inhibiting cell proliferation [44,45]. These effects were prevented under hypoxic conditions or by blocking phosphoglycolate release from damaged DNA. Moreover, as shown, EF-TU plays an important role in the reproduction, development, and response to environmental stress of higher plants [46]. Here, the increased synthesis of PGP (LOC_Os09g08660) and EF-TU (LOC_Os02g33850) may be the response of DXWR to low P stress. Furthermore, chitinase leads to the separation of parallel chitin microfibrils connected by β-1,6-branched chain β-1,3glucans in the cell wall, thus increasing the interval between the insertion of newly synthesized chitin and β-1,3glucans under swelling in vivo [47]. In this study, it was worth noting that the expression of chitinase (LOC_Os01g49320) and glycosyl hydrolase (LOC_Os01g47070) increased when the root of DXWR was under low P stress. The differential expression of these genes may imply a strategy for DXWR to respond to low P stress by improving P absorption efficiency.

Differential Expression of Variable Splicing-Related Genes May Contribute to the Low P Resistance for DXWR
Studies have shown that the absence of SR or SPF protein can lead to changes in splice sites [48,49]. LSm8 is essential for the assembly of the LSM nuclear complex (LSm2-8) and this complex acts in pre-mRNA splicing through U6 snRNA stabilization, thus allowing the formation of the U6 snRNP [50]. The Arabidopsis LSM2-8 complex differentially regulates plant tolerance to abiotic stresses by controlling the constitutive and alternative splicing of specific introns from selected abiotic stress-related pre-mRNAs [51]. In this study, the up-regulation of LSm8 (LOC_Os05g51650) and U1A (LOC_Os05g06280) may alternate splicing of pre-mRNA, while the down-regulation of SPF (LOC_Os01g16010) and SR (LOC_Os01g06290) may change the pre-mRNA splicing site, thus contributing to low P resistance of DXWR.

Increase the Antioxidant Capacity of DXWR by Regulating the Expression of Related Genes
Studies have shown that plants under abiotic stress are coerced to increase the activity of reactive oxygen species (ROS) and antioxidant enzymes [52,53]. As shown, a key enzyme in the BCAAs biosynthesis pathway, ALS, was downregulated under stress [52]. In addition, by increasing the concentration of BCAAs to 100 mg/L in the culture medium, ROS was significantly reduced, thereby reducing the level of antioxidant enzymes in herbicide-stressed plants [53]. In present research, DXWR may enhance the biosynthesis of BCAAs through up-regulating BCAT (LOC_Os03g01600) and ALS (LOC_Os11g14950) to reduce ROS and consequently antioxidant enzyme levels. What s more, Glutathione is a co-substrate for glutathione-S-transferase (GST), which in rice participates in various functions such as phytohormone homeostasis, hydroxy peroxide detoxification, apoptosis regulation [54], and also has a key role in response to biotic and abiotic stresses [55]. In this study, two GST proteins (LOC_Os10g38740 and LOC_Os10g38360) enriched in glutathione metabolism pathway were up-regulated after P-deficiency treatment of DXWR which might increase the tolerance to cope with low P stress. Furthermore, epoxide hydrolase has been reported as an enzyme that reduces the content of epoxides in organisms by means of chemically catalyzed transformation, as well as metabolizes endogenous aliphatic and aromatic epoxides [56]. It is a key enzyme involved in metabolism, detoxification, and signaling regulation in the organism [57]. In this study, an epoxide hydrolase protein (LOC_Os10g35500) increased in protein level might increase the tolerance to low P stress, but the corresponding transcription level decreased, possibly by increasing the stability of the corresponding mRNA and reducing the number of transcriptions to reduce energy consumption.

DXWR May Exist a Low P Tolerance Mechanism Different from Cultivated Rice NP
PHR1 is a central regulator of P-deficiency stress response [6]. There are two homologous genes of PHR1 in rice, namely OsPHR1 and OsPHR2. Some studies have shown that OsPHR2 plays a major role in the P starvation signaling pathway that is induced by low phosphorus stress. [7,38]. However, in this study, both OsPHR2 and OsPHR1 were downregulated after receiving low P stress in DXWR, which was opposite to that in NP. After Pi enters the column cells of the root, PHO1 located in the xylem of the root vascular bundle is responsible for loading Pi into the xylem and then transporting it from the xylem to the shoot [58,59]. In Arabidopsis, mainly AtPHO1 and AtPHO1;H1 are involved in the P transportation from root to shoot [60]. There are three homologous genes of AtPHO1 and AtPHO1;H1 in rice, which are OsPHO1;1, OsPHO1;2, and OsPHO1; 3, and all three genes have NAT at the 5' end [38]. Studies among cultivated rice showed that OsPHO1;2 played a main role in the transport and distribution of P, while OsPHO1;2 NAT was induced by P starvation and could activate the expression of OsPHO1;2 [38], which was consistent with our results in NP. However, in DXWR, OsPHO1;2 and its NAT both downregulated by low P stress, whereas OsPHO1;1 downregulated and its NAT up-regulated, as well as OsPHO1;3 up-regulated and its NAT downregulated. In addition, OsPHO2 containing a E2 ubiquitin-binding domain could degrade OsPHO1 [61], of which the downregulated expression results consistent with previous studies were obtained in both DXWR and NP. Through the above analysis, we found that OsPHR1 OsPHR2, three OsPHO1, and the corresponding NATs showed different response trends in DXWR and NP, except for OsPHO2. These differences indicate that DXWR may have a unique resistance to low P regulation, and these candidate genes screened by transcriptome and proteome in the present study may also have unique functions in DXWR different from cultivated rice.

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.

Conflicts of Interest:
The authors declare no conflict of interest.