Pi-starvation induced transcriptional changes in barley revealed by a comprehensive RNA-Seq and degradome analyses

Background Small RNAs (sRNAs) are 20–30 nt regulatory elements which are responsible for plant development regulation and participate in many plant stress responses. Insufficient inorganic phosphate (Pi) concentration triggers plant responses to balance the internal Pi level. Results In this study, we describe Pi-starvation-responsive small RNAs and transcriptome changes in barley (Hordeum vulgare L.) using Next-Generation Sequencing (NGS) RNA-Seq data derived from three different types of NGS libraries: (i) small RNAs, (ii) degraded RNAs, and (iii) functional mRNAs. We find that differentially and significantly expressed miRNAs (DEMs, Bonferroni adjusted p-value < 0.05) are represented by 15 molecules in shoot and 13 in root; mainly various miR399 and miR827 isomiRs. The remaining small RNAs (i.e., those without perfect match to reference sequences deposited in miRBase) are considered as differentially expressed other sRNAs (DESs, p-value Bonferroni correction < 0.05). In roots, a more abundant and diverse set of other sRNAs (DESs, 1796 unique sequences, 0.13% from the average of the unique small RNA expressed under low-Pi) contributes more to the compensation of low-Pi stress than that in shoots (DESs, 199 unique sequences, 0.01%). More than 80% of differentially expressed other sRNAs are up-regulated in both organs. Additionally, in barley shoots, up-regulation of small RNAs is accompanied by strong induction of two nucleases (S1/P1 endonuclease and 3′-5′ exonuclease). This suggests that most small RNAs may be generated upon nucleolytic cleavage to increase the internal Pi pool. Transcriptomic profiling of Pi-starved barley shoots identifies 98 differentially expressed genes (DEGs). A majority of the DEGs possess characteristic Pi-responsive cis-regulatory elements (P1BS and/or PHO element), located mostly in the proximal promoter regions. GO analysis shows that the discovered DEGs primarily alter plant defense, plant stress response, nutrient mobilization, or pathways involved in the gathering and recycling of phosphorus from organic pools. Conclusions Our results provide comprehensive data to demonstrate complex responses at the RNA level in barley to maintain Pi homeostasis and indicate that barley adapts to Pi-starvation through elicitation of RNA degradation. Novel P-responsive genes were selected as putative candidates to overcome low-Pi stress in barley plants. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07481-w.

