Transcriptome analysis of phosphorus stress responsiveness in the seedlings of Dongxiang wild rice (Oryza rufipogon Griff.)

Low phosphorus availability is a major factor restricting rice growth. Dongxiang wild rice (Oryza rufipogon Griff.) has many useful genes lacking in cultivated rice, including stress resistance to phosphorus deficiency, cold, salt and drought, which is considered to be a precious germplasm resource for rice breeding. However, the molecular mechanism of regulation of phosphorus deficiency tolerance is not clear. In this study, cDNA libraries were constructed from the leaf and root tissues of phosphorus stressed and untreated Dongxiang wild rice seedlings, and transcriptome sequencing was performed with the goal of elucidating the molecular mechanisms involved in phosphorus stress response. The results indicated that 1184 transcripts were differentially expressed in the leaves (323 up-regulated and 861 down-regulated) and 986 transcripts were differentially expressed in the roots (756 up-regulated and 230 down-regulated). 43 genes were up-regulated both in leaves and roots, 38 genes were up-regulated in roots but down-regulated in leaves, and only 2 genes were down-regulated in roots but up-regulated in leaves. Among these differentially expressed genes, the detection of many transcription factors and functional genes demonstrated that multiple regulatory pathways were involved in phosphorus deficiency tolerance. Meanwhile, the differentially expressed genes were also annotated with gene ontology terms and key pathways via functional classification and Kyoto Encyclopedia of Gene and Genomes pathway mapping, respectively. A set of the most important candidate genes was then identified by combining the differentially expressed genes found in the present study with previously identified phosphorus deficiency tolerance quantitative trait loci. The present work provides abundant genomic information for functional dissection of the phosphorus deficiency resistance of Dongxiang wild rice, which will be help to understand the biological regulatory mechanisms of phosphorus deficiency tolerance in Dongxiang wild rice.


Background
Phosphorus (P) is an essential mineral element required for plant growth and development. It plays a crucial role in the processes of energy transfer, signal transduction, photosynthesis and respiration [1]. Rice is one of the most important food crops in the world, and it is also considered as a model organism for monocotyledon genomics research [2]. Although the total amount of P in the soil may be high, it is often present in unavailable forms or in forms that are only available outside the rhizosphere [3]. Thus, the deficiency of phosphorus in soil is a worldwide problem. According to the statistics, about 5.8 billion hm 2 of arable land worldwide is deficient in P. And about 67 million hm 2 of arable land in China is deficient in P, which resulted in yield reduction by 5-15% (about 25-75 billion Kg) [4]. Over the past several decades, P fertilizer application for crop growth has been increased rapidly, but P-use efficiency has decreased to a low level of 10-20% [5]. P fertilizer consumes non-renewable phosphate rock reserves, which are expected to be exhausted

