Integrated mRNA and miRNA expression profiling analyses of Pinus massoniana roots and shoots in response to phosphate deficiency

doi:10.21203/rs.3.rs-16720/v1

Download PDF

Research article

Integrated mRNA and miRNA expression profiling analyses of Pinus massoniana roots and shoots in response to phosphate deficiency

https://doi.org/10.21203/rs.3.rs-16720/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background Masson pine ( Pinus massoniana ) is primarily present in subtropical and tropical areas of China, which are severely deficient in inorganic phosphate (Pi). Although some studies identified transcriptomic and proteomic responses to Pi deficiency in Masson pine seedlings, different tissues, especially the roots, that exhibit primary responses to low-Pi stress, have not been well studied. To shed further light on the complex responses of Masson pine to Pi deficiency, a spatiotemporal experiment was performed to identify differentially expressed mRNAs and miRNAs under Pi-deficient conditions.

Results Spatiotemporal analyses of 72 RNA sequencing libraries provided a comprehensive overview of the dynamic responses of Masson pine to low-Pi stress. Differentially expressed gene analysis revealed several high-affinity phosphate transporters and nitrate transporters, reflecting the crosstalk between nitrate and Pi homeostasis in plants. The MYB family was the most abundant transcription factor family identified. miRNA differential expression analysis identified several families that were associated with Pi deficiency, such as miR399. In addition, some other families, including novel miRNA families in Masson pine, were dramatically changed in response to Pi starvation. GO and KEGG analyses of these mRNAs and targets of miRNAs indicated that metabolic processes were most enriched under Pi deficiency.

Conclusions This study provided abundant spatiotemporal transcriptomic information to functionally dissect the response of Masson pine seedlings to Pi deficiency, which will aid in further elucidation of the biological regulatory mechanisms of pines in response to low-Pi stress.

Plant Physiology and Morphology

Plant Molecular Biology and Genetics

Pinus massoniana

Phosphate deficiency

Transcript analysis

miRNA

Targets

Phosphorus (P) is an essential macronutrient and a major limiting factor of plant growth, with crucial role in metabolism, development and reproduction (Abel et al., 2002). The bioavailable form of phosphorus for direct absorption by plants is orthophosphate, which is often limited due to high fixation with other soil constituents and slow diffusion (Shen et al., 2011). Thus, a deficiency in available phosphorus in soils may constrain the development of the plants that grow in that soil (Vance et al, 2003). In the past few decades, the application of inorganic phosphate (Pi) fertilizers has been the main alternative to maintain or increase growth of crop and forest. However, plants can utilize only less than one third of Pi fertilizers because they exist in an inaccessible form (Jez et al., 2016). Furthermore, as a non-renewable resource, Pi fertilizer applied excessively can result in severe environmental pollution (Andersson et al., 2013). Therefore, to reduce the high cost of Pi fertilizers and render sustainable agriculture development, it is necessary to understand the strategies involved in regulating phosphorous homeostasis in plants and then breed new crop or forest varieties with high phosphorus acquisition and use efficiency (Veneklaas et al., 2012; Liu et al., 2016).

Plants have evolved a series of response strategies that enable them to adapt and manage a deficiency in Pi, including the modification of plant growth, alteration in root architecture, and secretion of organic acids, nucleases and phosphatases to increase Pi acquisition (Secco et al., 2013). In addition, some metabolic pathways are remodelled, and macromolecules containing phosphorus are degraded under Pi-deficient conditions (Ren et al., 2018). These adaptive strategies have been mediated synergistically by a series of genes, primarily through the regulation of signal transduction, phosphorus acquisition, internal remobilization, and metabolic changes (Fang et al, 2009). Extensive research on the response of plants to Pi starvation has considerably improved the knowledge of the complex mechanisms regulating phosphorus signalling and homeostasis (Hammond et al., 2003; Jain et al., 2012; Gho et al., 2018).

As major proteins for the uptake and translocation of phosphorus, phosphorus transporters were first identified in Arabidopsis thaliana (Muchhal et al., 1996), and an increasing number of PHT1 genes involved in transporting phosphorus have been found in this species (Shin et al., 2004; Lapis-Gaza et al., 2014; Abel, 2017). Previous studies identified the key role of the MYB transcription factor PHOSPHATE STARVATION RESPONSE1 (PHR1) gene in the transcriptional responses of plants to Pi starvation and found that the proteins containing SYG1/PHO81/XPR1 (SPX) domains were vital for plants to sense Pi and maintain homeostasis (Zhang et al., 2014; Wild et al., 2016; Zhao et al., 2018). Moreover, miR399 and miR827 can also positively regulate phosphorus-related genes (Chien et al., 2017). As one of the most abundant and highly upregulated microRNAs (miRNAs) by Pi deficiency, miR399 represses its target gene PHOSPHATE2 (PHO2), which encodes a ubiquitin-conjugating E2 enzyme that in turn negatively regulates phosphate transporters (PHTs) (Park et al., 2014). A number of other protein-coding genes involved in Pi homeostasis were induced by Pi deficiency, including purple acid phosphatases (PAPs), UDP-sulfoquinovose synthases, nonspecific phospholipases (Duff et al., 1994; Lan et al., 2012), and other miRNAs, such as miR156 and miR2111 (Hsieh et al., 2009).

Each species and even different tissues and genotypes exhibit individual response patterns to Pi deficiency (Ren et al., 2018). The Masson pine (Pinus massoniana) is a native gymnosperm in southern portions of China that suffer from a critical lack of Pi (Zhang et al, 2010). Previously, we identified a number of differentially expressed genes (DEGs) and proteins involved in multiple metabolic pathways under P-deficient conditions in Masson pine seedlings at the transcriptome and proteome levels (Fan et al., 2014; 2016). In addition, some genes with unknown functions were up- and downregulated. However, previous studies focused on whole seedlings, different tissues, especially the roots, which exhibit primary responses to low-Pi stress, have not been well studied. Furthermore, in addition to mRNA regulation, miRNAs function as key regulators of nutrient stress signalling (Chiou, 2007). However, there have been no reports of Masson pine under Pi-deficient stress to date, revealing the need for the identification of low-Pi responsive miRNAs in Masson pine. Therefore, functional analyses of candidate genes and miRNAs in different tissues are necessary.

In this study, to shed further light on the complex responses of Masson pine to Pi deficiency, next-generation RNA-seq analysis coupled with a spatiotemporal experiment was performed to identify differentially expressed mRNAs and miRNAs (DE mRNAs and DE miRNAs, respectively) under Pi-deficient conditions. Additionally, a comprehensive and integrated dataset was generated to demonstrate the mechanism by which the transcripts responsive to phosphate starvation bolster the efficiency of the acquisition, recycling, and remobilization of phosphorus in roots.

Analyses of the Pi concentration and mRNA and small RNA sequencing of Masson pine in response to Pi deficiency

To better understand the complex mechanisms of Pi homeostasis in Masson pine seedlings, time-course experiments were conducted to evaluate the Pi concentration and mRNA and small RNA expression. In terms of the morphological and Pi concentration analyses, there was no significant growth difference in the aerial parts of the control and low-Pi-treated seedlings. The primary root length and lateral root number of the low-Pi-treated seedlings significantly increased compared with those of the control (Figure 1A). The Pi concentration analysis showed a gradual decrease as the duration of the treatment increased. In addition, Pi deficiency for 48 d resulted in a rapid decrease in the internal Pi concentration (Figure 1B). To identify transcripts, including mRNAs and miRNAs, involved in Pi-deficiency responses in Masson pine roots and shoots and to assess the dynamic responses of these RNAs to Pi levels, we selected three long-term time points and used three biological replicates per Pi condition for sequencing, representing a total of 72 libraries (including 36 RNA-seq libraries and 36 sRNA-seq libraries). Sequencing libraries were generated and subjected to Illumina deep sequencing.

