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 growth and shoot branching (Fig. 1A). Shoot fresh weight of plants at 23rd day post-sowing (dps) was significantly reduced, in comparison with control plants, with average mass 8.8 g for stressed plants and 18.5 g for plants growing under Pi-sufficient conditions (P = 0.001) (Fig. 1B). We observed a significantly decreased concentration of Pi ions, with only 0.48 µmol Pi per g of fresh root weight and 4.2 µmol Pi per g of fresh shoot weight, 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-deficient 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 file 1). Taking validated plant material, we performed tripartite deep-sequencing analysis to: (i) identify Pi-responsive small RNAs of significantly changed expression level using sRNAs Next-Generation Sequencing (NGS), (ii) elucidate changes in the barley transcriptome upon low-Pi using RNA-Seq, and (iii) seek mRNA targets for Pi-responsive sRNAs through degradome sequencing (Fig. 2).
Barley plants express an organ-specific set of microRNAs in response to low-Pi conditions
In the first step, we performed small RNA deep-sequencing analysis to elucidate specific miRNA and other sRNA molecules dysregulated by low-Pi treatment in barley shoots and roots. The average of 30.4 mln for roots and 25.2 mln raw reads for shoots were generated from 50 nt single-read sequencing using an Illumina System (Additional file 2). After adapter and quality trimming, we mapped reads to the miRBase Sequence Database (release 22) to annotate miRNA-derived sequences [59]. A set of parameters were used to define 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) significance of fold change (p-value < 0.05) was additionally verified using a restricted Bonferroni p-value test (Fig. 2).
We found 138 and 162 differentially expressed miRNAs (DEMs) annotated to the miRBase (p < 0.05) in barley roots and shoots, respectively. Only 25 DEMs were expressed in both examined barley organs (Fig. 3A, left panel). However, restricted Bonferroni p-value correction narrowed down set of miRNAs to 15 in shoots and 13 in roots (Fig. 3A, right panel). Four out of 25 common miRNAs passed the Bonferroni p-value test in both organs: miR399b, miR399a, and two miR827 isomiRs (Fig. 3A, right panel). We focused further only on those that passed the Bonferroni test. In both organs, most of the DEMs were significantly upregulated. Interestingly, out of 15 miRNA, only miRNA166d was downregulated in shoot under low-Pi (log2(fold change) = -1.18). In our previous work, we showed that miRNA166 is expressed in barley during different developmental stages [60]. 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 [61]. Similarly, only miRNA319b out of 13 DEMs was downregulated in low-Pi treated roots (-1.28). In a previous study, we presented data that Arabidopsis miR319 is a multi-stress responsiveness miRNA [62]. For instance, MIR319b gene expression level was downregulated in response to drought, heat, and salinity, but upregulated in response to copper and sulfur deficiency stresses [62].
miRNA399 and mRNA827 are highly induced by low-Pi stress in barley
In plants, a limited number of miRNA have been shown to be specifically and strongly induced by Pi limitation, including miRNA399 [63], miRNA778 [54], miRNA827 [50], and miRNA2111 [50, 64]. We found that a specific set of miRNAs was expressed in barley shoot or root under low-Pi (Additional file 3) and, in both organs, the most abundant DEMs represented miRNA399 and miR827. In shoot, while only these two miRNA families were induced (after Bonferroni p-value correction), we observed a more diverse response in root (Fig. 3A, right panel). Apart from miRNA399/miRNA827 induction, we found the following additional miRNA to be upregulated: miRNA9779, miRNA5083, miRNA156, miRNA5083, miRNA1511, and miRNA5072. Of these six miRNAs, only miR156 has been reported before as Pi-responsive in Arabidopsis studies [37, 50]. The other five root DEMs remain uncharacterized in all available databases. The miR156 isomiRs were also found dysregulated in shoot, but none of them pass the Bonferroni test. The miR156 has been shown to guide the cleavage of SQUAMOSA-PROMOTER BINDING-LIKE (SPL) transcription factors in Arabidopsis [37] and rice [65, 66]. SPL genes are a key regulators of developmental transitions and control the growth and number of tillers. Thus, miR156 may integrate environmental stress signals with plant development [67] and trigger the reduced branching phenotype under growth conditions of low-Pi (Fig. 1A). Altogether, our results suggest that there is a more complex response to low-Pi deficiency 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. NGS data were then validated by complex analysis of mature miR827 in all samples taken for deep sequencing. The absolute expression level of miR827 is significantly upregulated in both shoots and roots under a low-Pi regime (Fig. 3B). These data were confirmed by northern blot hybridization (Fig. 3C) and degradome profiling showed that SPX-MFS1 transcripts have been recognized and cleaved by Ago protein associated with miR827 in barley (p = 0.014) (Fig. 3D). The cleavage site may differ in various SPX-MFS1 isoforms (Additional file 13, Additional file 18). In Arabidopsis, miR827 has been shown in multiple studies to target the 5’-UTR of the NITROGEN LIMITATION ADAPTATION (NLA) gene [68, 69]. In rice, the OsNLA mRNA has a ‘degenerate’ osa-miR827 potential cleavage site, but does not cleave the OsNLA transcript in vivo [51, 52, 70]. Likewise, we did not find NLA mRNAs to be targeted by any of the identified hvu-miR827 isomiRs in our degradome records. The NLA gene encodes a ubiquitin E3 ligase with RING and SPX domains, which interacts with the PHO2 to prevent the excessive accumulation of Pi [71].
Our results are consistent with sRNA sequencing data published for Arabidopsis [37, 50] and Nicotiana benthamiana L. [72]. 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. 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 [73]. 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 deficiency [74]. The level of PHT1 proteins is negatively controlled by the PHO2, which is suppressed by miR399 in a feedback loop (see model in Fig. 7) [75]. A high level of miR399 molecules was detected in Arabidopsis wild type rootstocks grafted with miR399-overexpressing scions [37, 76]. 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 [76].
Different classes of small RNAs in barley accumulate in an organ-specific manner under low-Pi regime
In the next step, the remaining reads of other sRNAs were annotated to particular classes of barley cDNAs derived from the Ensembl Plants database (release 40), both separately and to all of them (Fig. 2). All annotated sequences are listed in Additional file 3. We found that differentially expressed other small RNAs (DESs) affected by low-Pi conditions were represented by 199 unique sequences (0.01% from total unique reads obtained under low-Pi) in barley shoots (Additional file 4) and by 1796 (0.13%) unique reads in roots (Additional file 5). Among the unannotated reads of sRNAs in roots, the highest fold change was observed for a 19 nt 5’-ACCTACTCGACCTCGGCCG-3’ molecule (log2(fold change) = 8.02, induction) and a 22 nt 5’-CTAATACCGGATACGCGAACCG-3’ molecule (-5.87, repression). The first (19 nt) molecule showed a perfect match to either the intergenic region of barley chromosome no. 5, soil bacteria (mesorhizobium), or Linum usitatissimum L., while the second molecule (22 nt) encodes 16S rRNA. Furthermore, in roots, the most abundant small RNA was a 25 nt 5’-ACCGACCTACTCGACCCTTCGGCCG-3’ molecule (15847.7 and 65590.5 mean of normalized counts in barley root in control and low-Pi conditions, log2(fc) = 2.82). This small RNA matched several barley loci encoding SSU (small subunit) rRNAs (Additional file 5). In total, 1070 reads were successfully mapped to cDNA sequences annotated in the Ensembl Plants database, while 726 unique sequences remained without match.
In shoot, we found 199 DESs under the low-Pi regime. Among the other sRNA reads, the highest fold change was represented by a 24 nt 5’-AAGATTGGTTGGTTGGTTGGGTCT-3’ molecule (log2(fc) = 8.72, induction). This 24 nt molecule is a part of transcript encoding a putative pentatricopeptide repeat (PPR) protein. The PPR protein family facilitates the processing, splicing, editing, stability, and translation of RNAs in plants [77]. The most abundant small RNA was a 19 nt 5’-GGGCCTGTAGCTCAGAGGA-3’ molecule (9471.5 and 49914.1 normalized mean counts in barley shoot in control and low-Pi, respectively, log2(fc) = 2,45). This sRNA was mapped to the barley genomic loci (EPlHVUG00000039813), which encodes arginyl-tRNA (trnR-ACG) and a cDNA encoding uncharacterized protein (HORVU2Hr1G084630) which is likely involved in carbon fixation. Altogether, 116 out of 199 differentially expressed small RNAs (DESs) were annotated to the barley Ensembl Plants database, where 83 sequences remained without match (Additional file 4). Interestingly, the pool of DESs was selective, considering organ-specific expression change, providing only three unique sequences that were significantly changed in both barley organs under low-Pi regime (Fig. 4A, left panel). These molecules were: (i) 20 nt 5’-AGTAGAGGTCGCGAGAGAGC-3’ (log2(fc) = 2.01 in root and 1.16 in shoot, respectively) annotated to the 26S rRNAs, (ii) 24 nt 5’-ATTCTCCGCGTCGGATACCTGAGA-3’ (3.69 in root and 2.07 in shoot) encoding the barley MYB21 transcription factor, and (iii) 21 nt 5’-TGCCAAAGGAGAACTGCCCTG-3’ (4.64 in root and 6.27 in shoot) mapped to the intergenic region of barley chromosome no. 3 (Additional file 4). In addition, total numbers of 166 DESs (83%) in shoots and 1560 DESs (87%) in roots were significantly upregulated after exposure to low-Pi stress (Additional file 4, Additional file 5, Fig. 7).
Different lengths 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 fluctuated 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 (9.55%), 23 nt (8.54%), and 25 nt (3.52%) sequences (Fig. 4B, Additional file 6). The DESs obtained from low-Pi roots were mostly annotated to protein-coding mRNAs (38.54%), rRNAs (34.17%), and non-translating RNAs (19.49%). Below 5% of overall DESs, we found a number of remaining cDNA classes, such as snoRNAs (2.49%), tRNAs (2.47%), SRP-RNAs (1.17%), snRNAs (0.95%), and pseudogenes (0.65%). In the case of shoot samples, 85% of annotated DESs represented only protein-coding mRNAs (47.87%) and non-translating RNAs (36.49%) (Fig. 4A, right panel; Additional file 7). We did not find any DESs annotated to the snRNAs, SRP-RNAs, or tRNAs from barley shoot upon low-Pi.
Again, barley roots exhibited a more diverse pool of Pi-responsive small RNAs, which may trigger developmental adaptation of the root to Pi scarcity. Additionally, rRNA-derived small RNAs are highly accumulated in barley roots. We believe that such sRNA may be further processed, serving as a Pi source to compensate Pi deficiency.
Identification of barley genes responsive to Pi starvation
As we observed that most of the other sRNAs in shoot were derived from either protein-coding mRNAs or non-translating RNAs, we tested whether this observation correlated with gene expression changes of polyadenylated RNAs in barley shoot under low-Pi. Experiments were performed in three biological replicates of plants grown under low-Pi and control conditions. Paired-end sequencing reactions of the 150 nt reads were performed using an Illumina System. Total read numbers from six samples were mapped to the barley reference genome from Ensembl Plants Genes 42 (Hordeum vulgare IBSC v2). The library’s quality and sequencing accuracy were verified carefully (i) by adding Spike-in RNA Variant Control Mixes (Lexogen) (Additional file 8) and (ii) by quality trimming. Differentially expressed genes (DEGs) were selected based on fold change calculations (p < 0.05) and the selective Bonferroni test. Among 98 of identified DEGs, the transcripts of 56 annotated loci were significantly upregulated, while those derived from 42 loci were downregulated in Pi-starved barley shoots (Additional file 9). 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).
The highest enrichment of DEGs was found in the GO terms, either (i) belonging to the cellular components of the chloroplasts; (ii) showing catalytic activity, either ion or metal binding properties; and (iii) involved in the various biological and metabolic processes related to stress responses and plant defense (Fig. 5A). A major set of upregulated DEGs represented genes involved in the Pi signaling. Among them, we found genes encoding, such as IPS1 (log2(fc) = 5.89) [49], inorganic pyrophosphatase (PPase, 4.01) [78], SPX-domain containing protein 5 (SPX5, 3.44) [79], phosphate transporter PHOSPHATE 1–3 (PHO1-3, 2.97) [80], SPX-MFS2 (2.79) [51], haloacid dehalogenase-like hydrolase (HAD1, 1.95), [81] and five different purple acid phosphatases (PAPs) (Fig. 5C; Additional file 9) [82]. Interestingly, four genes were induced to a higher extent than the low-Pi stress marker, IPS1 gene. These genes encode ferredoxin (FD1, log2(fc) = 14.20), mitochondrial-processing peptidase (13.35), chlorophyll a/b binding protein (8.90), and alpha-amylase (7.30), and are engaged in photosynthesis, redox reactions, reactive oxygen species (ROS) homeostasis, and co-ordinated mobilization of nutrients. Chloroplasts and mitochondria are the organelles with the highest Pi requirements. Strong FD1 gene upregulation most likely reflects the accumulation of reduced ferredoxin in chloroplasts. Pi limitation lowers the capacity to process incoming light and enhances starch accumulation in chloroplasts, thereby leading to photoinhibition [83, 84]. Within the category of genes that were significantly downregulated, most of them were related to stress and defense responses (Additional file 9); for instance, uncharacterized protein (HORVU2Hr1G030090, -6.50), oxalate oxidase (-4.41) [85], beta-sesquiphellandrene synthase (-3.41), glutamate carboxypeptidase (-3.17), chalcone synthase (-3.05) [86], or caleosin-like protein (-2.95) (Fig. 5C). Only two repressed genes are known to be directly involved in Pi signaling and metabolism, SPX-MFS1 (-2.58), targeted by miR827 [52] and probable inactive purple acid phosphatase (-1.75). Additionally, two genes encoding laccases (LAC19-like, Additional file 9), cell wall-localized multi-copper oxidases, were significantly downregulated (-2.10 and − 2.44) in our RNA-Seq data. Laccases are involved in copper homeostasis (Fig. 7) and lignin biosynthesis, and have been shown to be targeted by miR397 in maize [87] and Arabidopsis [88]. Recently, it was shown that miRNA397b-LACCASE2 regulates root lignification under either water or Pi deficiency in Arabidopsis [88]. MIR397 are not P-responsive genes and, thus, we did not find significant expression level changes in any of the barley miR397 isomiRs in both examined low-Pi stressed barley organs.
Furthermore, key genes encoding proteins involved in the nitrate and phosphate cross-talk were affected by low-Pi conditions in barley shoots, such as NIGT1 (NITRATE-INDUCIBLE, GARP-TYPE TRANSCRIPTIONAL REPRESSOR 1) transcription factor (3.80) [89, 90] and nitrite reductase (1.98), as well as high-affinity nitrate transporter NRT2.1 (NITRATE TRANSPORTER 2.1) (-2.60) [91]. Absolute quantification of a few select 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 analysis (Fig. 6A). We confirmed statistically significant changes (p < 0.05) in normalized copy number (per 1000 copies of the ARF1 reference gene) of all genes taken for analysis.
In 2018, Ren et al. published RNA-Seq data describing the barley transcriptome under low-Pi stress [92]. The authors compared the transcriptomes of two barley genotypes with contrasting low-Pi stress tolerance. In roots, they observed 28 DEGs classified into the following functional groups: phosphate transport, transcription, lipid metabolism, metabolism, and phosphorylation/dephosphorylation [92]. Likewise, our RNA-Seq data from barley shoot discovered the DEGs involved in all mentioned functional groups. Of the 28 DEGs described before, we found the same four DEGs in our shoot transcriptome analysis: (i) GLYCEROPHOSPHODIESTER PHOSPHODIESTERASE 1 (GDPD1) gene (HORVU3Hr1G079900, log2(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 [92]; however, they appeared from different barley genome loci than the five PAPs we found in our study. Other related RNA-Seq data published for either wheat (Triticum aestivum L.) [93], rice [94], soybean [95], Plantago major L. [96], and maize [97] demonstrated similar molecular patterns to those in our study. Interestingly, an elevated abundance of sRNAs has been associated with the upregulation of two types of nucleases (endonuclease S1/P1 and 3’-5’ exonuclease), which may catalyze the degradation of RNA into smaller components [98, 99] and play a relevant role in nutrient mobilization under Pi scarcity. Taken together, the 98 DEGs identified in our RNA-Seq analysis have broadened the landscape of low-Pi responsive genes modulating stress tolerance in barley (Fig. 7, Additional file 9).
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 [100]. Previously, we have shown the importance of the P1BS motif (PHR1 binding sequence, consensus GnATATnC, [101]) and P-responsive PHO elements (consensus ATGCCAT, [102]) in the expression efficiency of the barley PHO2 gene [43]. 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 [103]. Likewise, we hypothesized that regulatory regions of the identified DEGs had Pi-responsive motifs, which may be bound by PHR TFs, causing gene expression dysregulation. To confirm this hypothesis, we extracted DNA sequences from the 2000 bp region upstream of the predicted transcription start sites from all 98 DEGs (Additional file 10). 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 [104]. We confirmed 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 file 11) and 17 DEGs with at least one P-responsive PHO element (Fig. 6B, Additional file 12). The most over-represented motifs were found in the promoters of genes encoding sulfoquinovosyl transferase SQD2-like [105], phosphoenolpyruvate carboxylase 1-like [106], and pyridoxal phosphate-dependent transferase [107]. Each of the genes harbor three P1BSs and one P-responsive PHO element, as well (Additional file 9).
The presence of crucial Pi-responsive cis-regulatory elements within the promoter regions of more than 50% of identified 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 [101] and PHR2 in rice [108], 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 [103]. 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 [43]. The PHO element can be bound by PHR-like transcription factors in barley plants, as well [43], and has been found in the promoters of many DEGs in independent Arabidopsis [109, 110] or soybean [95] studies.
Degradome profiling describes post-transcriptional regulatory network of identified DEMs
After identification 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 small RNAs directly involved in RNA degradation. The other sRNAs 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 files 13–24).
First, we analyzed the potential target mRNAs for a broader set of differentially expressed miRNAs with significant fold change (p-value without Bonferroni correction), taking 162 DEMs from shoot and 138 DEMs from root, respectively (Additional file 3). A total of 1964 scores were obtained for shoot DEMs (1303 using the TargetSeek approach and 661 using PAREsnip2) (Additional file 13, Additional file 17), while in root there were 1515 records (921 and 594, respectively) (Additional file 15, Additional 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 obtainable score (0). A majority of records corresponded to different miR399 and/or miR827 isomiR and their known targets PHO2 or SPX-MFS1/SPX-MFS2, respectively. Those post-transcriptional gene regulations are crucial for proper plant responsiveness to Pi scarcity [111]. In shoot, the best scoring miRNA:mRNA match was found for mature miRNA399d (21 nt), which guides cleavage within the 5’-UTR of PHO2 mRNA in position 857 (TargetSeek: score = 0, MFE = -41,1, PAREsnip2: score = 0, MFE = -41.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 fungus and anthocyanin biosynthesis [112, 113]. 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 file 17). The SPL genes are post-transcriptionally regulated by miR156 and control shoot branching in rice [65]. Crosstalk between N and P has been extensively studied in the last few years [114]. The N-responsive miR171d (21 nt) was found in our degradome data to guide cleavage of mRNA encoding scarecrow-like protein 6 (SCL6, HORVU6Hr1G063650, TargetSeek: score = 2, MFE = -37.1, PAREsnip2: score = 1, MFE = -37.1) in a Pi-dependent manner. The SCL6 TFs may be involved in the regulation of primary root elongation in response to N shortage through PTGS mediated by miR171 action in Arabidopsis [115, 116]. Another N-related example was found for the ARF8 gene (HORVU6Hr1G058890), which is involved in root initiation and emergence [117]. Both algorithms proved its cleavage was guided by the 21 nt mature miR160 (TargetSeek: score = 1, MFE = -45.3, PAREsnip2: score = 1.5, MFE = -43.8).
In contrast, in roots, where the pool of identified DEMs is more diverse, the PTGS regulator network involves many TFs with various functional domains and/or different ARFs (Additional file 15). For instance, the 21 nt mature miR444b targets mRNA encoding MADS-box transcription factor 57 (HORVU6Hr1G073040, TargetSeek: score = 0, MFE = -40.1, PAREsnip2: score = 0, MFE = -40.1), which controls the outgrowth of axillary buds in rice [118]. The well-known abiotic stress-responsive miR319c (21 nt) [6] was found to target various isoforms of GAMyb TF (HORVU3Hr1G079490, TargetSeek: score = 5, MFE = -34.6, PAREsnip2: score = 3, MFE = -35.4). The same predicted miR319 target has been previously found in sRNA deep sequencing of Pi-starved tobacco plants [72]. Its target, the plant GAMyb TFs, has been shown to activate gibberellin-responsive gene expression of α-amylase in barley [119, 120]. Our degradome data indicate the possible cleavage of various mRNA isoforms encoding nuclear transcription factor Y subunit A-4 (NFYA4, HORVU4Hr1G005670, TargetSeek: score = 3, MFE = -27.6, PAREsnip2: score = 2, MFE = -28.1) mediated by miR169d (20 nt). This result is consistent with previous studies in Cunninghamia lanceolata L. showing NFYA4 mRNA to be targeted by miR169d (authors used PAREsnip approach) [121]. The NFYA4 protein is involved in many abiotic stress responses and may regulate the timing of transition from vegetative to reproductive phase [122]. All predicted DEM targets for both organs are listed in Additional files 13, 15, 17, and 21.
Putative miRNA-like molecules identified in degradome data
Degradome profiling 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 profiles matching root DESs (Additional file 16, Additional file 23) and 160 records (87 and 73, respectively) matching shoot DESs (Additional file 14, Additional file 19). Taking only either the most upregulated or the most downregulated sRNAs for degradome screening, we found five promising targets in roots and six in shoots. For example, in roots, the highly upregulated 20 nt sequence 5’-ACCGACCTACTTGACCCTTC-3’ (log2(fc) = 6.46) 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) (Fig. 4C). RNA-Seq data for potato (Solanum tuberosum L.) proved that expression of the MYB44 gene is highly downregulated under low-Pi in roots [123], which may be the result of our identified putative miRNA’s guided cleavage. 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 [123]. 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) (Fig. 4C). Additionally, the potential cleavage activity guided by sRNAs was found in the best scoring records (based on the match and MFE scores); that is, the 22 nt molecule 5’- CGCATTTATTAGATAAAAGGCT-3’ (3.19) has two predicted targets. Such 22 nt sRNA may bind and guide cleavage of mRNAs encoding two uncharacterized proteins; one potential senescence-associated protein (HORVU0Hr1G023930.1, TargetSeek: score = 0, MFE = -31.3) and one potential E3 ubiquitin-protein ligase (HORVU5Hr1G015600, TargetSeek: score = 0, MFE = -31.3). Meanwhile, the downregulated 21 nt sequence 5’-TCGGATCGCGGCGACGGGGGC-3’ (-1.02) binds to the mRNA’s coding region of RING-type E3 ubiquitin transferase containing the U-box domain (HORVU4Hr1G066070.1, TargetSeek: score = 6.5, MFE = -39.3), which functions as an activator of brassinosteroid signaling [124, 125]. Brassinosteroids have been shown to modulate root adaptation to adverse Pi conditions in plants [126, 127].
Analogous degradome screening was done for shoot data. Among all identified DESs, we found that the most upregulated sequence, 24 nt 5’- AAGATTGGTTGGTTGGTTGGGTCT-3’ (long2(fc) = 8.72), targets the 3’-UTR of mRNA encoding multiple organellar RNA editing factor 9 (MORF9, HORVU7Hr1G073170, TargetSeek: score = 16.5, MFE = -28.1). MORF9 proteins are required for RNA editing in plastid mRNAs, which may contribute to stress adaptation in plants [128, 129]. In both approaches, we found that the 21 nt sequence 5'-TGCCAAAGGAGAACTGCCCTG-3’ (6.27) targeted the same isoform of PHO2 mRNA (HORVU1Hr1G085570.3, TargetSeek: score = 6, MFE = -29.9, PAREsnip2: score = 3.5, MFE = -33.4). When we browsed the miRBase using this 21 nt sRNA as a query, we found high similarity to the osa-miR99a, exhibiting only one mismatch. Thus, we suspect that such sRNA may function as another miR399 isomiR in barley. Most dysregulated DESs were also found to target mRNAs encoding methyltransferase type 11 domain-containing protein (MT11), AAA-ATPase (At3g50940-like), lysine-specific demethylase 5A (LSD), or an uncharacterized protein with a predicted transmembrane domain (HORVU3Hr1G036970). The best scoring records were found for the 21 nt sequence 5’- TTAGATGACCATCAGCAAACT-3’ (3.11), which targets mRNAs encoding SPX-MFS1 (PAREsnip2: score = 2.5, MFE = -27.5) and SPX-MFS2 (TargetSeek: score = 3.5, MFE = -26). This result suggests that such sRNA may exist as another miR827 isomiR. Moreover, this is consistent with the screening made for differentially expressed miRNAs, where miR827 targeted both SPX-MFS proteins, depending on the approach we used (Fig. 3A). In addition, the 20 nt sequence 5’- AATCGTCTTTACATCGGATG-3’ (3.24) was found to target various mRNA isoforms encoding glutaminyl-peptide cyclotransferase (HORVU7Hr1G003920, TargetSeek: score = 4.5, MFE = -20.7) (Fig. 4C), which may be involved in plant defense reactions [130].
Some other interesting Pi-related targets 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-affinity nitrate transporter-activating protein (HORVU5Hr1G115500), MYB-like TF (HORVU7Hr1G027370), and stress-induced TF NAC1 (HORVU5Hr1G111590) were found. Interestingly, among the 98 DEGs identified 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 identified IDs classified for differentially expressed small RNAs. These data suggest that the generated small RNAs may serve as an internal Pi source through nucleolytic degradation, rather than contributing to PTGS. All predicted DESs targets are listed in Additional files 14, 16, 19, and 23.