Investigation of the Dynamic Regulation Mechanisms of Small Non-Coding RNA in Response to Drought Stress in Qingke Barley


 Background Drought is a common abiotic stressor that exerts a great influence on grain security worldwide. The mechanisms underlying the small non-coding RNA regulation of plant drought responses remain unclear.Results In this study, drought-resistant (Xila-16, Xila) and drought-sensitive (DingqingHYG, Diqing) ecotypes were selected. A systematic analysis was performed on the microRNAs and tRNA-derived fragments in small RNA sequencing data. By predicting the target genes of differentially expressed miRNAs and analysing their pathway, this study identified HVUL2H01030.2, HVUL2H06276.2, HVUL2H00175.2 and other important target genes. The analysis also identified base excision repair and other adversity stress-related pathways. tRNA-derived fragments (tRFs), as a novel type of small noncoding RNA, exist widely in organisms, but no study has shown that tRFs play a role in the drought resistance of qingke barley. Based on systematic identification and characteristic analysis of tRFs in small RNA sequencing data, this study found that qingke barley had a tendency to produce more tRFs under drought stress. These tRFs were widely expressed, showed specific tRNA cleavage modes and had conservative intertreatment cleavage and other characteristics.Conclusions Our findings lay the foundation for further investigation of the action mechanism of such novel small noncoding RNAs in the drought resistance of qingke barley.


Investigation of the Dynamic Regulation
, suggesting that the expression of small RNAs is variety-specific.
Barley is a crop of major economic significance, so investigating the specific expression of microRNAs in drought-resistant barley varieties is of vital significance to further understanding the biological mechanism of drought resistance. A previous study has reported the identification of three microRNAs-miR1435, miR5024 and miR7714-that have specific expression in the roots of drought-resistant barley varieties [5]. The expression of microRNAs under drought stress is also found to be tissue-specific. Under drought stress, the specific expression of some differentially expressed microRNAs was observed on only the leaves or the roots of barley [6]. In a study of rapeseed, disease resistance protein (DIRP), drought-responsive family protein (DRRP), early responsive to dehydration stress protein (ERD) and stress-responsive alpha-beta barrel domain protein were identified as target genes of differentially expressed microRNAs, and these target genes play a role in regulating the responses of rapeseed to salt and drought stress [7].
Over the past few years, with the development and application of high-throughput sequencing, many other types of small noncoding RNA have been discovered and investigated. Among them, RNA-derived RNA fragments (tRFs) have been identified in numerous eukaryotes [8][9][10][11]. Several studies have shown that tRFs have many different functions. For instance, some studies have shown that some tRFs exist in small RNA pools where Argonaute and Piwi complexes interact with each other; additionally, some tRFs have been treated by Dicer, suggesting that some tRFs can function similarly to microRNA [12][13][14]. tRFs can be detected in many cells and tissues, but studies have also shown that the production of tRFs can be induced in various adverse circumstances, including drought, high temperature and high pH [9]. A previous study has also shown that tRFs may inhibit protein synthesis [15]. In addition, angiogenin-induced tRFs tend to inhibit protein translation [16].  [20] and RFAM [21] databases. The results predicted by cmsearch were annotated as the miRNA precursor, and further analysis and confirmation of mature microRNA sequences was necessary. miRdup [22] is a computational predictor used to identify the most likely miRNA position in the given pre-miRNA or to verify candidate miRNAs. Trained by experimentally verified miRNAs from miRBase [23], this classifier is designed with miRNA/miRNA* duplex. There are differences among miRNAs in the sequences and species, mainly involving bistability-related structural characteristics. Thus, when predicting mature microRNA sequences with miRdup, a microRNA sequence-trained model based on the species in the miRBase database was used.