Overall, for RNA-seq, more than 3.0 billion clean reads were generated after filtering low-quality raw reads and removing rRNAs. On average, each sample had 84.0 million high-quality paired-end reads with Q20 and Q30 values of more than 98.5% and 96.0%, respectively (Table S2). Through the Trinity de novo assembly method, a total of 98,472 unigenes were obtained, with an N50 length of 1,163 bp and an average length of 757 bp. The lengths of the unigenes ranged from 201 to 12,445 bp, and the distribution of their lengths is shown in Figure 1C. The resultant unigenes were queried in the Nr, Swiss-Prot, KEGG, and KOG databases according to the sequence similarity analysis. Among the 98,472 unigenes, approximately 45,386 unigenes were identified in plants, and the Nr database provided the largest number of annotations, with a value of 45,078 (45.78%). In total, 12,925 (13.13%) unigenes had shared annotations between the four databases (Figure S1).

sRNA deep sequencing from Masson pine roots and shoots under Pi-deficient or Pi-sufficient conditions yielded a total of 604,281,401 clean reads. After removing contaminants and reads smaller than 18 nt, over 451.96 million clean sRNAs were generated. Table S3 shows an overview of the raw sRNA-seq reads, high-quality and clean reads with quality filtering. The distributions of miRNA lengths were similar among the 36 sRNA-seq libraries, and 21 nt sRNAs were the most abundant (Figure 1D), which is consistent with previous study on Populus (Ren et al., 2015; Bao et al., 2019).

Expression profiles of mRNAs in response to Pi starvation

Transcript abundance was estimated by the RPKM approach. More than 26 thousand expressed genes were detected in this study, and several thousands of these genes showed tissue-specific expression under different growth conditions (Figure 2A). Differential expression analysis via DEG-Seq showed that many genes were significantly differentially expressed in either roots or shoots from 24 d to 48 d of Pi deficiency (P < 0.05, FDR < 0.05). In addition, the number of DEGs, especially the upregulated genes, increased over time under Pi starvation in either the roots or shoots (Figure 2B).

The longest duration (48 d) of Pi starvation led to a significant decrease in the internal Pi content by more than twofold (Figure 1B), but the number of differentially regulated genes increased substantially in the shoots and roots. At this time point, many DEGs related to phosphorus homeostasis were induced in the roots and shoots, and the fold-change values are displayed in Table 1. These DEGs were involved in transport, transcription, phosphorylation/dephosphorylation, lipid metabolism, metabolism, and miscellaneous processes (Lan et al., 2015; Ren et al., 2018). In terms of transporters, DEGs defined as PHTs were identified in both the roots and shoots, and more PHTs were upregulated in the shoots than in the roots. Moreover, one nitrate transporter gene (Unigene0075150) was upregulated in the shoots. In the transcription factor group, all five SPX domain-containing proteins (SPXs) were differentially expressed in the roots and shoots. One ethylene responsive factor (ERF, Unigene0018015), which is crucial for root development under Pi-deficient conditions, was induced in Masson pine roots under low-Pi stress. PAPs and other acid phosphatases of the phosphorylation/dephosphorylation group, which can increase the excretion of organic acids and Pi-releasing enzymes, were upregulated in both the roots and shoots. In lipid metabolism, monogalactosyldiacylglycerol synthase (Unigene0022807), which can result in the degradation of phospholipids, was induced in the shoots under Pi-starvation conditions. In addition, other genes related to metabolism and disease resistance were sharply upregulated by Pi deficiency.

The DEGs of the Masson pine roots and shoots were further clustered using Short Time-series Expression Miner (STEM) with default parameters. Eight significant gene profiles (P-value ≤ 0.01) were identified in response to Pi deficiency in both the roots and shoots (Figure 2C, D). The results indicated that STEM profile 0 (48 genes in the roots and 49 genes in the shoots) and profile 7 (131 genes in the roots and 103 genes in the shoots) were down- and upregulated, respectively, as the duration of stress increased. STEM profile 2 (264 genes in the roots and 312 genes in the shoots) and profile 5 (66 genes in the roots and 58 genes in the shoots) exhibited a reversal in the trend of expression between the stress time points. Profile 3 (153 genes in the roots and 171 genes in the shoots) and profile 4 (194 genes in the roots and 278 genes in the shoots) gene clusters were down- and upregulated, respectively, at 48 d, and no changes at earlier stages were detected. Profile 1 (51 genes in the roots and 41 genes in the shoots) and profile 6 (277 genes in the roots and 205 genes in the shoots) showed no change between 36 d and 48 d.

GO and KEGG enrichment analyses of the DEGs

To illustrate the function of the persistent DEGs in the three stress time points, GO term enrichment analysis was conducted, and all differentially regulated transcripts in the Pi-deficient roots and shoots were assessed based on the categories biological process, cellular component, and molecular function (Figure 3). For the biological process ontology, the dominant terms in this study included “metabolic process”, “single-organism process”, “cellular process”, and “response to stimulus” in both the roots and shoots. The dominant terms for cellular component were “cell”, “cell part”, “organelle”, and “macromolecular complex” in the roots, while “cell”, “cell part”, “membrane”, and “organelle” were dominant in the shoots. The processes represented by the terms “catalytic activity”, “binding”, and “transporter activity” accounted for the majority of the molecular function category, followed by “transporter activity” in both the roots and shoots.

To obtain a precise understanding of the DEG functions in the Pi-limited roots and shoots, a KEGG enrichment analysis was conducted. Among the top 20 KEGG enrichments, four pathways (biosynthesis of secondary metabolites, ribosome, phenylpropanoid biosynthesis, and tyrosine metabolism) were the most abundant in the Pi-deficient roots, while metabolic pathways, biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, and plant hormone signal transduction were the four most abundant pathways in the shoots (Figure 4). However, as shown in Figure 4, the pathways tyrosine metabolism, flavonoid biosynthesis, and alpha-linolenic acid metabolism in the roots and the pathways phenylpropanoid biosynthesis, plant hormone signal transduction, and flavonoid biosynthesis in the shoots were over-represented among the DEGs, suggesting that these pathways could be involved in the response to low Pi. Other pathways that were closely linked to phosphorus were further detected, such as cutin, suberine and wax biosynthesis, betalain biosynthesis, isoquinoline alkaloid biosynthesis, pentose and glucuronate interconversions, and amino sugar and nucleotide sugar metabolism.

Identification of known and novel miRNAs

Through a BLAST search with the GenBank and Rfam databases, the obtained sRNAs, including rRNA, scRNA, sonRNA, snRNA, and tRNA, were removed. The remaining sequences were aligned with miRbase v. 21.0, and unique sRNAs matching known miRNAs were identified as shown in Table S4. Finally, the unannotated unique sRNA sequences of each sample were retained and utilized to predict novel miRNA candidates (Table S5). In total, we identified 88 known miRNAs belonging to 60 miRNA families, and 586 predicted novel miRNAs belonging to 511 families in the small RNA libraries (Table S6, Table S7). Details regarding the mature ID, hairpin information, sequence, and length of each known and novel miRNA are summarized in Table S6 and Table S7. Among the known miRNA families, approximately 19 families contained more than one member.