tolerance to heat stress [32], (ii) confer drought tolerance [33], (iii) regulate low-potassium tolerance [34], (iv) respond to aluminum stress [35], and (v) maintain inorganic phosphate (Pi) homeostasis [36]. On the other hand, siRNAs mostly function as a defenders of genome integrity in response to foreign nucleic acids [37]. The TAS3 gene expresses ta-siRNAs, which may negatively regulate auxin signaling by targeting AUXIN RESPONSE FACTOR 3 (ARF3) transcripts [38] and moderate oral architecture in response to drought stress in Arabidopsis thaliana L. [39]. The TAS-ARF pathway has been shown to be involved either in the development process of maize (Zea mays L.) [40] or regulating lateral root growth in Arabidopsis [41]. In addition, tRNA-derived small RNAs have been shown to accumulate in Arabidopsis roots under Pi-starvation [42], while rhizobial tRFs can regulate nodule formation in soybean (Glycine max L.) [13].
Changes in soil nutrient concentrations lead to aberrations in the set of sRNAs, with respect to the prevailing severe environmental conditions [6]. One of the most important macronutrients, which is indispensable for proper plant growth, is phosphorus (P) [43,44]. P is a component of DNA, RNA, phospholipids, and ATP, and is involved in several biochemical processes such as protein phosphorylation, energy storage and transfer, and regulation of protein synthesis [45]. From soil matrices, P is acquired by the root system in the form of inorganic phosphate (Pi) ions. Insu cient Pi supply leads to barley growth inhibition [46,47]. Plant transcriptome response to Pi-starvation involves protein coding genes, sRNAs, and long non-coding RNAs that form regulatory feedback loops. The most widely studied molecules in this context-miRNA399 molecules-are up-regulated in barley shoots and roots under low-Pi conditions [36]. MiRNA399 targets the 5'-UTR of the barley PHO2 (PHOSPHATE 2) transcripts [48], encoding an ubiquitin-conjugating E2 enzyme (UBC24), a negative regulator of Pi uptake and root-toshoot translocation. PHO2 is involved in ubiquitination of PHOSPHATE TRANSPORTER 1 (PHT1) family [49] and PHOSPHATE TRANSPORTER TRAFFIC FACILITATOR 1 (PHF1) [49]. Transgenic Arabidopsis thaliana plants overexpressing miR399 accumulate excessive Pi in shoots and display Pi overaccumulation toxic symptoms. Likewise, such a phenotype has been reported for the pho2 loss-offunction Arabidopsis mutant [50,51]. Thus, plants have developed a strategy to regulate the level of miR399 in the cytoplasm. The non-coding RNA molecule, IPS1 (INDUCED BY PHOSPHATE STARVATION 1), has been shown to be highly expressed in plants exposed to Pi-starvation [52][53][54]. IPS1 is a noncleavable miR399 target which inhibits miR399-mediated down-regulation of PHO2 mRNA by target mimicry [54]. Thus, the RNAi effect of miRNA activity may be counterbalanced by other RNAs, in a stressdependent manner.
Deep sequencing of sRNAs has uncovered up-regulation of miRNAs like miR156, miR778, miR827, and miR2111, and down-regulation of miR169, miR395, and miR398 in Arabidopsis plants upon Pi deprivation [42,55]. In rice (Oryza sativa L.), Pi-starvation induced the expression level of miR827 molecules, which dysregulate the transcript level of two genes encoding the SPX-MFS (named after proteins SYG¹⁄PHO ¹⁄XPR1 and the protein domain Major Facility Superfamily) protein family members SPX-MFS1 and SPX-MFS2 [56,57]. These two SPX-MFS membrane transporters mediate Pi transport and control Pi homeostasis in shoot [58]. In Arabidopsis, the level of mature miR778 was up-regulated in shoots and roots in low-Pi conditions, while its target gene expression SUVH6 (SU(VAR)3-9 HOMOLOG 6) was accordingly reduced [59]. The SUVH6 gene encodes a histone H3 lysine 9 (H3K9) methyltransferase, which may enable plants to adapt to environmental conditions by changing their chromatin structure [60]. miR2111 functions as an activator of rhizobial nodulation, which is strictly correlated with the balanced assimilation of nitrogen (N) and P in plants [61,62]. However, there is still a gap in understanding how Pistarvation affects the quantity and quality of sRNAs distributed in barley shoots and roots. What kind of sRNAs are preferentially induced? What is the role of sRNAs in responding to Pi-starvation? What are the mRNA targets recognized by those sRNAs in barley?
In this paper, we analyzed changes in the expression levels of RNAs in barley growing under Pi-starvation, as compared to control/Pi su cient conditions. Our results support the hypothesis that Pi-starvation triggers underlying molecular mechanisms and the expression level of key genes involved in maintaining proper barley growth and development. Combined deep sequencing data (sRNAs, degradome and mRNAs) reveals the widespread importance of low-Pi-dependent miRNAs and genes representing various biological pathways. Using degradome analysis, we identi ed mRNAs targeted by sRNAs identi ed in this study. Among these sRNAs, only a small fraction maps perfectly to miRNA sequences deposited in miRBase. Our degradome data show that most sRNAs produced upon Pi-starvation are not involved in gene silencing. In addition, we performed transcriptome analysis of the protein-coding gene expression in barley shoots upon Pi-starvation. Subsequent analyses were performed (GO analysis, chromosomal mapping, and Pi-responsive motifs localization) to characterize speci c stress responses in barley plants to accomplish Pi homeostasis.

Results
Barley plants display low-Pi symptoms at the morphological and molecular levels Severe low-Pi responses were induced in the barley plant line Rolap grown in the soil containing 8 mg P/kg. P undernourishment caused over 2-fold reduction of plant shoot biomass (Fig. 1A). Shoot fresh weight of plants at 23 rd day post-sowing (dps) was signi cantly reduced, in comparison with control plants, with average mass 8.8 g for stressed plants and 18.5 g for plants growing under Pi-su cient conditions (p = 0.001) (Fig. 1B). We observed a signi cantly decreased concentration of Pi ions, with only 0.48 µmol Pi per g of fresh root weight (FW) and 4.2 µmol Pi per g of shoot FW, when compared with the control plants having 3.84 (p = 0.0056) and 24.35 µmol Pi/g FW (p = 0.0001), respectively (Fig. 1C). To examine the induction of changes at a molecular level by low-Pi stress in barley plants, we measured the absolute gene expression of the low-Pi-responsive marker gene IPS1. The barley IPS1 gene is highly expressed under Pi-de cient conditions in the plant line Rolap. At the tillering stage (23 dps), we detected 4191 copies of IPS1 RNA for low-Pi treated roots, normalized per 1000 copies of ADP-RIBOSYLATION FACTOR 1-LIKE (ARF1) reference gene, in comparison to the control plants, with only 58 copies of IPS1 RNA (p = 0.00006) (Additional le 1). Taking validated plant material, we performed tripartite deepsequencing analysis to: (i) identify Pi-responsive sRNAs, (ii) elucidate changes in the barley transcriptome upon Pi starvation, and (iii) identify mRNA targets for Pi-responsive sRNAs through degradome sequencing (Fig. 2).
Barley plants express an organ-speci c set of microRNAs in response to low-Pi conditions We performed small RNA deep-sequencing to nd out which small RNAs are up-or down-regulated by Pi starvation in barley shoots and roots. The average of 30.4 mln reads for roots and 25.2 mln reads for shoots were generated in 50 nt single-read Illumina sequencing (Additional le 2). After adapter and quality trimming, we mapped reads to the miRBase Sequence Database (release 22) to annotate miRNAderived sequences [63]. A set of parameters were used to de ne the pool of differentially expressed miRNAs: (i) no mismatches with the reference sequences in the miRBase were allowed; (ii) different types of miRNA sequences were permitted, whether they were annotated as precursor, mature, or isomiR; (iii) miRNA sequences were named accordingly to the name of the assigned reference miRNA; and (iv) signi cance of fold change (p-value < 0.05) was additionally veri ed using a restricted Bonferroni p-value test (Fig. 2).
We found 138 and 162 differentially expressed miRNAs (DEMs) annotated to the miRBase (p-value < 0.05) in barley roots and shoots, respectively. Only 25 DEMs were expressed in both examined barley organs (Additional le 3). However, restricted Bonferroni p-value correction narrowed down set of miRNAs to 15 in shoots and 13 in roots (Fig. 3A, left panel). We focused further only on those that passed the Bonferroni test (Fig. 3A, right panel). In both organs, most of the DEMs were signi cantly up-regulated. Interestingly, out of 15 miRNA, only miRNA166d was down-regulated in shoot under low-Pi (log 2 (fold change) = -1.18). In our previous work, we showed that miRNA166 is expressed in barley during different developmental stages reaching the highest level in 2-week-old plants [64]. miRNA166 plays an important role in plant development, including root and leaf patterning, by targeting mRNA encoding HOMEODOMAIN LEUCINE-ZIPPER CLASS III (HD-ZIP III) transcription factors [65]. Similarly, only miRNA319b out of 13 DEMs was down-regulated in low-Pi treated roots (log 2 (fold change) =-1.28). In a previous study, we presented data that Arabidopsis miR319 is a multi-stress responsiveness miRNA [22]. For example, MIR319b gene expression was down-regulated in response to drought, heat, and salinity, but up-regulated in response to copper and sulfur de ciency stresses [22].