Open Access
Biological Research *Correspondence: xdluolf@163.com; xiejiankun@yahoo.com 1 College of Life Science, Jiangxi Normal University, Nanchang 330022, China Full list of author information is available at the end of the article in the near future [6]. Furthermore, much of the applied P has caused serious environmental pollutions [7]. Therefore, it is a significant development direction to improve the status of P-deficiency in crops to explore the absorption and utilization of high P-use efficiency of crops.
Previous studies have found that there are significant genotypic differences in P-deficiency tolerance among different genotypes [8] which make it possible to breed varieties with improved P-deficiency tolerance. During the course of domestication from wild rice to cultivated rice, however, the number of alleles of cultivated rice reduced by 50-60% compared to wild rice [9]. Moreover, the wide application of commercial hybrids in recent decades has caused the loss of many excellent local varieties, and the genetic diversity of cultivated rice has become more and more narrow.
Dongxiang wild rice (O. rufipogon Griff., Hereinafter referred to as DXWR), a Chinese type of wild rice grown in Dongxiang County, Jiangxi Province (28°14′ N latitude and 116°30′ E longitude), is considered to be the northernmost region in the world where O. rufipogon is found, which is one of unique wild resources in Jiangxi province and one of the second class national protected wild plants in China [10]. Previous studies confirmed that DXWR has many useful genes lacking in cultivated rice, including P-deficiency tolerance [11], strong cold tolerance [12], high grain yield [13] and drought resistance [14]. Our previous results and relative reports indicated that DXWR was more resistant to low-P stress than low-P tolerant cultivated rice 'Dalidao' and 'Liantangzao3' , suggesting that DXWR has strong resistance to low-phosphorus stress [15]. Therefore, it is considered to be a valuable resource for the exploitation and utilization of P-deficiency tolerance genes in rice. So far, there has been a great deal of quantitative trait loci (QTLs) for abiotic stress tolerance in DXWR, including P-deficiency tolerance [14,[16][17][18][19][20]. However, the regulatory mechanisms of P-deficiency tolerance have not been fully understood. Thus, a better understanding of P-deficiency tolerance mechanisms would be helpful for breeding P-deficiency tolerance rice cultivars.
In this study, the transcriptome of DXWR under P-deficiency stress was obtained by experiment. Then some candidate genes were found by combining the DEGs interval of this study with the previously identified QTLs interval associated with P-deficiency tolerance. The results will provide a basis for explaining the molecular mechanism of resistance, as well as cloning and utilizing the P-deficiency genes from wild rice.

Material and treatment
Dongxiang wild rice (O. rufipogon Griff., hereafter referred to as DXWR) was used as test material. Disinfection of rice seeds with 10% sodium hypochlorite, then soak the seeds with water at room temperature for 30 h. Selection of the uniformly germinated seeds of DXWR were grown in a plastic pot in a plant growth chamber at day/night temperature of 30 °C/26 °C (14 h day/10 h night) with relative humidity of 70% and 3000 lx of light intensity. Germinated seeds with coleoptiles 8-12 mm in length, adding Yoshida culture medium [21]. At two and half leaves stage (about 15 days), the seedlings were treated with P-deficiency. P-deficiency (0.016 mM NaH 2 PO 4 ) and P-sufficiency (0.32 mM NaH 2 PO 4 ) were as treatment and control, respectively. There were 10 seedlings per treatment with three replications. The growth culture solution was renewed every 3 days. After 9 days treatment under P-deficiency stress, leaves and roots under low P-stressed tissues (LLP and RLP) and control group stressed tissues (LCK and RCK) were collected and immediately frozen in liquid nitrogen, then stored at − 80 °C until RNA extraction.

RNA extraction, cDNA library preparation, and transcriptome sequencing
Total RNA was extracted for three biological replicates from the sampled leaf or root tissues using the TRIzol kit following the instructions (Invitrogen) of manufacturer. The quality and quantity of the resulting RNA were examined using agarose gel electrophoresis and an ND-1000 spectrophotometer (NanoDrop Technology, USA). Magnetic beads with Oligo (dT) were used to isolate mRNA from the total RNA. Mixed with the fragmentation buffer, the mRNA is fragmented into short fragments. Then cDNA was synthesized using the mRNA fragments as templates. Short fragments were purified and resolved with EB buffer for end reparation and poly (A) addition. After that, the short fragments were connected with adapters. After agarose gel electrophoresis, the suitable fragments (200 bp) were selected for the PCR amplification as templates. The library was sequenced using the Illumina HiSeq ™ 2000 platform.

Reads filtration and assessment of differential gene expression
Before assembly, adaptor sequences, empty reads, low quality sequences with 'N' percentage over 10% and those containing more than 50% bases with a Q < 20 were removed using the Perl program written according to the custom method of program editing. After filtering, the remaining reads were called clean reads and used for downstream bioinformatics analysis. The retained highquality reads were mapped to the Nipponbare reference genome [22] by TopHat. And then the resulting aligned reads were used to create a RABT (Reference Annotation Based Transcript) assembly using Cufflinks [23].
Expression levels for each gene were calculated by quantifying the reads according to the RPKM (reads per kilobase per million reads) method [24]. We used 'FDR (false discovery rate) ≤ 0.001 and the absolute value of log 2 RPKM ratio ≥ 1' as the threshold to judge the significance of gene expression difference [25].

Gene ontology (GO) term analysis
Blast2GO program was used to classify unigenes to GO terms based on molecular function, biological processes and cellular components [26] for leaf and root tissues, at p < 0.05.

Validation of transcriptome sequencing
qRT-PCR was performed to confirm 15 randomly selected differentially expressed genes (DEGs) (10 upregulated and 5 down-regulated genes) from each group of P-deficiency-induced genes identified from RNA sequencing using the SYBR premix Ex Taq kit on a Ste-pOnePlus ™ Real-Time PCR System. Diluted cDNA was amplified using gene specific primers (Additional file 1: Table S1) and SYBR Green real-time PCR master mix (Toyobo). All reactions were performed using one biological sample and three technical replicates, and each sample was conducted at least in triplicate and normalized using OsActin1 as an internal control [27].

Transcriptome sequencing statistics
To understand the molecular mechanism of DXWR under P-deficiency stress, we investigated the gene expression of DXWR in response to P-deficiency stress by transcriptome sequencing. The total RNA was extracted from the leaves of DXWR at seedling stage to construct cDNA library for transcriptome sequencing analysis. In total, 46.78, 48.82, 43.59, and 46.76 million high-quality paired end reads were generated by Illumina-sequencing the LCK, LLP, RCK, and RLP cDNA libraries, respectively ( Table 1). The O. sativa ssp. japonica cv. Nipponbare genome has been completely sequenced through Sanger sequencing technology which is considered to be the best tool for assembling and annotating the rice genomes [28][29][30]. Therefore, in this study, we used the Nipponbare genome as reference for reads matching. The alignment results indicated that 73.37-75.20% (71.43-73.72% uniquely matched) of the total reads were mapped to the Nipponbare reference genome and 61.35-63.52% (38.08-39.57% uniquely matched) were mapped to the gene regions (Tables 1, 2). Meanwhile, there was a significant difference in the percentage of reads matched to the genome and the gene region, especially the unique matching reads number, which was similar to previous findings. The previous study reported that about 72% of the total reads mapped to the genome and 46% to the gene regions with 68 and 38% uniquely matched for deep transcriptome sequencing of rhizome and aerial-shoot in Sorghum propinquum [31]. This may be due to the fact that reads match the intergenic spacer or alternatively spliced region of mRNA. The 24.80-26.63% of reads remained unmapped, mainly attributable to gene intervals and the differences between DXWR and reference genome sequence (Table 1). Among more than 60% of the mapped genes, on average, at least 50% were covered by the uniquely mapped reads, and only approximately 15% of the genes had gene coverage of 20% or lower (Additional file 2: Figure S1), which suggested that a high quality of transcriptome data was obtained.
Previous studies have shown that Asian cultivated rice was domesticated from wild varieties (O. rufipogon) [32]. Moreover, the sequence of wild rice W1943 has a great similarity with that of Nipponbare [33]. However, some of the sequences in the W1943 cDNAs that could not be matched to the genome may be located in the gap of genomic sequences or may be related to the W1943 specific gene [33], which suggests that the use of the Nipponbare genome to match wild rice (O. rufipogon) is limited.

Analysis of differentially expressed genes
The raw data obtained from Illumina sequencing can be used to assess the level of gene expression [34]. Putative DEGs from the P-deficiency-stressed and control samples (LLP vs. LCK and RLP vs. RCK) were identified. In the LLP sample, 323 and 861 transcripts were up-(Additional file 3: Table S2) and down-regulated (Additional file 4: Table S3), respectively, when compared to the LCK sample. In the RLP sample, 756 and 230 transcripts were up-(Additional file 5: Table S4) and down-regulated (Additional file 6: Table S5), respectively, when compared to the RCK sample. Among these DEGs, 43 and 35 transcripts were up-(Additional file 7: Table S6) and downregulated (Additional file 8: Furthermore, many of the detected DEGs represented genes that have been previously identified and demonstrated to play roles in responses to P-deficiency stress in cultivated rice (Table 3).
In this study, we found that some genes have been shown to be associated with P-deficiency tolerance, such as MYB family TF (LOC_Os01g16810.1), WRKY TFs (LOC_ Os01g60640.1, LOC_Os01g61080.1, LOC_Os01g53040.1),  [45]. PHR1, a MYB TF, involved in response to P-deficiency stress [46], which was not very sensitive to P-deficiency stress, and acted a role in the downstream of the phosphorus signal transduction pathway. PHR1 consists of a MYB domain and a coiled coil domain (perhaps to form a dimer with an imperfect palindromic sequence at specific promoters). It plays a role in maintaining the balance of P under the condition of adequate nutrition. It was also found that the increase of sulfur transporter and iron transporter was regulated by P stress [45]. This suggested that the transcription factors played an important role in the transcriptional regulation of downstream genes in plants under P-deficiency stress. Previous studies have found some zinc finger proteins responded to P-deficiency stress, for example, two C2H2-type zinc finger protein gene ZOS3-12 (LOC_ Os03g32230) and ZOS5-08 (LOC_Os05g37190): the former was proved to be related to nitrogen stress in rice, and the latter was related to defoliation, which suggested that zinc finger proteins played an important role in the regulation of multiple stress tolerance.
Uptake, transport and translocation of P in plants are performed by P transporters (Pht). Under P-deficiency stress, lipid changed from phospholipid to galactose and sulfanilamide, phosphatase activity increased, and phosphate monoester decomposed into Pi and the related fatty acid and phosphate transporter gene will be induced [47]. In this study, we also found some transporter genes that have been shown to be associated with P-deficiency tolerance. In the LLP sample, we found three genes editing transmembrane amino acid transporter (LOC_ Os01g41420. 1  Pht1;1 ~ Pht1;4 were related to plant uptake of Pi from soil. The Pht1 gene family plays an important role in the process of P transport in plants, and express in flowers, cotyledons, pollen, leaf vascular tissues and shoots [48,49]. Pht2 and Pht3 family encode Pi transporter associated with organelle. Pht2;1 is the only member of the Pht2 family and is a low affinity phosphate transporter, which is present in the chloroplast membrane [50,51]. Pht3 family has 3 members, and the protein is located in the mitochondria. The Pht4 gene family consists of 6 members, and its structure is similar to that of P transporter SLC17/type1. The Pht4 gene expressed in both roots and leaves, 5 of which were present in plastids, and the other in Golgi, suggesting that the Pht4 gene family is associated with the transport of Pi in the cytosol, plastids and Golgi [52]. Although the function of a low affinity Pi transporter system is present in the root of the plant, the genes encoding these transporters remain to be identified. The transcriptome changes of cultivated rice to P starvation have been reported in previous study [53]. In this study, we compared the DEGs between DXWR and cultivated rice. We found two genes encoding metallothioneins (LOC_Os12g38270.2 and LOC_Os12g38290.1) were up-redulated in roots, which was same with previous research [53]. Metallothioneins affect metal tolerance and homeostasis and scavenge reactive oxygen species [54], which could be a mechanism to overcome the increase in certain ion concentration, such as iron, upon Pi starvation. Under Pi starvation, plants overaccumulated some ions, including iron [45,55,56]. We also found one gene encoded iron transporter (LOC_Os09g23300.1) was up-regulated both in roots and leaves and two genes encoded vacuolar iron transporters (LOC_Os04g45520.1 and LOC_Os09g23300.1) were up-regulated in roots and leaves, respectively (LOC_Os04g45520.1 was up-regulated in roots and LOC_Os09g23300.1 was up-regulated in leaves). In addition, LOC_Os08g06010.1, a putative glycerol 3-phosphate permease, and LOC_Os03g40670.1, a putative glycerophosphoryl diester phosphodiesterase suggested to be involved in Pi remobilization [57], were also up-regulated both in roots and leaves after P-deficiency stress. One gene encoded MYB family TF (LOC_Os02g22020.1), one encoded transporter (LOC_Os08g31670.1) and one encoded glycerophosphoryl diester phosphodiesterase family protein (LOC_ Os02g31030.1) which were up-regulated both in roots and leaves of DXWR but only up-regulated in leaves of cultivated rice.
To confirm the validity of the DEG data, quantitative RT-PCR was performed to investigate the expression patterns of 15 randomly selected genes under the same conditions. Expression trends were consistent for all transcripts in RNA-Seq analysis and quantitative RT-PCR analysis, with a correlation coefficient of (R 2 = 0.9195) (Fig. 2). Thus, the DEGs detected in this study can be considered to be a high accuracy.

Functional classification by gene ontology
The GO gene ontology mapping software (WEGO) was used to classify the function and draw gene ontology tree, and the down and up regulated transcripts in roots and leaves were classified into the functional groups. In the present study, a total of 12,180 root transcripts in RLP vs. RCK and 6803 leaf transcripts in LLP vs. LCK were assigned GO terms. Among the 12,180 root transcripts, 1829 were annotated for their molecular function, 5865 transcripts were annotated for their cellular component,

OsPHR2 LOC_ Os07g25710
None None [35] OsUPS LOC_ Os03g13740 None None [36] OsPI1 LOC_ Os05g34940 None None [37] OsSPX1 LOC_ Os06g40120 Up ( OsRCI2-9 LOC_ Os06g44220 Up (4.96) None [43] OsbHLH172 LOC_ Os06g12210 Up (2.87) None [44] and 4486 were annotated for their biological process (Fig. 3). Among the 6803 leaf transcripts, 1241 were annotated for their molecular function, 2816 transcripts were annotated for their cellular component, and 2746 were annotated for their biological process (Fig. 4). In the classification of biological processes, cellular processes and metabolic processes were the most functional groups, which indicated that the DXWR had a wide range of metabolic activities under the P-deficiency stress. In the classification of cell components, the transcripts of cells, cell components and organelles were the most abundant. In the classification of molecular function, the transcription of immobilized and active were the most highly represented groups. We also identified biological process, cellular component and molecular function GO terms that were overrepresented (p < 0.05) among the DEGs of LLP vs. LCK and RLP vs. RCK, respectively (Additional file 10: Table  S9, Additional file 11: Table S10). It was found that the expression of transcripts in the extracellular matrix increased in roots and leaves for the cellular component category (GO: 0005576, annotated as extracellular region), suggesting that the same cell components in the root and leaf are involved in the P-deficiency response.

Kyoto encyclopedia of genes and genomes (KEGG) pathway mapping
KEGG pathway analysis showed that 700 of the 1184 leaf DEGs and 622 of the 986 root DEGs could be classified into 20 functional categories and 114 and 117 subcategories, respectively. Furthermore, the over-represented KEGG Orthology (KO) terms (Q < 0.05) could be classified into 7 and 4 categories, respectively (Fig. 5). As shown in Fig. 5, the most common KO terms represented by both the leaf and root DEGs was the translation KEGG pathways. The over-represented KO terms for the leaf and root DEGs were further classified into 16 and 17 subcategories, respectively (Additional file 12: Table S11,  Table S12). Among these subcategories, 2 subcategories were over-represented among the leaf DEGs, which are RNA transporters and mRNA monitoring pathways, suggesting that these pathway may regulate the expression of P-deficiency inducible genes. Leaf and root did not have the same subcategory. Moreover, 12 KO terms were exclusively enriched among the leaf DEGs and 7 KO terms were exclusively enriched among the root DEGs. This finding suggests that there could be considerable differences in the biochemical and physiological processes involved in the P-deficiency responses of leaves and roots, and these annotations provide a valuable resource for investigating the specific processes, functions, and pathways involved in such differences.

DEGs mapped to previously identified P-deficiency responses related QTL intervals
Map based cloning is a method for identifying potential QTLs sites. Many QTLs have been identified in DXWR were associated with yield and multiple resistance. Based on the Gramene QTL database, a total of 57 QTLs related to P-deficiency stress have been found. We located 278 genes differentially regulated by P-deficiency stress on 10 of these identified QTL intervals in rice (Additional file 14: Table S13). Among them, the QTL Accession AQBD003, AQCI008, AQCI013, and AQCI011 had the greatest number of co-localized DEGs with 86, 33, 31 and 29 genes, respectively (Additional file 15: Table S14).
In the present study, 13 DEGs were co-localized within the qFWS-4 interval (Additional file 16: auxin-responsive genes [62]. And another one gene encoded homeobox domain containing protein (LOC_ Os11g06020.1), whose products are transcription factors sharing a characteristic protein fold structure that bound DNA [63][64][65], regulating gene transcription. Furthermore, another gene, LOC_Os11g44310.1, which encoded calmodulin binding protein, regulated the expression of calmodulin to regulate the intracellular calcium concentration to control cells in many important biochemical reactions, which might play an important role in stress responses [66]. 11 DEGs were co-localized within the qFWR-4 interval (Additional file 17: Table S16), including one gene encoding phosphoglycerate mutase (LOC_Os11g05260.1), which was the key enzyme in sugar metabolism process, regulating the adaptation of plants to environment. Another one encoded nucleoside-triphosphatase (LOC_Os11g03290.1) [67], which maintained the connection between the nucleus cytoplasm and cytoplasm, including transcription and regulating biosynthesis to adapt to the environment. The interval also included LOC_Os11g05400.1, which encoded Ser/Thr protein phosphatase family protein, modified other proteins by chemically adding phosphate groups to them and, thus, regulated cellular pathways, signal transduction, and responses to biotic and abiotic stresses [68]. Meanwhile, LOC_Os11g04300.1 and LOC_Os11g03940.1 encoded retrotransposon proteins, whose replicative mode of transposition by means of an RNA intermediate rapidly increased the copy numbers of elements and thereby could increase genome size, which could induce mutations by inserting near or within genes to response to biotic and abiotic stresses. LOC_Os11g04290.1 and LOC_ Os11g05380.1, which encoded cytochrome P450, modified other proteins by chemically adding OH group to them or as an enzyme catalyzed reaction and [69], thus, regulated cellular pathways, signal transduction, and responses to biotic and abiotic stresses.

Conclusion
DXWR has a lot of useful agronomic traits we need, such as resistance of cold, drought, salt and P-deficiency. Therefore, it is considered to be an important resource for rice breeding. In this study, we analyzed the transcriptome of roots and leaves of DXWR seedlings under P-deficiency stress. A large number of DEGs and some critical paths were identified, such as RNA transport and mRNA monitoring path. By combining the DEGs identified in the present study with previously identified P-deficiency resistance QTLs from rice, important candidate genes were identified, including a variety of transcription factors and some functional protein genes. These findings will be useful in the future studies of molecular adaptations to P-deficiency stress and will facilitate the genetic manipulation of rice to improve its P-deficiency resistance.