Expression profiles of miRNAs in response to Pi starvation

The DE miRNAs were identified via the formula log₂ (treatment/control) with an absolute value ≥ 1 and a P-value ≤ 0.05. The results showed that some known and novel miRNAs were significantly differentially expressed in either the roots or shoots from 24 d to 48 d of Pi deficiency (Figure 5A, Table S8). In addition, the expression patterns of the DE miRNAs were similar to those of the DEGs, especially in the roots, and the number of DE miRNAs increased with the duration of Pi starvation. The longest (48 d) Pi-starvation treatment led to a significant increase in the number of DE miRNAs in the roots. At this time point, several known miRNA families related to the response to phosphorus conditions were differentially expressed in the roots, such as miRNA399, miRNA169, miRNA398, and miRNA408, and the fold-change values are displayed in Table S8.

To directly compare the expression patterns of these miRNAs at the three time points under low-Pi stress, eight specific profiles were identified as significant via the STEM cluster analysis (Figure 5B). The results indicated that 42 miRNAs in the roots and 20 miRNAs in the shoots were clustered in different profiles (Figure 5C). Profile 0 (downregulated, 3 miRNAs in the roots and 3 miRNAs in the shoots) and profile 7 (upregulated, 8 miRNAs in the roots and 0 miRNAs in the shoots) followed the stress-time trend. Profile 2 (7 miRNAs in the roots and 6 miRNAs in the shoots) and profile 5 (4 miRNAs in the roots and 3 miRNAs in the shoots) exhibited a reversal in expression profiling trends.

Integrated analysis of DE miRNAs and mRNAs in response to Pi deficiency

To define as many low-Pi-induced miRNA-mRNA interactions as possible, expression correlation between DE miRNAs and DE mRNAs at 48 d, at which the greatest number of differentially expressed transcripts were identified, was evaluated using the PCC. Target prediction results showed that 921 targets were identified for 26 known miRNAs, and 442 targets were found for 46 novel miRNAs (Table S9). KEGG enrichment analysis for the target mRNAs showed the top 20 pathways in response to Pi deficiency, with “spliceosome” the most frequently represented subclass (Figure 6A). In general, miRNAs negatively regulate the expression of their target mRNAs by either mRNA cleavage or translational repression. Therefore, the simple expectation is that when miRNAs are induced by Pi deficiency, their target mRNAs are reduced and vice versa. In total, twelve negative miRNA-mRNA interactions and five unique novel miRNAs and their 11 target mRNAs were identified in response to 48 d of Pi deficiency (Figure 6B, Table S10). Pathway enrichment analysis for the negative miRNA-mRNA pairs revealed that the “carbohydrate metabolism” pathway significantly changed (P-value < 0.05) (Figure 6C).

qRT-PCR validation of the DE mRNAs and DE miRNAs

The expression levels of eight DE miRNAs and the 11 DEGs (Table S1) from the miRNA-mRNA interaction network (Figure 6B), together with other randomly selected miRNAs and mRNAs, were further validated using qRT-PCR. The results revealed that most of these miRNAs and mRNAs shared similar expression levels with the miRNA-seq and mRNA-seq data (Figure 7), suggesting that the RNA-seq data were reproducible and reliable.

Phosphorus is an increasingly limiting factor for plant growth and yield worldwide as a result of a deficiency in available phosphate in the soil. The use of phosphorus fertilizers has alleviated the problem, but it is rapidly depleted as a non-renewable natural resource, and lead to serious environmental pollution (Tong et al., 2017). Therefore, a substantial number of studies based on next-generation sequencing have documented the molecular responses of different plant species under Pi-deficient conditions (O′Rourke et al., 2013; Gho et al., 2018; Huen et al., 2018). Previously, we reported temporal transcriptome responses to Pi starvation in Masson pine seedlings and found that most transcripts were stably and gradually up- or downregulated within 48 d of Pi starvation (Fan et al., 2014). However, relatively few functional analyses of the transcripts of different tissues, especially candidate root genes and miRNAs exhibiting primary responses to Pi-deficient conditions, have been performed. In this study, three time points (24, 36, and 48 d) were selected for the investigation of all transcript (mRNA and miRNA) responses to Pi deficiency in Masson pine shoots and roots using next-generation sequencing methods. The results showed that the morphological changes were consistent with previous observations, and a large number of DEGs and DE miRNAs were identified. Some reports have demonstrated that the miRNA-mRNA regulatory network responds to abiotic stress (Xie et al., 2017). In this work, we attempted to construct a low-Pi-induced miRNA-mRNA regulatory network according to the DEG and DE miRNA datasets. qRT-PCR validation results showed that the expression changes generally agreed with the results of RNA-seq and confirmed that the identified DEGs and DE miRNAs were reliable.

Analysis of the morphological and physiological responses

Over the past two decades, various studies have reported that low Pi concentrations inhibit primary root growth and enhance lateral roots and root hairs. Recently, Zhang et al. (2019) demonstrated that the direct illumination of the root surface with blue light is critical and sufficient for inhibiting the growth of primary roots under Pi-deficient conditions. However, root growth in soil could be quite different from transparent Petri dishes because of the difference in light penetration. Some investigators have found that the plant primary root length increased in Pi-deficient soil (Anuradha and Narayanan, 1991; Foti et al., 2014; Strock et al., 2018). In this study, the primary root length and lateral root number of Masson pine seedlings increased under Pi-deficient conditions, which is similar to the results of our previous work (Fan et al., 2016). The Pi concentration can reflect the absorption ability of plants under different Pi conditions (Ren et al., 2018). In the present study, the Pi concentration of seedlings decreased gradually as the duration of the treatment increased, whereas the aerial parts of Masson pine seedlings exhibited no observable growth differences to the naked eye between the low-Pi treatment and control seedlings (Figure 1A). In addition, the degree of the decrease in the Pi concentration revealed no significant difference under Pi-deficient conditions compared with the control conditions (Figure 1B). These results indicate the acclimation of the Masson pine to Pi-deficient conditions.

Sensing Pi-deficiency signals in Masson pine

A previous study reported that plants are able to sense external Pi concentration changes by a root cell membrane-localized sensor and sense internal nutrient status by an intracellular sensor. In addition, plasma membrane-localized nutrient transporters and related proteins are considered nutrient sensors (Zhang et al., 2014). Based on this concept, Pi transporters (PHT1 family) can work as external Pi sensors (Zhang et al., 2014). In our study, nine upregulated PHT1s were identified under Pi-deficient conditions for 48 d, three of which were differentially expressed in the roots (Table 1). Based on these findings, one of these PHT1s may sense the external Pi concentration. Future studies are required to identify which of these PHT1s function as Pi sensors and the mechanism underlying the Pi sensor function.

SPXs (including SPX domain-containing protein 1 (SPX1) to SPX domain-containing protein 1 (SPX4)) play an important role in Pi sensing, signalling, and transport in plants. As an internal phosphate sensor, SPX4 was reported to be the only SPX weakly repressed under Pi-starvation conditions (Osorio et al., 2019). Ruan et al. (2019) identified and demonstrated two RING-finger ubiquitin E3 ligases facilitating the degradation of SPX4 to modulate PHR2 activity and regulate Pi homeostasis in rice. In the present study, five SPXs were differentially expressed under Pi starvation for 48 d, and only one N-terminal SPX was suppressed. Ubiquitin E3 ligase, which can facilitate the degradation of SPX4, was upregulated in this study. In contrast, the gene described as SPX4 was mildly induced in both the roots and shoots (Table 1). This result was different from the results of other investigations (Lv et al., 2014; Osorio et al., 2019; Ruan et al., 2019).

