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

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. The high-affinity phosphate transporter GmPT5 regulates phosphate transport to nodules and nodulation in

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 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 ubiquitinconjugating 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), UDPsulfoquinovose 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-Pitreated 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 lowquality raw reads and removing rRNAs. On average, each sample had 84.0 million highquality 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 Pisufficient 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 Pideficient 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 Timeseries Expression Miner (STEM) with default parameters. Eight significant gene profiles (Pvalue ≤ 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 Pideficient 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 2 (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.

Discussion
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,  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 . Based on this concept, Pi transporters (PHT1 family) can work as external Pi sensors . 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 . Sucrose and starch have been reported to be increased under Pi starvation in plants (Park et al., 2012;. 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  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 2 (AP 2 ) domain-containing transcription factors, ERFs involve in ethylenemediated 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 Pistarvation-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 nitrogenstarvation 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 . 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.

Conclusions
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 2 PO 4 ) and sufficient in Pi (0.5 mM KH 2 PO 4 ) 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 2 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 6 . The DE miRNAs between the deficient-Pi and sufficient-Pi (control) treatments were identified via the formula log 2 (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)  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  Figure 1 Morphological and physiological responses of Masson pine seedlings to Pi deficiency, and length distribution of transcripts using next-generation sequencing.

Figure 2
Differentially expressed profiling of the transcriptome response to Pi deficiency.  KEGG enrichment analysis diagram of DEGs. Q values correspond to P-values corrected for multiple hypothesis testing for KEGG analyses, with values from 0 -1, with 0 indicating a more significant enrichment than 1.