Uncharacterized miRNAs are highly induced by low-Pi stress in barley roots
A speci c set of miRNAs was expressed in barley shoot or root under low-Pi (Additional le 3). In shoot, only two miRNA families, miRNA399 and miRNA827, were induced, while in root we observed a more diverse response (Fig. 3A, right panel). Apart from miRNA399/miRNA827 induction, we found the following additional miRNA to be up-regulated in root: two miRNA5083 (id = 3, and id = 4), miRNA1511 (id = 6), two miRNA9779 (id = 16, and id = 17), two miRNA156 (id = 65, and id = 69), and miRNA5072 (id = 118). Among these eight miRNAs, only miR156 has been reported before as Pi-responsive in Arabidopsis [42,55]. The miR156 isomiRs were also found dysregulated in shoot, but none of them pass the Bonferroni test. Our results suggest that there is a more complex response to low-Pi stress regarding miRNA expression in roots than in shoots, where the miRNA action is directed to control the transcript level of either PHO2, SPX-MFS1, or SPX-MFS2 by just two miRNA families.
In order to determine the potential cleavage activity of miRNAs identi ed in root we performed degradome analysis. None of these eight miRNAs were found to recognize a speci c target from the pool of barley mRNAs. Our result raises the hypothesis that plants may release Pi from the pool of stress-responsive small RNAs. In addition, most of them match to the pre-miRNAs, not mature molecules (Fig. 3A, right panel).
NGS (Next Generation Sequencing) data were validated by complex analysis of mature miR827 (from 3' arm) in all samples taken for deep sequencing. The absolute expression level of miR827 is signi cantly up-regulated in both shoots and roots under a low-Pi regime (Fig. 3B). The log 2 fold change of miR827 molecules de ned by deep-sequencing in shoot (id: 2073) was found on the same level in root (id: 114), log 2 (fc) = 3.05 and 3.01, respectively (Fig. 3A). The ddPCR results were consistent with NGS data showing up-regulation of mature miR827 molecule in both tested organs (Fig. 3B). These data were con rmed by northern blot hybridization (Fig. 3C) and degradome analysis (Fig. 3D). The results showed that SPX-MFS1 transcripts have been recognized and cleaved by Ago protein associated with miR827 in barley (p = 0.014) (Fig. 3D).

Different classes of small RNAs in barley accumulate in an organ-speci c manner under low-Pi regime
The small RNA reads which did not map to miRBase were mapped to particular classes of barley cDNAs derived from the Ensembl Plants database (release 40), both separately and to all of them (Fig. 2). All sequences mapped to barley cDNAs are listed in Additional le 3. We found that small RNAs, other than We analyzed whether different lengths (taking sequences from 18 to 25 nt in lenght) and classes of small RNAs contributed to either root or shoot response to low-Pi conditions. In roots, the length distribution of DESs remained balanced, from 10.91 % for the representation of 24 nt sequences to 15.26 % for the 18 nt sequences, which were the most abundant (including 274 DESs) (Fig. 4B). In shoots, the representations of DES lengths uctuated more than in roots. The 19 nt sequences were the most visible (21.11 %), while three representations did not score more than 10%: the 22 nt . We did not nd any DESs annotated to the snRNAs, SRP-RNAs, or tRNAs from barley shoot upon low-Pi. In addition, total numbers of 166 DESs (83 %) in shoots and 1560 DESs (87 %) in roots were signi cantly upregulated after exposure to low-Pi stress (Additional le 4, Additional le 5, Fig. 7).
Among the unannotated reads of sRNAs in roots, the highest fold change was observed for a 19 nt The results obtained in this study show again that barley roots exhibit a more diverse pool of Piresponsive small RNAs which may trigger developmental adaptation of the root to Pi-starvation. Additionally, 613 rRNA-derived sRNAs are up-regulated, whereas 176 rRNA-derived sRNAs are downregulated in barley roots (Additional le 5). We believe that such sRNA may be further processed, serving as a Pi source to compensate Pi de ciency.

Identi cation of barley genes responsive to Pi-starvation
Since we observed, that most of the other sRNAs in shoot were derived from either protein-coding mRNAs or non-translating RNAs, we checked whether this observation is correlated with gene expression changes of polyadenylated RNAs in barley shoot under Pi-starvation. Among 98 of identi ed DEGs, the transcripts of 56 annotated loci were signi cantly up-regulated, while those derived from 42 loci were down-regulated in Pi-starved barley shoots (Additional le 8). Repressed loci were found to be preferentially located at barley chromosome no. 2, while induced loci were found mostly at barley chromosomes no. 3 and 5 (Fig.  5B).
Absolute quanti cation of a few selected transcripts was performed to validate RNA-Seq data obtained in this study. Two genes which were highly induced (encoding endonuclease S1/P1 and 3'-5'-exonuclease) and two which were severely repressed (encoding oxalate oxidases) under the low-Pi regime were taken for ddPCR (droplet digital PCR) analysis (Fig. 6A). We con rmed statistically signi cant changes (p < 0.05) in normalized copy number (per 1000 copies of the ARF1 reference gene) of all genes taken for analysis.

Pi-responsive motifs found in the promoters of DEGs
In general, genes that are affected by Pi status possess characteristic cis-regulatory elements within either promoter or 5'-UTR regions [81]. Previously, we have shown the importance of the P1BS motif (PHR1 binding sequence, consensus GnATATnC, [82]) and P-responsive PHO elements (consensus ATGCCAT, [83]) in the expression e ciency of the barley PHO2 gene [48]. Both motifs may bind PHR-like (PHOSPHATE STARVATION RESPONSE) transcription factors (TFs) and act as activators or repressors of downstream gene expression in a Pi-dependent manner [84]. Likewise, we hypothesized that regulatory regions of the identi ed DEGs had Pi-responsive motifs, which may be bound by PHR TFs, causing gene expression dysregulation. To con rm this hypothesis, we analyzed DNA sequences from the 2000 bp region upstream of the predicted transcription start sites from all 98 DEGs (Additional le 9). In the next step, promoter data were directly screened for P1BS and P-responsive PHO element consensus sequences by multiple promoter analysis using the PlantPAN3.0 tool [85]. We con rmed the presence of Pi-dependent motif in 55 out of 98 DEGs promoters. An in silico approach detected 46 DEGs having at least one P1BS motif (Additional le 10) and 17 DEGs with at least one P-responsive PHO element (Fig.  6B, Additional le 11). The most over-represented motifs were found in the promoters of genes encoding sulfoquinovosyl transferase SQD2-like (log 2 (fc) = 3.74) [86], phosphoenolpyruvate carboxylase 1-like (log 2 (fc) = 1.73) [87], and pyridoxal phosphate-dependent transferase (log 2 (fc) = -2.61) [88]. Each of the genes harbor three P1BSs and one P-responsive PHO element, as well (Additional le 8).
Degradome pro ling describes post-transcriptional regulatory network of identi ed DEMs After identi cation of (i) differentially expressed miRNAs (DEMs), (ii) other sRNAs (DESs), and (iii) mRNAs (DEGs), we used this comprehensive data together with cDNAs annotated in the Ensembl Plants database to identify the sRNAs directly involved in RNA degradation. The DESs were also examined, because we assumed that there may have been putative miRNAs that were not mapped to the miRbase, due to restricted query settings allowing no mismatch. Molecules which exhibited a single mismatch (or more) may still function as miRNA in barley. Degradome libraries were carried out for root, as well as for shoot, and sequenced using an Illumina System. The received data were analyzed using two independent in silico approaches (Fig. 2). At times, the different algorithms used elicited different miRNA targets; however, the general degradome pattern was equivalent for both approaches (Fig. 3A, right panel, Additional les 12-24).
First, we searched for the potential target mRNAs for a broader set of differentially expressed miRNAs with signi cant fold change (p-value without Bonferroni correction), taking 162 DEMs from shoot and 138 DEMs from root, respectively (Additional le 3). A total of 1964 scores were obtained for shoot DEMs (1303 using the TargetSeek approach and 661 using PAREsnip2) (Additional le 12, Additional le 16), while in root there were 1515 records (921 and 594, respectively) (Additional le 14, Additional le 20, Additional le 21). Accordingly, shoot examination proved the proper selection of miRNAs from all differentially expressed small RNAs, as all 661 records predicted by the PAREsnip2 approach were categorized with the best score (0). A majority of records corresponded to different miR399 and/or miR827 isomiRs and their known targets PHO2 or SPX-MFS1/SPX-MFS2, respectively. In shoot, the best scoring miRNA:mRNA match was found for mature miRNA399d (21 nt, id = 2066), which guides cleavage within the 5'-UTR of PHO2 mRNA (isoform no. 3) in position 857 (TargetSeek: score = 0, MFE = -41,1, PAREsnip2: score= 0, MFE = -41.1) ( Table 1). The miR5049e (21 nt) targets various isoforms of mRNA encoding ENT domain-containing proteins with high probability (TargetSeek: score = 1.5, MFE = -27.4), which may be involved in the plant defense response to fungal infection and anthocyanin biosynthesis [89,90]. Moreover, we found two target mRNAs (HORVU3Hr1G094730, HORVU6Hr1G031450) encoding SPL-like TFs guided for cleavage by two miR156 isomiRs (id: 2053 and 2054) in barley shoot (Additional le 16, Additional le 17). The SPL genes are post-transcriptionally regulated by miR156 and control shoot branching in rice [91].
Putative regulatory small RNAs identi ed in degradome data Degradome pro ling was performed to test whether any of the sequences from the 1796 DESs found in roots or 199 DESs found in shoots contribute to the complexity of gene regulation during low-Pi stress. A total of 759 records (245 using the TargetSeek approach and 514 using PAREsnip2) were found in the degradome pro les matching root DESs (Additional le 15, Additional le 22, Additional le 23) and 160 records (87 and 73, respectively) matching shoot DESs (Additional le 13, Additional le 18, Additional le 19). Taking only either the most up-regulated or the most down-regulated sRNAs for degradome screening, we found six promising target genes in shoot and ve in root (Table 2). For example, in roots, the highly up-regulated 20 nt sequence 5'-ACCGACCTACTTGACCCTTC-3' (log 2 (fc) = 6.46, id = 348) binds to the 3'-UTR region of the MYB44 TF's mRNA and guides/promotes cleavage in the 1037 position (PAREsnip2: score = 4; MFE = -33.3) ( Table 2). RNA-Seq data for potato (Solanum tuberosum L.) proved that expression of the MYB44 gene is highly downregulated under low-Pi in roots [98], which may be the result of miRNA-guided PTGS. Studies in potato have indicated that MYB44 TF may form a regulatory complex together with WRKY6 TF, which negatively regulates Pi transport by suppressing PHO1 expression [98]. Other degradome records, among the most differentially expressed sRNAs, were found to target mRNAs of the V-ATPase assembly factor (VMA21-like) and three barley genomic loci encoding uncharacterized proteins (HORVU7Hr1G053570, HORVU1Hr1G027340, and HORVU0Hr1G023910) ( Table   2). For example, the potential cleavage activity was predicted for 24 nt sequence 5'-AGAGGAAACTCTGGTGGAGGCTCG-3' (log 2 (fc) = -3.58, id = 463), which may target the mRNA encoding uncharacterized protein with unknown PTHR47188 domain (Fig. 4C).
Some other interesting Pi-related targets which are recognized by DESs were found in our root degradome data, but the prediction scores were weaker than those in the examples described above. For instance: nitrate reductase (HORVU6Hr1G003300), high-a nity nitrate transporter-activating protein (HORVU5Hr1G115500), MYB-like TF (HORVU7Hr1G027370), and stress-induced TF NAC1 (HORVU5Hr1G111590) were found. Interestingly, among the 98 DEGs identi ed in this study, only two of them (SPX-MFS1 and SPX-MFS2) were found as putative targets of miRNA guided activity. In addition, none of the DEGs gene IDs were found to match with any of the identi ed IDs classi ed for differentially expressed small RNAs. All predicted DESs targets are listed in Additional les 14, 16, 19, and 23.

Discussion
In this study, we used a tripartite approach (sRNA-Seq, mRNA RNA-Seq, and degradome-seq) to describe the set of small RNAs differentially expressed in barley roots and shoots under low-Pi stress. We detailed the sophisticated responses of barley shoots and roots involved in the maintaining of Pi homeostasis (Fig. 6). Integrated deep-sequencing data were used to describe organ-speci c adaptations to low-Pi through either activation or repression of different classes of 18-25 nt small RNAs. Additionally, the mRNA-Seq analysis of low-Pi treated barley shoot was performed to analyze the correlation between shoot-derived small RNAs, annotated to either protein-coding mRNAs (47.87%) or non-translating RNAs (36.49%), and gene expression changes of polyadenylated RNAs. We identi ed a total of 28 differentially expressed miRNAs (Bonferroni tested p-value) annotated to miRBase (release 22) without mismatches and a total of 1989 differentially expressed other small RNAs (Bonferroni tested p-value).
In plants, a limited number of miRNA have been shown to be speci cally and strongly induced by Pi limitation, including miRNA399 [102], miRNA778 [59], miRNA827 [55], and miRNA2111 [55,103]. In this work, the majority of DEMs represent various miR399 and miR827 isomiRs in both tested organs. Our results are consistent with sRNA sequencing data published for Arabidopsis [42,55] and Nicotiana benthamiana L. [93]. In both plant species, authors have shown that the number of various miR399 isomiRs was the most abundant in shoots and roots under low-Pi. Eight of the 15 DEMs (after Bonferroni p-value correction) we found in barley shoots belonged to the miR399 family. However, in root, miR399 was represented only by one DEM; the miR399b isomiR (Fig. 3A, right panel). Previously, our absolute copy number analysis of mature miR399 demonstrated that its normalized expression level is 4-fold downregulated in barley roots, as compared to in shoots, under a low-Pi regime [104]. The long-distance movement of signal molecules is known to be crucial for Pi recycling and allocation from root to shoot.
The root system is responsible for Pi acquisition conducted by phosphate transporters belonging to PHT1 family, which saturate cell membranes during Pi de ciency [105]. The level of PHT1 proteins is negatively controlled by the PHO2, which is suppressed by miR399 (see model in Fig. 7) [106]. A high level of miR399 molecules was detected in Arabidopsis wild type rootstocks grafted with miR399-overexpressing scions [42,107]. Thus, miR399 is involved in a plant's systemic response to low-Pi conditions and acts as a long-distance signal, moving from shoot to root to control Pi homeostasis [107]. In Arabidopsis, miR827 has been shown in multiple studies to target the 5'-UTR of the NITROGEN LIMITATION ADAPTATION (NLA) gene [108,109]. In rice, the OsNLA mRNA has a 'degenerate' osa-miR827 potential cleavage site, that is why miR827 does not cleave the OsNLA transcript in vivo [56,57,110]. Likewise, we did not nd NLA mRNAs to be targeted by any of the identi ed hvu-miR827 isomiRs in our barley degradome records. The NLA gene encodes an E3 ubiquitin-protein ligase with RING and SPX domains, which interacts with the PHO2 to prevent the excessive accumulation of Pi [111]. In roots, a more diverse set of miRNAs contributed to the compensation of low-Pi stress, compared to that in shoots. We found six up-regulated miRNA molecules (DEMs) in roots mapped to pre-miRNAs, such as: miR9779, miR5083, miR9779, miR5083, miR1511 and miR5072. In addition, none of them was found in our degradome analysis. The differentially expressed other small RNAs in roots (DESs, 1796 molecules) were represented by 90 % of the total set of other sRNAs (DESs from both organs), annotated to all classes of cDNAs taken for analysis. Among the identi ed set of DESs, many putative miRNA-like molecules were predicted to target various mRNAs involved in plant adaptations to abiotic stresses, plant defense, and/or transcription (Table 1). Further analysis will be performed to experimentally validate the in silico predicted PTGS role of Pi-responsive small RNAs found in this study, as well.
In this paper, we detailed the particular shoot differentially expressed genes (DEGs) harboring Piresponsive cis-regulatory elements, involving various molecular pathways and biological processes. These DEGs were mostly engaged in Pi mobilization and utilization upon Pi-starvation in barley shoots.
Other sRNAs selected from shoots were much less abundant and represent sequences belonging mostly to non-translating and/or protein-coding mRNAs. None of the sRNAs mapped to the differentially expressed mRNAs found in the transcriptomic analysis, suggesting that they may inhibit gene expression through translational repression or may serve as a Pi source for developing plant organs. Plants are adapted to recycle nutrients from senescing organs. For example, class II RNases are involved in the degradation of housekeeping rRNAs before cell death occurs [112]. During senescence extracellular class I RNases were shown to degrade RNA during Pi-starvation in Arabidopsis as well [113]. In 2018, Ren et al. published RNA-Seq data describing the barley transcriptome under low-Pi stress [114]. The authors compared the transcriptomes of two barley genotypes with contrasting low-Pi stress tolerance. In roots, they observed 28 DEGs classi ed into the following functional groups: Pi transport, transcription, lipid metabolism, metabolism, and phosphorylation/dephosphorylation [114]. Likewise, our RNA-Seq data from barley shoot discovered the DEGs involved in all mentioned functional groups. In our shoot transcriptome analysis, we found the same four DEGs: (i) GLYCEROPHOSPHODIESTER PHOSPHODIESTERASE 1 (GDPD1) gene (HORVU3Hr1G079900, log 2 (fc) = 5.78), (ii) MONOGALACTOSYLDIACYLGLYCEROL SYNTHASE 2 gene (MGD2, HORVU4Hr1G044140, 5.27), (iii) SPX5 (HORVU2Hr1G031400, 3.44), and (iv) SPX1 (HORVU7Hr1G089910, 1.78). Furthermore, Ren et al. found three genes encoding purple acid phosphatases (PAPs) [114]; however, they appeared from different barley genome loci than the ve PAPs we found in our study. It was shown that vacuolar and secreted PAPs are involved in Pi scavenging and remobilization during Pi-starvation and leaf senescence. Other related RNA-Seq data published for either wheat (Triticum aestivum L.) [115], rice [116], soybean [117], Plantago major L. [118], and maize [119] demonstrated similar molecular patterns to those in our study (Additional le 8). Based on our work we propose a model of barley adaptations to Pi-starvation (Fig.7).
Interestingly, the presence of crucial Pi-responsive cis-regulatory elements within the promoter regions of more than 50 % of identi ed DEGs may indicate their essential and direct role in conditioning low-Pi tolerance (Fig. 6B). The most widely studied PHR TFs, such as PHR1 in Arabidopsis [82] and PHR2 in rice [120], bind to P1BS elements present in the promoter of a broad range of Pi-related genes. Moreover, the PHR protein family exhibits high functional redundancy and its protein members may co-operatively form a regulatory network to maintain Pi homeostasis in plants [84]. In our previous paper, we showed that, within the 5'-UTR of the PHO2 gene, there is another Pi-responsive motif called the PHO element in close proximity to the P1BS [48]. The PHO element can be bound by PHR-like transcription factors in barley plants, as well [48], and has been found in the promoters of many DEGs in independent Arabidopsis [121,122] or soybean [117] studies.
The elevated abundance of sRNAs has been associated with the up-regulation of two types of nucleases (endonuclease S1/P1 and 3'-5' exonuclease), which may catalyze the degradation of RNA into shorter fragments [123,124] and play a relevant role in nutrient mobilization under Pi-starvation. It seems likely that sRNA production upon Pi-starvation is an effect of RNA degradation by different types of nucleases. Thus, degraded RNA may serve as a source of Pi necessary in emerging plant organs. We found also two genes encoding oxalate oxidases, the expression of which was signi cantly downregulated during Pistarvation. This class of genes is responsible for the inactivation of oxalic acid, which mediates fungalplant pathogenesis in barley [125].

Conclusion
To conclude, our studies provide comprehensive data sets, which may serve as a rich platform for the characterization of barley responses to Pi-starvation at an RNA level. Furthermore, our data may be used as a reference tool for parallel studies in other crop plants.

Plant material
Three biological replicates of barley root and shoot samples were analyzed. One replicate consisted of three plants growing in a single pot containing 1.5 kg of soil mixed with sand in a 7:2 ratio. Material was collected from the barley line Rolap (obtained from the Institute of Plant Genetics of the Polish Academy of Sciences, Poznań, Poland [126]) growing under low-Pi (8 mg P/kg soil) and Pi-su cient conditions (after addition of 60 mg P/kg soil), as described before [48]. On the 23 rd day after sowing plant shoots (Zadoks decimal code 22-23 [127]), they were cut off and fresh tissue weight was measured.
Immediately afterwards, shoots and roots were collected and frozen in liquid nitrogen to be kept at -80 °C until use.

Pi concentration measurements
Measurements of inorganic phosphate level were performed according to the protocol we have described before [46]. The samples were measured in two technical and three biological replicates using an In nite F200 Pro (TECAN, Switzerland).

RNA isolation
Four procedures of RNA isolation were used, depending on the following experiments: (i) small RNA expression level analysis (ddPCR using TaqMan TM MicroRNA assays, NGS of small RNAs, Northern hybridization); (ii) transcriptome analysis; (iii) degradome -PARE (Parallel Analysis of RNA Ends) [128] analysis for mRNA cleaved by miRNA; or (iv) validation of RNA-Seq data using ddPCR.
(i) Small RNA expression level analysis RNA isolation was performed using a modi ed method allowing enrichment of small RNAs, according to the detailed protocol we published before [64].
(ii) RNA for RNA-Seq RNA was extracted from a 100 mg sample using RNA extraction buffer [104] and a Direct-Zol RNA MiniPrep Kit (Zymo Research). According to the Lexogen's SENSE mRNA-Seq Library prep kit v2 user guide DNase treatment step was omitted to avoid RNA hydrolysis.
(iii) RNA for degradome analysis Procedure of RNA isolation from barley root and shoot (growing in low-Pi conditions) used for degradome pro ling was performed using a method described by German et al. using RNA extraction buffer [128], along with some modi cations that we have described previously [64,129].
(iv) mRNA-Seq data validation To validate the transcript level of signi cantly changed genes, we used precise dd-PCR analysis. To isolate RNA for these analyses, we used a Direct-Zol RNA MiniPrep Kit (Zymo Research) with some modi cations that we have described in detail previously [104]. The RNA material was treated using DNase I enzyme from the above kit (Zymo Research).

Preparation of NGS libraries
We  according to the manufacturer's protocol and as previously described [48]. A 420 pg of Spike-In RNA Variants SIRV-set3 (Lexogen) was added to 1500 ng of total RNA. ERCC mix was used for Spike-in analysis.

Data analysis
Differences in small RNAs, RNAs levels, and preliminary degradome data analysis were performed using a CLC Genomics Workbench (Qiagen Aarhus A/S).
Small RNA data analysis The trimming procedure was used with default settings for quality trimming (quality score limit 0.02), adapter trimming, and for removal of small RNAs longer than 25 nt and shorter than 18 nt. Reads were extracted, counted, and normalized per 1,000,000 reads. Then, we set up the Experiment analysis for samples derived from roots and shoots separately (two-group comparison, unpaired

Degradome data analysis
Degradome construction was performed using two different approaches, which allowed for a more indepth analysis (Fig. 2). In the rst method, the raw sequencing reads were processed by Cutadpt program (https://cutadapt.readthedocs.io/en/stable/) to trim low-quality and adapter sequences. Only sequences of length 15 nt and above were selected for further analyses. The processed sequencing reads were aligned to the reference sequences using bowtie. The count of 5'-end marked cleavage sites was scored by Perl script and normalized to the depth of sequencing and total signal for each of the reference transcripts. The putative miRNA:target pairs were predicted by a custom program (targetSeek) which included the following steps: (i) calculation of perfect match MFE (minimum free energy); (ii) RNAplexbased (Vienna package) screening for sRNA:transcript pairs; (iii) ltering number of bulges and length of sequence overhangs by MFE (percent of the perfect MFE match); and (iv) calculation of prediction score using a penalty schema for loops, bulges, and G:U wobble pairing. In the second approach, we used the PAREsnip2 software [130] to generate t-plots in conjunction with ve databases (Fig. 2). Potential miRNA targets are classi ed into one of ve categories, where category 0 indicates the best miRNA-target match.
The lower the alignment score, the better the alignment between the sRNA and the target site [130]. During PAREsnip2 analysis, we set the Fahlgren and Carrington targeting rules to permit a mismatch or G:U wobble at position 10 [131]. the CLC Genomics Workbench (QIAGEN) software, as previously described [48]. Among potential differentially expressed transcripts only those which went through restricted Bonferroni p-value correction were considered as differentially expressed genes (DEGs), which are listed in Additional le 8.

GO analysis
Gene ontology (GO) analyses were performed using the PANTHER in silico tool [132]. The overrepresentation binomial tests classi ed DEGs within GO domains (cellular component, biological process, and molecular function) with Bonferroni-corrected p-value < 0.05.

cis-regulatory motif localization within DEG promoters
To analyze the enrichment of Pi-related cis-regulatory motifs, we extracted 2000 bps upstream of transcription start site from each identi ed DEG. Such data were directly screened to look for any either P1BS-or P-responsive PHO element consensus sequences using multiple promoter analysis with the ddPCR To determine the absolute copy number of genes encoding IPS1, SPX-MFS1, endonuclease S1/P1, 3'-5' exonuclease, oxalate oxidase, and oxalate oxidase 2, we performed ddPCR using either EvaGreen Supermix (Bio-Rad) or TaqMan Assay (Bio-Rad) for mature miR827, according to the protocols previously described [48]. To normalize the copy number of miR827, we ran ddPCR for the ARF1 reference gene using the TaqMan Assay id: AIMSIL4 (Thermo Fisher Scienti c). Absolute gene expression was shown as normalized copy number per 1000 copies of the barley ARF1 reference gene. All speci c primers and probes (mature miR827, U6 snRNA) used in this paper are listed in Additional le 25.
Northern blot of mature miR827 To determine the mature miR827 expression level, we performed northern blot hybridization using a speci c probe for analysis. All steps of these experiments were done according to a detailed protocol as described previously [64]. 10 µg of each RNA sample was run alongside a radioactively labelled Decade Marker (Invitrogen, Thermo Scienti c) on a 15% polyacrylamide gel with 8M urea. The miR827and U6 probe sequences are available in Additional le 25. The Decade Marker (Ambion) was loaded to control the length of the tested RNAs. To calculate band intensity, we used the ImageQuant TL 8.1 software (GE Healthcare Life Sciences).

Consent for publication
Not applicable.

Availability of data and materials
The datasets generated and/or analyzed during the current study have been submitted to GEO database (accession numbers: GSE145423, GSE145425, GSE145426, GSE145427).

Competing interests
The authors declare that they have no competing interests.  identi ed in this study. Green color = DEMs from shoot, brown color = DEMs from root. Table 2. List of genes predicted in degradome analysis to be guided for cleavage by the most up-and down-regulated sRNAs identi ed in this study. Table 3. List of genes predicted in degradome analysis to be guided for cleavage by putative regulatory sRNAs (identi ed as DES) with best scoring matches.

Additional Files
Additional le 1. Normalized copy numbers of barley IPS1 gene transcript in low-Pi treated root material. DdPCR was performed to examine the absolute gene expression of the barley IPS1 gene. Obtained copy numbers were normalized per 1000 copies of the ARF1 reference gene transcript. Asterisks indicate a signi cant differences (*p-value < 0.05) calculated using two-tailed Student's t-tests.
Additional le 2. Characteristic of reads obtained from small RNA deep sequencing.