Transporters related to phosphorus metabolism

Specific membrane transport systems are needed for plants to acquire Pi from the soil (Schroeder et al., 2013). Previous studies identified a number of phosphorus transporters with functions in the uptake, transport, and translocation of phosphorus in plants. Most transporters that are in the high-affinity Pi transporter (PHT1) family are known to be differentially expressed in wheat, soybean, rice, and Brachypodium distachyon (Miao et al., 2009; Qin et al., 2012; Gho et al., 2018; Zhao et al., 2018). In this study, differentially expressed PHT1s were identified in both the roots and shoots, suggesting the function of phosphorus transport from the roots to the shoots. In addition to PHTs, other transporter families were found to be significantly changed under Pi-deficient conditions. Nitrate transporters have been shown to play an important role in regulating Pi-deficiency responses in Arabidopsis (Cui et al., 2019). As a nitrate transporter, NRT1.1B was indicated to interact with and degrade the phosphate signalling repressor SPX4 and then activate phosphate utilization in plants (Hu et al., 2019). In the present study, one nitrate transporter (Unigene0075150) was upregulated in the shoots under Pi-deficient conditions; however, whether this enables the absorption of Pi by Masson pine merits investigation.

Metabolic adaptations in response to low-Pi stress

Metabolic changes are essential for the acclimation of plants to Pi-deficient conditions. Among these, sugar and lipid metabolic pathways are two well-studied adaptations to Pi deficiency in plants (Zhang et al., 2014). Sucrose and starch have been reported to be increased under Pi starvation in plants (Park et al., 2012; Zhang et al., 2014). Carbon fixation is critical for the synthesis of sucrose and starch. As a biosynthesis-related enzyme, phosphoenolpyruvate carboxylase (PEPC) plays a major role in fixing carbon for organic acid synthesis in response to Pi deficiency in plants (Peñaloza et al., 2005). Later, PEPC was speculated to play a role in producing sucrose (Daloso et al., 2016). In this study, PEPCs were upregulated in both the roots and shoots, suggesting a function in regulating sugar metabolism in response to Pi deficiency.

Pi-deficiency may result in dramatical altering of lipid composition of the plant membranes, with a decrease of phospholipids and an increase of phosphorus-free lipids (Nakamura, 2013). Generally, Pi-deficiency induces the replacement of membrane phospholipids by galactolipids, sulfolipids, and glycolipids (Shemi et al., 2016; Tawaraya et al., 2018). Glycerophosphodiester phosphodiesterases (GDPDs) can hydrolyse intermediate products of phospholipid catabolism and have been reported to be induced by Pi starvation in plants (Mehra and Giri, 2016; Mehra et al., 2019). As an important sulfolipid biosynthesis enzyme, UDP-sulfoquinovose synthase (SQD1) was induced under Pi-starvation conditions in plants (Nakamura, 2013). In the present study, GDPDs and SQD1 were significantly upregulated in both the roots and shoots. The induction of these genes suggested that the membrane lipid composition might change under low-Pi stress. Previous studies have reported that type-B monogalactosyldiacylglycerol synthases (MGDs) are involved in lipid remodelling in response to low Pi (Kobayashi et al., 2009). The type-A MGD1 gene was upregulated in Masson pine shoots in this study, revealing that this class of MGDs may play a role in the plant response to Pi-deficient conditions.

Transcriptional regulation in response to Pi starvation

Our previous study demonstrated that the MYB transcription factor family represented the largest class of transcription factors induced by a Pi deficiency in Masson pine seedlings (Fan et al., 2014); the same result was obtained in this study. As an MYB transcription factor, PHR1 was previously reported to be a central player in Pi-starvation transcriptional responses (Zhang et al., 2014) and has been shown to be involved in Pi deficiency by activating a series of genes, including those that contain SPX domains, and miRNA399 (Bari et al., 2006; Secco et al, 2012). SPX proteins have been reported to negatively or positively regulate phosphate-starvation-inducible genes (Duan et al., 2008; Jung et al., 2018). In this study, SPXs were differentially expressed in both the roots and shoots, indicating that a series of related genes would be induced by Pi deficiency. As a group of APETALA₂ (AP₂) domain-containing transcription factors, ERFs involve in ethylene-mediated transcription as either repressors or activators, contributing to the regulation of root system architecture (RSA) and both the biosynthesis and accumulation of anthocyanins in plants under low Pi conditions (Chen et al., 2018). This study revealed that ERFs were upregulated in response to low-Pi stress in both the roots and shoots; however, whether these confer Pi-deficiency tolerance to Masson pine merits investigation.

Pi-starvation-responsive miRNAs, and the integrated analysis of miRNAs and mRNAs in Masson pine

Previous studies identified several miRNAs as new participants in the regulation of Pi-starvation-responsive genes at the posttranscriptional level. miR399 was the first miRNA identified to be upregulated by Pi deficiency and rapidly decreased in the present of a sufficient Pi supply (Fujii et al., 2005). It has been reported that miR399 could direct the cleavage of PHO2, which negatively regulates Pi transporters to inhibit Pi absorption, to regulate phosphate homeostasis in plants (Bari et al., 2006; Ouyang et al., 2016). Following the identification of miR399, other Pi-starvation-responsive miRNA families have been identified, including miR827 and miR211, which are specifically induced by Pi deficiency but not by other nutrient deficiencies in Arabidopsis (Kuo and Chiou, 2011). As expected, miR399 also showed regulatory effects in Masson pine roots under Pi-starvation conditions, suggesting the conserved roles of miRNAs in regulating Pi homeostasis in plants. Several miRNAs also showed dramatic decreases in Pi-deficient Masson pine roots, such as miR169, miR3695, miR5076, and miR6192. Notably, Arabidopsis miR169 was downregulated by nitrogen starvation, thereby aiding plants in adapting to nitrogen-starvation conditions (Zhao et al., 2011). miR169 has been speculated to be involved in the acclimation to Pi deficiency because of crosstalk between nitrate and Pi homeostasis in plants (Secco et al., 2013).

A transcriptome-wide integrated analysis of low-Pi-induced miRNA and mRNA interactions was performed to understand which pathways were most likely regulated by miRNAs under Pi-deficiency stress. In general, miRNAs negatively regulate the expression of their target mRNAs by either mRNA cleavage or translational repression. In total, five unique novel miRNAs and their 11 target mRNAs were identified in response to 48 d of Pi deficiency (Figure 6B, Table S10). Pathway enrichment analysis indicated that the pathway “carbohydrate metabolism” was significantly changed (P-value < 0.05), which was proven to be essential for plant acclimation to Pi-deficient conditions (Zhang et al., 2014). The differentially expressed targets, with respect to their cognate miRNAs, could be clustered into NBS-LRR, ubiquitination, pectate lyase, and so on. NBS-LRR genes were documented to enhance tolerance to abiotic stress (Li et al., 2017) and dramatically change under nitrogen deficiency in poplar (Song et al., 2019). In this study, the predicted novel miR0145-NBS-LRR interaction was identified, and its functions in response to Pi deficiency require further analysis.