Analysis of MicroRNAs in Small RNA Sequencing Data
The alignment analysis software Bowtie2 [24] (v2-2.3.5) was used to map the clean sRNA sequencing reads onto the reference genome; one mispairing was allowed during alignment. miRDeep2 [25] was employed to analyse the expression level and secondary structure of known microRNAs and predict novel microRNAs. DESeq2 [26] was used to perform a differential expression analysis of microRNAs (screening conditions of differential miRNAs: padj ≤ 0.05 and |log2(Fold_change)| ≥ 1. For plant samples, target genes were predicted using psrobot [27] software.

Differential Expression and Specific Analysis of tRFs
In the differential expression analysis of tRFs, the criteria of differential expression were defined as P-value < 0.01 (P-value was obtained using the Wilcoxon rank sum test and permutation test) and fold change between samples more than 2 or less than 0.5.
First, based on the expression levels of individual tRFs, the Jensen-Shannon divergence (i.e., JS score) under experimental conditions was calculated according to the information entropy method [29]; the maximum JS score of a tRF was adopted as its condition-specific score. A higher condition-specific score indicated that the gene concerned had a heterogeneous expression under various experimental conditions and a stronger expression bias. A lower condition-specific score suggested that the gene concerned had less obvious differences in expression.

MicroRNA Identification
A total of 3,259 microRNAs were predicted using previously published genomes. miRDeep2 was employed to analyse the expression level and secondary structure of microRNAs predicted. MicroRNAs having different mature miRNA sequences, precursor sequences and positions were regarded as different novel miRNAs. Based on sequencing data, a total of 844 novel microRNAs were predicted by miRDeep2.

Preference Analysis and Differential Expression Analysis of MiRNAs
The course of processing from the miRNA precursor to mature miRNAs requires Dicer cleavage, and because of the specificity of cleavage sites, mature miRNA sequences have some base preferences. The first base of microRNAs strongly prefer base T. By analysing the base preferences of microRNAs predicted on the basis of genomes, the accuracy of this prediction was further verified in this study. As shown in Fig. 2A, the preferences of microRNAs predicted in this project were relatively consistent and conformed to general preference laws.
miRNAs are tissue-specific and time-dependent, and their expression may vary across different tissues or different stages of the same tissue. These differentially expressed miRNA play a very important role in gene regulation. In this study, DESeq2 was used to conduct a differential analysis, and the screening conditions of differential miRNAs were defined as padj ≤ 0.05 and |log2(Fold_change)| ≥ 1.
In this study, there were relatively few differentially expressed miRNAs, and differentially expressed miRNAs fell mainly within DQ-S-4 versus XL-S-4. That is, there was an abundance of differential miRNAs between the drought-resistant variety and drought-

Prediction of the Target Genes of Differentially Expressed MiRNAs
For plant samples, the psrobot software was used to predict target genes and a total of 45 target genes of differentially expressed miRNAs were predicted. As indicated by an analysis on the expression of these target genes, differential expression was observed 4 h after drought stress treatment in four genes, three of them homologous to the 32 kDa protein in the NR library. Two of the target genes (HVUL2H01030. 2

KEGG Analysis on the Target Genes of Differentially Expressed MiRNAs
For each KEGG pathway, the hypergeometric test was used for the enrichment analysis.
Given that there were relatively few significantly different miRNAs, the p-value was relatively high in the enrichment results of target genes. In this study, we focused on the five highest-ranked pathways in terms of the p-value (from low to high). Many studies [31,32] have mentioned the important role of the most significant pathway, base excision repair, in stress response, especially under oxidative stress conditions. The protein coded by the gene HVUL5H44389.2 in this pathway is poly(ADP-ribose) polymerase, which influences energy stability, cell death and adversity tolerance in plants.  To study the expression abundance of different lengths of tRFs, a cluster analysis was performed on tRFs of different cleavage lengths ( Fig. 4C and D). The fragments produced by different tRNAs had unique high-expression fragment lengths, which further demonstrated that different tRNAs had specific modes of cleavage.

Inter-treatment Dysregulation of tRF Expression
Empirical cumulative distribution plots exhibit the empirical cumulative distribution of the expression levels of tRF-3 and tRF-5 in XL and DQ (Fig. 5). For tRF-3, the curve of the drought-resistant variety XL showed a more obvious right average than the curve of the drought-sensitive variety under drought treatment, suggesting that it was easier for XL to produce more tRF-3 under drought treatment. tRF-5 presented similar trends in the two varieties, which is to say that they both showed the right average of curve under drought treatment, suggesting that it was easier for both varieties to produce more tRFs under drought stress.
In the analysis of tRFs, the criteria of differential expression were defined as P-value < 0.01 (P-value was obtained using Wilcoxon rank sum test and permutation test) and fold change between samples of more than 2 or less than 0.5. Based on differential analysis of DQ and XL under drought and control treatments, DQ had no differentially expressed tRFs between the treatment and control, whereas XL had four differentially expressed tRFs between them. analysing related RNAseq data, it was found that a 32 kDa protein also showed differential expression as a target gene of differentially expressed microRNAs 4 h after drought treatment. Another study has also found a 32 kDa protein was as an important gene for dealing with stress [35]. Further studying this gene can lay a foundation for a deeper understanding of the mechanism of drought resistance in qingke barley in the future.

Discussion
tRNA fragments, as a novel type of small noncoding RNAs, are rarely explored in plants.
According to studies on animals and humans, the production of many tRFs is related to diseases and stress. For the purpose of revealing whether this type of small noncoding RNAs are produced and play a role in qingke barley, 18-30 nt tRFs at the 3′ and 5′ rRNAs were systematically identified in this study. It was found that the production of these tRFs was not random; instead, a very high degree of conservativeness in cleavage was observed, and they were produced in organisms by corresponding enzymes through the regular cleavage of tRNAs, so that tRFs must have played a role in biological regulation.
By analysing the cumulative distribution of the expression levels of tRFs in different samples, it was found that both varieties tended to produce more tRF-5 under drought stress treatment. In contrast, the drought-resistant variety produced more tRF-3 under drought stress, whereas the drought-susceptible variety showed no such pattern, suggesting that the two tRFs have different functions and production mechanisms, as well as that further studies should be conducted to explore their functions.

Conclusion
In this study, microRNAs and tRNA-derived fragments were analyed. There were a few differentially expressed miRNAs were identified. Some target genes of these miRNAs are stress-related. In these target genes, the gene HVUL5H44389.2 (poly(ADP-ribose) polymerase) involved in base excision repair pathway, which influences energy stability, cell death and adversity tolerance in plants. The tRFs were found had a tendency to produce more under drought stress in our study. These tRFs were widely expressed, showed specific tRNA cleavage modes and had conservative intertreatment cleavage and other characteristics. This study has offered helpful data resources and a basis for further study of the mechanisms related to drought resistance. Table   Table 1 KEGG pathway analysis on the target genes of DQ-S-4_XL-S-4 differentially expressed

Competing Interests
The authors declare no competing interests

Availability of data
The data was presented in the manuscript and the supporting materials. The raw reads data was submitted to the Short Read Archive (SRA

Supplementary Table Legends
Table S1 Result of microRNA differential expression analysis of DQ-S-1_XL-S-1.

Table S2
Result of microRNA differential expression analysis of DQ-S-4_XL-S-4.

Table S3
Result of microRNA differential expression analysis of DQ-S-24_XL-S-24.

Table S4
Result of microRNA differential expression analysis of DQ-S-48_XL-S-48.

Table S7
Result of microRNA prediction using miRdup.