In this study, transcripts analysis coupled with a spatiotemporal experiment was performed, and many DE mRNAs and DE miRNAs were identified to correlate with low-Pi stress. GO and KEGG analyses of these mRNAs and targets of miRNAs indicated that metabolic processes were most enriched under Pi deficiency. This global analysis of the DE mRNAs and miRNAs in Masson pine under Pi-deficient conditions revealed candidate genes and miRNAs that can be further characterized to elucidate the complex responses of Masson pine to adapt to Pi deficiency.

Plant materials and growth conditions

Masson pine seeds were collected from 72-family grown in a 1.5-generation seed orchard located in the Duyun Forest Farm in Guizhou Province, China; the seedlings from these seeds were previously shown to survive well under Pi-deficient conditions. Filed studies did not involve protected species, and no specific permissions were required for this location and activities. The seeds were disinfected by treatment with 5% sodium hypochlorite for 5 minutes and 70% ethanol for 10 minutes. Next, the seeds were soaked overnight in sterile water at 30°C. The germinated seeds were planted in plastic pots and grown at a day/night temperature of 28°C/18°C (14-h day/10-h night) in a plant growth chamber with a light intensity of 2000 lux and 70% relative humidity. At 15 d after emergence, the seedlings were treated with a modified Hoagland nutrient solution (Fan et al., 2014) that was deficient in Pi (0.01 mM KH₂PO₄) and sufficient in Pi (0.5 mM KH₂PO₄) as the treatment and control, respectively. The nutrient solution was added at three-day intervals. The experiment included six seedlings per treatment and was replicated three times. The roots and shoots were harvested separately at the following time points: 24 d, 36 d, and 48 d after transfer to low-Pi treatment. The harvested samples were subsequently snap frozen with liquid nitrogen and stored at –80°C. Furthermore, the sample collection was performed at a similar time of day (approximately 10:00 am) to minimize possible circadian effects.

Quantification of the Pi concentration

The cellular contents of Masson pine seedlings were released into water by incubation in strong sulfuric acid and digested for 1 h at 360℃, and the concentration of Pi was quantified by the molybdate assay according to the procedure of Ames (1996).

Total RNA isolation, library construction, and sequencing

Total RNA from the roots and shoots was extracted with the Invitrogen Plant RNA Isolation Kit based on the manufacturer’s instructions. The quality and quantity of total RNA were assessed via agarose gel electrophoresis and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). For RNA-seq library construction, mRNA was specifically enriched from the extracted RNA with oligo (dT) magnetic beads and fragmented and subsequently subjected to cDNA generation. Then, the cDNA fragments were purified with the QiaQuick PCR extraction Kit, end repaired, supplemented with poly(A), and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using the HiSeq^TM 4000 platform (Illumina, San Diego, CA, USA). The platform software was used to conduct the primary analysis of the data and base calling.

For small RNA-seq library generation, the RNA molecules in the size range of 18-30 nt were enriched by polyacrylamide gel electrophoresis (PAGE). Then, specific 3’ and 5’ Illumina RNA adapters were sequentially ligated to the small RNAs. Adapter-ligated molecules were then reverse transcribed and PCR amplified, and the PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeq^TM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Transcriptome assembly, annotation, and expression assessment

To generate a set of transcripts that were non-redundant, adapter sequences, low-quality bases and empty reads were removed, and high-quality clean reads were selected using Q20 and Q30. Trinity software was used for clean read assembly after filtering to construct the unigenes. For homology annotation, non-redundant sequences were used to search against the NCBI non-redundant protein (Nr), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), and euKaryotic Ortholog Groups (KOG) databases using the BLASTx program with a threshold of E-value < 1E-5. Protein functional annotations were then obtained according to the best alignment results.

The levels of expression for each unique transcript were calculated by quantifying the reads based on the reads per kilobase per million reads (RPKM) approach (Mortazavi et al., 2008). DEG-Seq was used to judge the DEGs with a false discovery rate (FDR) < 0.05 and the absolute value of log₂ RPKM ratio ≥1 as the threshold (Anders and Huber, 2010).

Small RNA sequence data analysis and expression profiles of miRNA

To obtain clean reads of 18-30 nt, raw reads were filtered to exclude sequences containing more than one low-quality (Q-value ≤20) base or unknown nucleotides, adapters, poly(A) tails, and reads shorter than 18 nt (not including adapters). All of the clean tags were aligned with small RNAs in the GenBank (Release 209.0) and Rfam (11.0) databases to identify and remove rRNA, scRNA, sonRNA, snRNA and tRNA. Then, the sRNA sequences were aligned with miRbase (release 21.0) to identify known miRNAs. The unannotated sequences were utilized to predict novel miRNA candidates with MIREAP-v0.2 software.

The miRNA expression level was calculated and normalized to transcripts per million (TPM): TPM=actual miRNA counts/total counts of clean tags × 10⁶. The DE miRNAs between the deficient-Pi and sufficient-Pi (control) treatments were identified via the formula log₂ (treatment/control), with an absolute value ≥ 1 and a P-value ≤ 0.05.

Gene Ontology (GO) enrichment and KEGG pathway analyses

The unigenes or targets of miRNAs with significantly altered abundance under Pi-deficient conditions were subjected to GO and KEGG analyses. GO enrichment analysis provides all GO terms that are significantly enriched in DEGs and filters the DEGs that correspond to biological functions. The targets were mapped to the terms in the GO database (http://www.geneontology.org/), and the gene numbers were calculated for every term. Then, a hypergeometric test was performed to categorize enriched GO terms. As a major public pathway-related database, the KEGG algorithm was utilized to identify significantly enriched metabolic pathways or signal transduction pathways. The adjusted standard of the significantly substantiated GO terms and enriched pathways was a P-value less than 0.05.

Integrated analysis of miRNA-seq and mRNA-seq data

To predict the genes targeted by DE miRNAs, Patmatch (V 1.2) software was used to identify miRNA binding sites with default parameters. To define the possible miRNA-mRNA interactions, the expression correlation between miRNA targets was evaluated using the Pearson correlation coefficient (PCC). We normalized all sample-matched miRNA and mRNA sequencing data and integrated the expression profiles of DE miRNAs and DE mRNAs with the DE miRNA-target information. For the construction of the miRNA-target network, pairs with PCC < -0.7 and p < 0.05 were selected as negatively co-expressed miRNA-target pairs and then visualized using Cytoscape software (v3.6.0) (http://www.cytoscape.org/).

Gene expression validation of 48-d-old roots using quantitative real-time PCR (qRT-PCR)

Small RNAs with lengths no longer than 200 nt were extracted using the miRcute Plant miRNA Isolation Kit (Tiangen, Beijing, China) based on the manufacturer’s instructions. First-strand miRNA was synthesized using the miRcute Plus miRNA First-Strand cDNA Kit, and the relative expression of miRNAs was quantified using Reverse Primer (Tiangen, Beijing, China) with miRNA-specific primers listed in Table S1. qRT-PCR was performed on a CFX96 Touch Real-Time PCR System (Bio-Rad, USA) with SYBR Green Real-time PCR Master reagents (Applied Biosystems, CA, USA) according to the provided instructions. Each sample was tested in triplicate. The expression level was calculated by the 2^-ΔΔCt method, and 5.8S rRNA was used as an internal control.

To verify the RNA-seq results, the target and non-target unigenes were randomly selected and quantitatively analysed using qRT-PCR to confirm the RNA sequencing expression profiles. Primers were designed using Primer Premier 5.0 (Table S1). SYBR Select Master Mix Kit (Applied Biosystems, CA, USA) was used for triplicate qRT-PCRs with a CFX96 Touch Real-Time PCR System (Bio-Rad, USA) according to the manufacturer’s protocol. The expression level was calculated by the 2^-ΔΔCt method, which was normalized using 18s-rRNA and UBQ as internal controls.

P, phosphorus; Pi, inorganic phosphate; PHR1, phosphate starvation response1; SPX, SYG1/PHO81/XPR1; PHO2, phosphate2; PHTs, phosphate transporters; PAPs, purple acid phosphatases; DEGs, differentially expressed genes; ERF, ethylene responsive factor; PEPC, phosphoenolpyruvate carboxylase; GDPDs, glycerophosphodiester phosphodiesterases; SQD1, UDP-sulfoquinovose synthase; MGDs, monogalactosyldiacylglycerol synthases; AP2, APETALA2; RSA, root system architecture; RPKM, reads per kilobase per million reads; PCC, Pearson correlation coefficient

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

The datasets generated and analysed during the current study are available in the supplemental materials and from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

This study was funded by the National Natural Science Foundation of China (31660185), the Natural Science Foundation of Guizhou Province, P. R. China (20161051), Innovation Driven Development Special Fund of Guangxi, P. R. China (AA17204087-4), Post-National Key Research and Development Project, P. R. China (20185261), the First Class Discipline Construction Project in Guizhou Province, P. R. China (GNYL2017007).

Authors’ Contributions

FF conceived, designed the experiments and wrote the paper. XS performed experiments. GD contributed reagents and materials. ZZ performed experiments. JT collected data. All authors have read and approved the manuscript.

Acknowledgements

Thanks English native speaker professional with a science background at AJE for revising our manuscript.

Abel S. Phosphate scouting by root tips. Curr Opin Plant Biol. 2017; 39: 168-177.

Abel S, Ticconi CA, Delatorre CA. Phosphate sensing in higher plants. Physiol Plantarum. 2002; 115: 1-8.

Ames B. Assay of inorganic phosphate, total phosphate and phosphatases. Methods Enzymol. 1966; 8: 115-118.

Anders S, Huber W. Differential expression analysis for sequence count data, J. Genome biology. 2010; 11: R106

Andersson H, Bergström L, Djodjic F, Ulén B, Kirchmann H. Topsoil and subsoil properties influence phosphorus leaching from four agricultural soils. J Environ Qual. 2013; 42: 455-463.

Anuradha M and Narayanan A. Promotion of root elongation by phosphorus deficiency. Plant Soil. 1991; 136: 273-275.

Bao H, Chen H, Chen M, Xu H, Huo X, Xu Q, Wang Y. Transcriptome-wide identification and characterization of microRNAs responsive to phosphate starvation in Populus tomentosa. Funct Integr Genomics. 2019; 19: 953-972.

Bari R, Pant BD, Stitt M, Scheible WR. PHO2, microRNA399, and PHR1 define a phosphate-signaling pathway in plants. Plant Physiol. 2006; 141: 988-999.

Chen Y, Wu P, Zhao Q, Tang Y, Chen Y, Li M, Jiang H, Wu G. Overexpression of a phosphate starvation response AP2/ERF gene from physic nut in Arabidopsis alters root morphological traits and phosphate starvation-induced anthocyanin accumulation. Front Plant Sci. 2018; 9: 1186.

Chien PS, Chiang CB, Wang Z, Chiou TJ. MicroRNA-mediated signaling and regulation of nutrient transport and utilization. Curr Opin Plant Biol. 2017; 39: 73-79.

Chiou TJ. The role of microRNAs in sensing nutrient stress. Plant Cell Environ. 2007; 30: 323-332.

Cui YN, Li XT, Yuan JZ, Wang FZ, Wang SM, Ma Q. Nitrate transporter NPF7.3/NRT1.5 plays an essential role in regulating phosphate deficiency responses in Arabidopsis. Biochem Bioph Res Co. 2019; 508: 314-319.

Daloso DM, dos Anjos L, Fernie AR. Roles of sucrose in guard cell regulation. New Phytol. 2016; 211: 809-818.

Duan K, Yi K, Dang L, Huang H, Wu W, Wu P. Characterization of a sub-family of Arabidopsis genes with the SPX domain reveals their diverse functions in plant tolerance to phosphorus starvation. Plant J. 2008; 54: 965-975.

Duff SMG, Sarath G, Plaxton WC. The role of acid-phosphatases in plant phosphorus-metabolism. Physiol Plantarum. 1994; 4: 791-800.

Fan F, Cui B, Zhang T, Qiao G, Ding G, Wen X. The temporal transcriptomic response of Pinus massoniana seedlings to phosphorus deificiency. PLoS ONE. 2014; 9: e105068.

Fan F, Ging G, Wen X. Proteomic analyses provide new insights into the response of Pinus massoniana seedlings to phosphorus deficiency. Proteomics. 2016; 16: 504-515.

Fang ZY, Shao C, Meng YJ, Wu P, Chen M. Phosphate signaling in Arabidopsis and Oryza sativa. Plant Sci. 2009; 176: 170-180.

Foti M, Dho S, Fusconi A. Root plasticity of Nicotiana tabacum in response to phosphorus starvation. Plant Biosystems-An International Journal Dealing with all Aspects of Plant Biology. 2016; 150: 429-435.

Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK. A miRNA involved in phosphate-starvation response in Arabidopsis. Curr Biol. 2005; 15: 2038-2043.

Gho YS, An G, Park HM, Jung KH. A systemic view of phosphate starvation-responsive genes in rice roots to enhance phosphate use efficiency in rice. Plant Biot Rep. 2018; 12: 249-264.

Hammond JP, Bennett MJ, Bowen HC, Broadley MR. Changes in gene expression in Arabidopsis shoots during phosphate starvation and the potential for developing smart plants. Plant Physiol. 2003; 132: 578-596.

Hsieh LC, Lin SI, Shih AC, Chen JW, Lin WY, Tseng CY, Li WH, Chiou TJ. Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiol. 2009; 151: 2120-2132.

Hu B, Jiang Z, Wang W, Qiu Y, Zhang Z, Liu Y, Li A, Gao X, Liu L, Qian Y, Huang X, Yu F, Kang S, Wang Y, Xie J, Cao S, Zhang L, Wang Y, Xie Q, Kopriva S, Chu C. Nitrate-NRT1.1B-SPX4 cascade integrates nitrogen and phosphorus signaling networks in plants. Nat Plants. 2019; 5: 401-413.

Huen A, Bally J, Smith P. Identification and characterisation of microRNAs and their target genes in phosphate-starved Nicotiana benthamiana by small RNA deep sequencing and 5’RACE analysis. BMC Genomics. 2018; 19: 940.

Jain A, Nagarajan VK, Raghothama KG. Transcriptional regulation of phosphate acquisition by higher plants. Cell Mol Life Sci. 2012; 69: 3207-3224.

Jez JM, Lee SG, Sherp AM. The next green movement: plant biology for the environment and sustainability. Science. 2016; 353: 1241-1244.

Jung JY, Ried MK, Hothorn M, Poirier Y. Control of plant phosphate homeostasis by inositol pyrophosphates and the SPX domain. Curr Opin Biotechnol. 2018; 49: 156-162.

Kobayashi K, Awai K, Nakamura M, Nagatani A, Masuda T, Ohta H. Type-B monogalactosyldiacylglycerol synthases are involved in phosphate starvation-induced lipid remodeling, and are crucial for low-phosphate adaptation. Plant J. 2009; 57: 322-331.

Kuo HF, Chiou TJ. The role of microRNAs in phosphorus deficiency signaling. Plant Physiol. 2011; 156: 1016-1024.

Lan P, Li W, Schmidt W. Complementary proteome and transcriptome profiling in phosphate-deficient Arabidopsis roots reveals multiple levels of gene regulation. Mol Cell Proteomics. 2012; 11: 1156-1166.

Lan P, Li W, Schmidt W. ‘Omics’ approaches towards understanding plant phosphorus acquisition and use. Ann Plant Rev. 2015; 48: 65-97.

Lapis-Gaza HR, Jost R, Finnegan PM. Arabidopsis PHOSPHATE TRANSPORTER1 genes PHT1; 8 and PHT1; 9 are involved in root-to-shoot translocation of orthophosphate. BMC Plant Biol. 2014; 14: 334.

Li X, Zhang Y, Yin L, Lu J. Overexpression of pathogen-induced grapevine TIR-NB-LRR gene VaRGA1 enhances disease resistance and drought and salt tolerance in Nicotiana benthamiana. Protoplasma. 2017; 254: 957-969.

Liu TY, Huang TK, Yang SY, Hong YT, Huang SM, Wang FN, Chiang SF, Tsai SY, Lu WC, Chiou TJ. Identification of plant vacuolar transporters mediating phosphate storage. Nat Commun. 2016; 7: 11095.

Lv Q, Zhong Y, Wang Y, Wang Z, Zhang L, Shi J, Wu Z, Liu Y, Mao C, Yi K, Wu P. SPX4 negatively regulates phosphate signaling and homeostasis through its interaction with PHR2 in rice. Plant Cell. 2014; 26: 1586-1597.

Mehra P, Giri J. Rice and chickpea GDPDs are preferentially influenced by low phosphate and CaGDPD1 encodes an active glycerophosphodiester phosphodiesterase enzyme. Plant Cell Rep. 2016; 35: 1699-1717.

Mehra P, Pandey BK, Verma L, Giri J. A novel glycerophosphodiester phosphodiesterase improves phosphate deficiency tolerance in rice. Plant Cell Environ. 2019; 42: 1167-1179.

Miao J, Sun J, Liu D, Li B, Zhang A, Li Z, Tong Y. Characterization of the promoter of phosphate transporter TaPHT1.2 differentially expressed in wheat varieties. J Genet Genomics. 2009; 36: 455-466.

Mortazavi A, Williams BA, Mccue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008; 5:621- 628.

Muchhal US, Pardo JM, Raghothama K. Phosphate transporters from the higher plant Arabidopsis thaliana. Proc Natl Acad Sci USA. 1996; 93: 10519-10523.

Nakamura Y. Phosphate starvation and membrane lipid remodeling in seed plants. Prog Lipid Res. 2013; 52: 43-50.

O′Rourke JA, Yang SS, Miller SS, Bucciarelli B, Liu J, Rydeen A, Bozsoki Z, Uhde-Stone C, Tu ZJ, Allan D, Gronwald JW, Vance CP. An RNASeq transcriptome analysis of orthophosphate-deficient white lupin reveals novel insights into phosphorus acclimation in plants. Plant Physiol. 2013; 161: 705–724.

Osorio MB, Ng S, Berkowitz O, De Clercq I, Mao C, Shou H, Whelan J, Jost R. SPX4 acts on PHR1-dependent and independent regulation of shoot phosphorus status in Arabidopsis. Plant Physiol. 2019; 181: 332-352.

Ouyang X, Hong X, Zhao X, Zhang W, He X, Ma W, Teng W, Tong Y. Knock out of the PHOSPHATE 2 gene TaPHO2-A1 improves phosphorus uptake and grain yield under low phosphorus conditions in common wheat. Sci Rep. 2016; 6: 29850.

Park BS, Seo JS, Chua NH. NITROGEN LIMITATION ADAPTATION recruits PHOSPHATE2 to target the phosphate transporter PT2 for degradation during the regulation of Arabidopsis phosphate homeostasis. Plant Cell. 2014; 26: 454-464.

Park MR, Baek SH, De los Reyes BG, Yun SJ, Hasenstein KH.Transcriptome profiling characterizes phosphate deficiency effects on carbohydrate metabolism in rice leaves. J Plant Physiol. 2012; 169: 193-205.

Peñaloza E, Muñoz G, Salvo-Garrido H, Silva H, Corcuera LJ. Phosphate deficiency regulates phosphoenolpyruvate carboxylase expression in proteoid root clusters of white lupin. J Exp Bot. 2005; 56: 145-153.

Qin L, Zhao J, Tian J, Chen L, Sun Z, Guo Y, Lu X, Gu M, Xu G, Liao H. The high-affinity phosphate transporter GmPT5 regulates phosphate transport to nodules and nodulation in soybean. Plant Physiol. 2012; 159: 1634-1643.

Ren P, Meng Y, Li B, Ma X, Si E, Lai Y, Wang J, Yao L, Yang K, Shang X, Wang H. Molecular mechanisms of acclimatization to phosphorus starvation and recovery underlying full-length transcriptome profiling in barley (Hordeum vulgare L.). Front Plant Sci. 2018; 9: 500.

Ren Y, Sun F, Hou J, Chen L, Zhang Y, Kang X, Wang Y. Differential profiling analysis of miRNAs reveals a regulatory role in low N stress response of Populus. Funct Integr Genomics. 2015; 15: 93-105.

Ruan W, Guo M, Wang X, Guo Z, Xu Z, Xu L, Zhao H, Sun H, Yan C, Yi K. Two RING-finger ubiquitin E3 ligases regulate the degradation of SPX4, an internal phosphate sensor, for phosphate homeostasis and signaling in rice. Mol Plant. 2019; 12: 1060-1074.

Schroeder JI, Delhaize E, Frommer WB, Guerinot ML, Harrison MJ, Herrera-Estrella L, Horie T, Kochian LV, Munns R, Nishizawa NK, Tsay Y, Sanders D. Using membrane transporters to improve crops for sustainable food production. Nature. 2013; 497: 60-66.

Secco D, Jabnoune M, Walker H, Shou H, Wu P, Poirier Y, Whelan J. Spatio-temporal transcript profiling of rice roots and shoots in response to phosphate starvation and recovery. Plant Cell. 2013; 25: 4285-4303.

Secco D, Wang C, Arpat BA, Wang Z, Poirier Y, Tyerman SD, Wu P, Shou H, Whelan J. The emerging importance of the SPX domain-containing proteins in phosphate homeostasis. New Phytol. 2012; 193: 842-851.

Shemi A, Schatz D, Fredricks HF, Van Mooy BA, Porat Z, Vardi A. Phosphorus starvation induces membrane remodeling and recycling in Emiliania huxleyi. New Phytol. 2016; 211: 886-898.

Shen J, Yan L, Zhang J, Li H, Bai Z, Chen X, Zhang W, Zhang F. Phosphorus dynamics: from soil to plant. Plant Physiol. 2011; 156: 997-1005.

Shin H, Shin HS, Dewbre GR, Harrison MJ. Phosphate transport in Arabidopsis: Pht1; 1 and Pht1; 4 play a major role in phosphate acquisition from both low- and high- phosphate environments. Plant J. 2004; 39: 629-642.

Song H, Cai Z, Liao J, Tang D, Zhang S. Sexually differential gene expressions in poplar roots in response to nitrogen deficiency. Tree Physiol. 2019; 39: 1614-1629.

Strock CF, De La Riva LM, Lynch JP. Reduction in root secondary growth as a strategy for phosphorus acquisition. Plant Physiol. 2018; 176: 691-703.

Tawaraya K, Honda S, Cheng W, Chuba M, Okazaki Y, Saito K, Oikawa A, Maruyama H, Wasaki J, Wagatsuma T. Ancient rice cultivar extensively replaces phospholipids with non-phosphorus glycolipid under phosphorus deficiency. Physiol Plantarum. 2018; 163: 297-305.

Tong Y, Zhang W, Wang X, Couture RM, Larssen T, Zhao Y, Li J, Liang H, Liu X, Bu X, He W, Zhang Q, Lin Y. Decline in Chinese lake phosphorus concentration accompanied by shift in source since 2006. Nature Geosci. 2017; 10: 507-511.

Vance CP, Uhde-Stone C, Allan DLP. Phosphorus acquisition and use: critical adaptations by plants for securing a nonrenewable resource. New Phytol. 2003; 157: 423-447

Veneklaas EJ, Lambers H, Bragg J, Finnegan PM, Lovelock CE, Plaxton WC, Price CA, Scheible WR, Shane MW, White PJ, Raven JA. Opportunities for improving phosphorus-use efficiency in crop plants. New Phytol. 2012; 195: 306-320.

Wild R, Gerasimaite R, Jung JY, Truffault V, Pavlovic I, Schmidt A, Saiardi A, Jessen HJ, Poirier Y, Hothorn M, Mayer A. Control of eukaryotic phosphate homeostasis by inositol polyphosphate sensor domains. Science. 2016; 352: 986-990.

Xie R, Zhang J, Ma Y, Pan X, Dong C, Pang S, He S, Deng L, Yi S, Zheng Y, Lv Q. Combined analysis of mRNA and miRNA identifies dehydration and salinity responsive key molecular players in citrus roots. Sci Rep. 2017; 7: 42094.

Zhang Y, Zhou Z, Ma X, Jin G. Foraging ability and growth performance of four subtropical tree species in response to heterogeneous nutrient environments. J For Res. 2010; 15: 91-98.

Zhang Z, Liao H, Lucas WJ. Molecular mechanisms underlying phosphate sensing, signaling, and adaptation in plants. J Integr Plant Biol. 2014; 56: 192-220.

Zhang Z, Wang Z, Wang X, Liu D. Blue light-triggered chemical reactions underlie phosphate deficiency-induced inhibition of root elongation of Arabidopsis seedlings grown in Petri dishes. Mol Plant. 2019; https://doi.org/10.1016/j.molp.2019.08.001.

Zhao M, Ding H, Zhu JK, Zhang FS, Li WX. Involvement of miR169 in the nitrogen-starvation response in Arabidopsis. New Phytol. 2011; 190: 906-915.

Zhao P, Wang L, Yin H. Transcriptional responses to phosphate starvation in Brachypodium distachyon roots. Plant Physiol Bioch. 2018; 122: 113-120.

Table 1. List of the DEGs related to phosphorus metabolism after 48 d of Pi-deficiency in the roots and shoots

Gene ID	Description	Fold Change: Log2(Pi/CK)
Gene ID	Description	Root	Shoot
Transporter
Unigene0012504	Phosphate transporter, PHT1-5		2.71
Unigene0002389	Phosphate transporter, PHT1-4		2.57
Unigene0072510	Phosphate transporter, PHT1-12		2.54
Unigene0072512	Phosphate transporter, PHT1-12	1.98	2.05
Unigene0066726	Phosphate transporter, PHT1-4		1.73
Unigene0066725	Phosphate transporter, PHT1-4		1.56
Unigene0075955	Phosphate transporter, PHT1-4	1.36	1.51
Unigene0007376	Phosphate transporter, PHT1-5	1.13
Unigene0075954	Phosphate transporter, PHT1-1		1.03
Unigene0075150	Nitrate transporter		1.86
Transcription
Unigene0025956	SPX domain-containing protein 1	2.70	2.80
Unigene0025957	SPX domain-containing protein 1	2.70	2.70
Unigene0055609	SPX domain-containing protein 2	1.19	1.73
Unigene0004869	SPX domain-containing protein 4	1.07	1.02
Unigene0027298	SPX, N-terminal	-1.22	-1.35
Unigene0011607	Myb family transcription factor	2.68
Unigene0071681	Myb-related protein 308	1.52
Unigene0006787	Transcription factor MYB30		2.36
Unigene0067897	MYB-like protein		1.81
Unigene0018015	Ethylene responsive factor 1-like protein	2.87
Unigene0080776	Ethylene responsive factor 1-like protein		2.75
Unigene0071017	Ethylene responsive factor 1-like protein		2.38
Phosphorylation/Dephosphorylation
Unigene0061850	Purple acid phosphatase 1	6.14
Unigene0047568	Purple acid phosphatase-like protein	2.11	5.52
Unigene0069941	Purple acid phosphatase 8-like		1.65
Unigene0056388	Purple acid phosphatase 2-like isoform X2	1.57	1.63
Unigene0074941	Acid phosphatase, partial	1.56	2.17
Unigene0017024	Acid phosphatase, partial	1.51	2.10
Unigene0055606	Purple acid phosphatase 23 isoform X1		1.38
Unigene0066113	Purple acid phosphatase 3-like		1.24
Unigene0055621	Purple acid phosphatase 15-like	1.10
Unigene0055622	Purple acid phosphatase 23 isoform X2		1.10
Unigene0010082	Acid phosphatase 1		1.08
Unigene0056389	Purple acid phosphatase 2		1.02
Lipid metabolism
Unigene0068798	Glycerophosphodiester phosphodiesterase, GDPD1	2.08	2.21
Unigene0068803	Glycerophosphodiester phosphodiesterase, GDPD2	2.07	2.33
Unigene0022807	Monogalactosyldiacylglycerol synthase, MGD1		1.75
Unigene0055927	UDP-sulfoquinovose synthase, SQD1	1.33	1.02
Metabolism
Unigene0016253	Phosphoenolpyruvate carboxylase	1.80	1.34
Unigene0036016	Phosphoenolpyruvate carboxylase	1.35
Unigene0070466	Phosphoenolpyruvate carboxylase		1.45
Unigene0079423	Haloacid dehalogenase-like hydrolase	1.55
Miscellaneous
Unigene0075158	Disease resistance protein	3.46
Unigene0073312	Disease resistance protein	3.43
Unigene0014755	Disease resistance protein	3.00
Unigene0096400	Disease resistance protein		8.48
Unigene0002124	Disease resistance protein		8.40
Unigene0071934	Disease resistance protein		5.22
Unigene0067522	U-box domain-containing protein kinase		6.28
Unigene0071447	U-box domain-containing protein kinase		2.42
Unigene0068986	E3 ubiquitin-protein ligase RING1	2.21	3.68

Download PDF

Editorial decision: Major revision
10 Jun, 2020
Review #2 received at journal
12 May, 2020
Reviewer #2 agreed at journal
29 Apr, 2020
Review #1 received at journal
26 Mar, 2020
Reviewer #1 agreed at journal
16 Mar, 2020
Reviewers invited by journal
15 Mar, 2020
Editor assigned by journal
12 Mar, 2020
Submission checks completed at journal
06 Mar, 2020
Editor invited by journal
06 Mar, 2020

You are reading this latest preprint version

Integrated mRNA and miRNA expression profiling analyses of Pinus massoniana roots and shoots in response to phosphate deficiency

Status:

Version 1

Abstract

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Table

Supplementary Files

Status:

Version 1