Identication and Proling of microRNAs and Their Target Genes in Tibetan Hulless Barley in Response to Barley Leaf Stripe (BLS) Infection

Tibetan hulless on where Here, we compared the miRNA proles before and after BLS in two Tibetan hulless barley genotypes: Z1141, a BLS-sensitive wild variety, and Kunlun14, a BLS-tolerant hybrid variety. A total of 36 conserved and 56 novel miRNAs were identied. Of these, 10 conserved and 10 novel miRNAs exhibited signicantly changed expression between the normal and infected leaves of Kunlun14, respectively, while 3 conserved and 5 novel miRNAs exhibited signicantly changed expression between the normal and infected leaves of Z1141, respectively. A total of 24 miRNAs were found in Z1141 and Kunlun14, and a further 546 putative target genes were predicted. Transcriptome sequencing analysis showed that among the 546 candidate genes, 131 had signicant differences in expression between the normal and infected leaves of Kunlun14 and Z1141. Gene ontology, pathway, and Blast analyses indicated 10 candidate target genes that were involved in the barley disease resistance. These 10 candidate target genes may be regulated by 7 miRNAs. According to quantitative real-time PCR results, the 10 targets were negatively correlated with their corresponding miRNAs after infection with BLS. study may improve our understanding of miRNA regulation in the response of Tibetan hulless barley to BLS infection in crop elds.


Background
Tibetan hulless barley (Hordeum vulgare L. var. nudum Hook. f.) is a self-reproducing annual species that produces naked grains. This plant is widely cultivated on the Qinghai-Tibet Plateau and has served as a staple food for the Tibetan people since the 5th century CE [1,2]. Grains of Tibetan hulless barley not only contain abundant nutritional components and physiologically-active compounds, but also are rich in edible pigments necessary the food manufacturing industry [3,4].
Barley leaf stripe (BLS) is a serious fungal disease affecting the yield and quality of Tibetan hulless barley in the Qinghai-Tibet plateau [5]. BLS is generally carried by grains and infects the leaves, sheaths, and spikes of host plants [6]. BLS is common in the main cultivation areas of Tibetan hulless barley in China, including Tibet, Qinghai, Gansu, Sichuan, and Yunnan [7]. Plant resistance of BLS is considered a complex trait, involving many disease resistance genes [8,9,10]. Quantitative trait locus (QTL) analysis is an e cient way to analyze some key traits when breeding for barley disease resistance [11,12]. For example, applying QTL analysis to progeny crossed by BLS-resistant and -susceptible varieties identi ed two disease resistant genes, Rdg1a and Rdg2a, which are located on the 2HL and 7HS chromosomes, respectively [8,9,13,14]. Overexpression of these two genes markedly enhanced the resistance of barley to BLS [15,16,17]. However, only a few genes have been discovered via QTL analysis. This suggests that we may have reached a bottleneck in the study of BLS through traditional gene mapping analysis, and thus new methods and perspectives are needed to explore novel plant genes directly or indirectly related to BLS resistance and resistance mechanisms.
MicroRNA (miRNA) is a type of non-coding small RNA consisting of 21-24 nt that plays various critical roles in the growth and development of plants, including how plants respond to abiotic and biotic stressors, such as pests and diseases [18,19,20]. High-throughput sequencing of microRNA has become a fast, e cient, informative, and cost-effective strategy for gene discovery [21,22,23]. Some miRNAs and their target genes contribute critically to regulating plant disease. For example, miR393 responds to bacterial disease by negatively regulating the auxin receptor genes TIR1, AFB2, and AFB4 [24]. Some arti cial miRNAs (amiRNAs) in barley participate in regulating yellow dwarf disease under low temperature stress conditions [25], while in tomato plants, miR403 can regulate leaf-shortening disease by targeting the AGO2 gene [26]. Researchers analyzed miRNA the changes in cotton roots infected by verticillium wilt and found that the expression levels of at least 65 miRNAs changed, most of which were non-conservative miRNAs [27]. It was discovered that viral-inducible OsAGO18 sequesters miR168 to ease the repression of rice OsAGO1 by miR168 to enable antiviral defense in infected rice [28] .
With the development of high-throughput sequencing, increasing numbers of plant miRNAs are being discovered. We are currently aware of 47,618 mature miRNAs covering 271 species, of which 10,414 mature miRNAs from 82 plant species have been identi ed in the miRNA database (http://www.mirbase.org/V21.0). Among the latter, 812 miRNAs originate from Arabidopsis, 325 from maize, 738 from rice, 241 from sorghum, and 126 from wheat. By contrast, only 71 miRNAs have been found in barley, and not a single miRNA has yet been reported from Tibetan hulless barley. Here, to ll this knowledge gap and explore disease-resistant miRNAs, high-throughput sequencing was used to construct a miRNA database of disease resistance and susceptible Tibetan hulless barley varieties. From our results, eight miRNAs and 10 candidate target genes related to BLS resistance were found and their expression patterns were analyzed in detail. Our study results not only provide a data foundation for the functional veri cation of BLS resistance genes, but also provide information for the future selection and breeding of new disease-resistant Tibetan hulless barley varieties.

Materials And Methods
Plant materials and BLS treatment.
Wild Tibetan hulless barley accession 'Z1141' and the hybrid variety 'Kunlun14' were used to investigate BLS resistance. The proportion of eld area (per square meter) of Kunlun14 and Z1141 was calculated from 2015 to 2017. Kunlun14 is a cultivar that can tolerate BLS, while Z1141 is susceptible to BLS ( Fig. 1). In this study, we used ZS, ZN, KL14S, and KL14N, as experimental materials. First, the infected leaves of Z1141 were collected to construct the infection system of BLS in the laboratory. To complete this process, ve pieces (5 mm × 5 mm) of one typical infected leaf from Z1141 were isolated, rinsed with 75% alcohol for 30 s and 2% sodium hypochlorite for 30 s, and then washed with distilled water three times. Finally, they were cultivated in a sealed Petri dish lled with potato dextrose agar (PDA) culture medium (200 g of dissolved peeled potato, 17 g of glucose, 16 g of agar to 1000 mL with distilled water) after drying off. The dish was placed in a biochemical incubator at 22-25 °C for 3 d. Non-contaminated mycelia were transferred to a new PDA culture dish four times until the colony was puri ed. First, three ensuing fungal cakes were removed with a perforator and cultivated in another new PDA culture medium, which was then placed in a biochemical incubator at 25 °C for 5-7 d until the whole Petri dish was lled with the BLS colony. Second, 300 grains of each Tibetan hulless barley were disinfected with 75% alcohol for 1 min, and then 150 grains were equally placed into three Petri dishes covered with the BLS colony, and another 150 grains were equally placed in three Petri dishes covered with the PDA culture medium without BLS colony. The above six dishes were then placed into a low-temperature incubator for two weeks at 6 °C. Finally, 300 grains of each variety were respectively planted in six owerpots with vegetative soil and cultured in an arti cial climate incubator (daytime: 14 h, at 20 °C; nighttime: 10 h, at 10 °C) for approximately 10 d. The second infected or normal leaf from the bottom of the seedling was removed, and leaves with at least three biological repeats were refrigerated at − 80 °C after freezing in liquid nitrogen.
Small RNA library construction and sequencing.
The total RNA of each leaf was extracted using TRIzol reagent (Invitrogen, CA, USA) according to the manufacturer's instructions, and the RNA purity was assessed using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). The RNA concentrations were then measured using the Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, MD, USA). The RNA integrity of each sample was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The RNA integrity was further checked through a 1% agarose gel analysis. All RNA samples were then sent to Novogene (Beijing, China) to construct libraries for small RNA sequencing, where single-end reads (SE50) were sequenced on an Illumina HiSeq 2500 platform (Illumina Inc., CA, USA).
A total amount of 3 µg total RNA per sample was used as the input material for the small RNA library. Sequencing libraries were generated using the NEBNext®Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA), following the manufacturer's recommendations, with index codes added to label each sample with its corresponding sequences. PCR ampli cations were performed using the Long Amp Taq 2X Master Mix, SR Primer for Illumina, and index (X) primer. The ensuing PCR products were then puri ed by 8% polyacrylamide gel analysis. DNA fragments corresponding to 140-160 bp-the length of small non-coding RNA plus the 3' and 5' adaptors-were recovered and dissolved in 8 µL elution buffer. Lastly, library quality was assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. The clustering of index-coded samples was performed on a cBot Cluster Generation System that used the TruSeq SR Cluster Kit v3-cBot-HS (Illumina), following the manufacturer's instructions. After completing the cluster generation, the library preparations were sequenced on an Illumina HiSeq 2500/2000 platform, and 50-bp single-end reads were generated.
Prediction of miRNA.
The small RNA tags were mapped to reference sequences by Bowtie [29], with mismatching not allowed, in order to analyze their expression and distribution vis-à-vis the barley reference genome (ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/hordeum_vulgare/cds/Hordeum_vulgare.Hv_IBSC_ PGSB_v2.cds.all.fa.gz) [30]. These reference-mapped small RNA tags were then used to search for known miRNA. The miRBase 22 served as the reference, with modi ed mirdeep2 [31] and srna-tools-cli software [32] used to obtain the potential miRNA and to draw the secondary structures. Custom scripts were used to obtain the miRNA counts as well as the base bias on the rst position of an identi ed miRNA with a certain length and on each position of all identi ed miRNA, respectively.
Hairpin structural characteristics of the miRNA precursor may be used to predict novel miRNA. The available software, miREvo [33] and mirdeep2 [31], were integrated to predict novel miRNAs by exploring their secondary structure, dicer cleavage site, and minimum free energy of the small RNA tags nonannotated in the prior steps above. The occurrence of miRNA families identi ed from the samples in other species was then explored. In this analysis pipeline, for a known miRNA, miFam.dat 22 (http://www.mirbase.org/ftp.shtml) was used to search for families; likewise, a novel miRNA precursor was submitted to Rfam 32.0 (http://pfam.xfam.org/) to identify Rfam families.
Identi cation of differentially expressed miRNAs.
The miRNA expression levels were quanti ed as transcripts per million (TPM) using these criteria [34,35]: normalization formula: normalized expression = mapped readcount/total reads*1000000. Differential expression analysis of the four treatment groups was performed in the R software environment using the DESeq package (v1.8.3), for which the P-values were adjusted using the Benjamini & Hochberg method [36]. A corrected P-value of 0.05 was set, by default, as the threshold for a statistically signi cant level of differential gene expression.
Prediction of miRNA target genes, GO and KEGG pathway analysis.
Predicting the target gene of a given miRNA was carried out using psRNATarget for plants (http://plantgrn.noble.org/psRNATarget/) [37]. GO enrichment was conducted to analyze the likely functions of genes targeted by their miRNAs using the GO enrichment tool in PlantRegMap (http://plantregmap.cbi.pku.edu.cn/) [38]. Pathway analysis was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/) [38], for which KOBAS (KEGG Orthology Based Annotation System) (http://kobas.cbi.pku.edu.cn) was used to test the statistical enrichment of the target gene candidates identi ed in the KEGG pathways [39].
Validating the expression of BLS-responsive miRNAs and their target genes by qRT-PCR.
The expression of 10 target genes was validated by qRT-PCR. Total RNA per sample was extracted with the TaKaRa MiniBEST Universal RNA Extraction Kit (Takara, Tokyo, Japan) and this RNA was then reverse-transcribed to cDNA using the PrimeScript 1st Strand cDNA Synthesis Kit (Takara, Japan). The qRT-PCR was run using the TB Green premix Ex Taq II Kit in a Roche LC480 in order to detect the expression of the target genes. Total small RNA per sample was extracted using RNAiso for Small RNA (Takara, Tokyo, Japan), and this RNA was then reverse-transcribed into cDNA and validated using the Mir-X™ miRNA First-Strand Synthesis and TB Green™ qRT-PCR User Manual (Takara, Japan) in a Roche LC480. The primers used for this are listed in Table 4; here, TC139057 and 5SrRNA served as internal controls. The cycling parameters were 95 °C for 30 s (1 cycle); 95 °C for 5 s, 60 °C for 30 s (40 cycles); 95 °C for 5 s, 60 °C for 1 min (1 cycle); 40 °C for 5 s (1 cycle). Gene expression was calculated using the Results Difference in BLS resistance between Kunlun14 and Z1141 barley.
The proportion of total leaf area that was diseased by BLS in KL14 was 12.80%, while it was nearly four times that in Z1141 (47.43%) ( Fig. 1A and B). The proportion of eld area (per square meter) diseased by BLS of Kunlun14 was 2.67%, while it was 92.67% for Z1141, based on eld experiments conducted from 2015 to 2017 ( Fig. 1C and D). The percentages differed signi cantly between the two barley cultivars (P < 0.01), con rming that Kunlun14 is a BLS-resistant cultivar, whereas Z1141 is BLS-susceptible.
Basic statistics of the small RNA data from Tibetan hulless barley libraries.
To identify BLS-related miRNAs in Tibetan hulless barley, small RNA libraries were constructed for Z1141 with (ZS) and without (ZN) BLS, and Kunlun14 with (KL14S) and without (KL14N) BLS. Totals of 15,141,660, 14,167,060, 13,591,121, and 12,992,954 raw reads were respectively generated from ZS, ZN, KL14S, and KL14N by high-throughput sequencing (Table 1). After removing the adaptors and junk reads, this left corresponding totals of 14,608,407, 13,573,623, 13,113,086, and 12,567,992 clean reads. Sequences within the range of 18-30 nt were then further analyzed. The length distributions from the libraries tested were similar overall, with most of the small RNAs being 19 to 24 nt in length. Among them, small RNAs with 20 nt of ZN were the most abundant, comprising 16.14% of the total reads (Fig. 2).

Identi cation of conserved and novel miRNAs in Tibetan hulless barley.
To identify conserved miRNAs from the small RNA libraries tested, reads with lengths of 18-30 nt were mapped onto H. vulgare-conserved miRNAs in the miRBase 22. This detected 36 conserved miRNAs in all the libraries tested, which belonged to 29 miRNA families (Table 2). Different families thus had different numbers of miRNAs. The family of hvu-miR5049 had the most members, with ve miRNAs. The expression levels of different miRNAs differed signi cantly, with read numbers ranging from zero to several thousand. The 10 most highly-expressed miRNAs were hvu-miR159b, hvu-miR156a, hvu-miR5048a, hvu-miR5051, hvu-miR168, hvu-miR171, hvu-miR397a, hvu-miR171, hvu-miR5049e, and hvu-miR6198; their accumulated reads accounted for 99.50% of all total conserved reads. Among them, hvu-miR159b had the highest number of reads, with an average of 4476 per library, amounting to 68.81% of the total conserved reads (Fig. 3A).
Identi cation of miRNAs expressed in two Tibetan hulless barley genotypes infected with BLS.
Target gene prediction and analysis of 24 BLS-responsive miRNAs.
The psRNATarget server predicted 546 target genes of BLS-responsive miRNAs in Tibetan hulless barley that were identi able for the 24 miRNAs (Additional le 6). Predicted target genes were classi ed and enriched via gene ontology (GO) and KEGG analyses. GO was used to describe the functional categories of the annotated unigenes and classify the transcripts with known annotated proteins. P ≤ 0.05 was used as the threshold to determine signi cantly enriched GO terms. GO classi cation based on sequence homology revealed that 427 out of all the assembled genes were categorized (Additional le 7). The top 30 enriched GO terms are indicated in Fig. 6, which shows that the target genes of the BLS-responsive miRNAs were included in all three categories: biological process, cellular component, and molecular function. Speci cally, the two most enriched categories were heterocyclic compound binding (207 genes) and organic cyclic compound binding (207 genes). Also, a term of copper ion binding (GO:0005507) was signi cantly enriched in 14 genes, most of which were laccase gene family members (Additional le 8). This indicated that laccase genes may play an important role in the response of Tibetan hulless barley to BLS infection.
A pathway-based analysis was performed using the KEGG pathway database to further elucidate the biological functions and interactions of the gene. We analyzed 546 target genes predicted by the KEGG database. A total of 75 genes were assigned to 64 pathways (Additional le 9). The three main pathways with the highest number of genes were glycerophospholipid metabolism, amino sugar and nucleotide sugar metabolism, and starch and sucrose metabolism. Top 20 enriched KEGG pathways are indicated in Fig. 7, one pathway of which was plant-pathogen interaction, which included four genes (Gene ID: HORVU4Hr1G051540, HORVU3Hr1G079210, HORVU7Hr1G008170, and HORVU7Hr1G055790). These annotations provide valuable information for investigating the disease-resistant pathway involved in BLS.
Expression levels of candidate miRNAs and their target genes in the barley response to BLS.
To validate the differential expression patterns of the BLS-responsive miRNAs, eight miRNAs and 10 target genes were randomly selected and their expression patterns were analyzed by quantitative realtime PCR (qRT-PCR) (Fig. 8). It was observed that the expression patterns of eight miRNAs were correlated with results obtained in the small RNA sequencing (RNA-Seq) study (Additional le 4). We also validated the expression of 10 target genes identi ed as BLS-responsive in this study. It was observed that the expression patterns of 10 target genes were correlated with results obtained in the RNA-Seq study (Additional le 4). Comparison of the expression patterns of miRNAs with their respective target genes showed that the steady-state levels of the miRNA transcripts were inversely related to the expression patterns of their target genes. For instance, the expression of the target gene (Gene ID: HORVU4Hr1G027270) was observed to be downregulated in contrast to the expression of the target miR168-3p, which was upregulated in the leaves after infection with BLS (Fig. 8A). On the contrary, the expression of the target genes (Gene ID: HORVU3Hr1G067760, HORVU4Hr1G088470, HORVU1Hr1G090190, HORVU7Hr1G083270, HORVU4Hr1G011590, HORVU5Hr1G096940, HORVU5Hr1G096950, and HORVU7Hr1G073050) was observed to be upregulated in contrast to the expression of their target miRNAs (hvu-miR171-5p, hvu-miR159b, hvu-novel-91, hvu-novel-46, hvu-novel-52, and hvu-novel-11) that were downregulated in the leaves after infection with BLS ( Fig. 8C-J).
The miR168 of barley was rst isolated in 2010 [45], and it was later isolated from rice [46], barley [47], and cucumber [48] under drought, cadmium, aluminum, and salt stress. A promoter analysis using transgenic A. thaliana plants revealed the transcriptional activation of miR168 by fungal elicitors [49], which was upregulated to modulate leaf morphology and plant growth, signi ed its positive role during leaf curl virus in tomato [26]. In this study, the predicted target gene of hvu-miR168-3p was IBR3 (Gene ID: HORVU4Hr1G027270), which was located on chromosome 4H and annotated to encode an acyl-CoA dehydrogenase of A. thaliana in Swiss-Prot. In NCBI, HORVU4Hr1G027270 was annotated to encode a clone HV_Mba184-F11 of H. vulgare (GenBank accession no. AC256280.1) and LZ-NBS-LRR class RGA of Aegilops tauschii (GenBank accession no. AF446141.1). Nucleotide sequence consistency of them was 90.61% and 78.98%, respectively. Plants evolved a complex innate immune system that recognized speci c pathogens via a speci c group of gene products called resistance (R) proteins [50]. RGAs were involved in plant disease resistance [51,52]. Thus, RGA appears to be crucial to the response of plant disease resistance. In this study, the expression of miR168-3p and its target gene HORVU4Hr1G027270 had opposite trends in the leaves after infection with BLS, which suggested that HORVU4Hr1G027270 might be regulated by hvu-miR168-3p during BLS-resistant of Tibetan hulless barley.
Kantar et al. suggested that miRNA156a in barley leaves responded to dehydration stress and played an important role in stress adaptation [45] (Kantar et al. 2010). miR156a targeted the squamosa promoterbinding protein-like genes, which played an important role in plant growth, development and abiotic stress [53]. In citrus, the overexpression of miR156a or knockout of its target genes SPL3 and SPL14 signi cantly increased the activity of wild kumquat callus [54]. In this study, the predicted target gene of hvu-miR156a was LIN (Gene ID: HORVU2Hr1G104640), which was located on chromosome 2H and annotated to encode an E3 ubiquitin-protein ligase of Medicago truncatula in Swiss-Prot. In NCBI, HORVU2Hr1G104640 was annotated to encode a predicted protein (NIASHv2044K07) of H. vulgare (GenBank accession no. AK366598.1), an E3 ubiquitin-protein ligase LIN-2 of Zea mays, and an E3 ubiquitin-protein ligase LIN-1 of Brachypodium distachyon (GenBank accession no. XM_010242010.3). Nucleotide sequence consistency of them was 99.95%, 84.88%, and 78.67%, respectively. In legumes, LIN controlled early infection by rhizobia and worked exclusively during rhizobial colonization [55]. In this study, the expression of hvu-miR156a and its target gene HORVU2Hr1G104640 had opposite trends, which suggested that HORVU2Hr1G104640 might be regulated by hvu-miR156a during BLS-resistant of Tibetan hulless barley. miR171 was reported in pear and tomato and proved to vital for plant growth and development [56,57]. In barley, overexpressing miR171 would lead to phase transitions and oral meristem determinacy [58]. In Arabidopsis, the expression of miR171 increased but the target gene AtSCL6-II decreased when both miR171 and AtSCL6-II were co-express transiently [59]. In addition, another study showed that miR171 overexpression altered the transition from vegetative to reproductive stage by activating the miR156 pathway [59]. Thus, the relationship between miR156 and miR171 in the regulation of BLS remains to be further studied. In this study, the predicted target gene of hvu-miR171-5p was DDB_G0268948 (Gene ID: HORVU3Hr1G067760), which was located on chromosome 3H and annotated to encode a methyltransferase of Dictyostelium discoideum in Swiss-Prot. In NCBI, HORVU3Hr1G067760 was annotated to encode a predicted protein (NIASHv2113O06) of H. vulgare (GenBank accession no. AK370649.1) and an S-adenosylmethionine-dependent methyltransferase of Aegilops speltoides (GenBank accession no. FJ236274.1) and Triticum urartu (GenBank accession no.FJ236273.1). Nucleotide sequence consistency of them was 99.73%, 93.76%, and 94.11%, respectively. S-adenosyl methionine (SAM), a widespread biological cofactor which was found in all aspects of life where it played very important roles in the transfer of methyl groups to DNA, proteins and small-molecule secondary metabolites [60]. The methylation process played important roles in various disease processes and industrial chemical processing [60]. SAM-dependent methylation affected mating through altering the pheromone signaling pathway of the maize smut pathogens (Ustilago maydis) [61]. The second target gene of miR171 was HORVU4Hr1G088470, which was located on chromosome 4H and annotated to encode a phosphatidylserine decarboxylase proenzyme 1 of Oryza sativa L. subsp. japonica in Swiss-Prot. In NCBI, HORVU4Hr1G088470 was annotated to encode a BAC 615K1 of H. vulgare (GenBank accession no. AY485643.1) and a phosphatidylserine decarboxylase (PSD) of Triticum monococcum (GenBank accession no. AY485644.1). Nucleotide sequence consistency of them was 100% and 94.44%, respectively. PSD activity in oat (Avena sativa) root cells was involved in the changes of phospholipid composition of cell membrane during domestication under drought stress [62]. In this study, the expression of hvu-miR171 and its target genes HORVU3Hr1G067760 and HORVU4Hr1G088470 had opposite trends in the leaves after infection with BLS, which suggested that HORVU3Hr1G067760 and HORVU4Hr1G088470 might be regulated by hvu-miR171-5p during BLS-resistant of Tibetan hulless barley.
In this study, the targets of three miRNAs (hvu-novel-91, hvu-novel-46 and hvu-miR159b) were WRKYs. WRKY proteins was a large family of transcription factors in plants, which was involved in the regulation of plant development, aging and stress resistance [63,64]. In rice, WRKY genes were involved in the resistance of the rice blast fungus Pyricularia oryzae [65]. Yokotani et al. showed that an OsWRKY24 gene was upregulated in response to rice blast [66]. In Arabidopsis, the overexpression of WRKY45 or WRKY1 enhanced disease resistance and drought tolerance [44,67]. In wheat, the overexpression of WRKY45 enhanced disease resistance to multiple fungi [68]. In Swiss-Prot, the target of hvu-novel-91 was a WRKY (Gene ID: HORVU7Hr1G083270), which was located on chromosome 7H and annotated to encode a WRKY46 transcription factor of A. thaliana. The target of hvu-novel-46 was a WRKY (Gene ID: HORVU4Hr1G011590), which was located on chromosome 4H and annotated to encode a WRKY1 transcription factor of Zea mays. Furthermore, the target of hvu-miR159b was also a WRKY (Gene ID: HORVU1Hr1G090190), which was located on chromosome 1H and annotated to encode a WRKY53 transcription factor of A. thaliana. In NCBI, HORVU7Hr1G083270 was annotated to encode WRKY transcription factor 146 in Triticum (GenBank accession no. MF770640.1), WRKY transcription factor 70 in Aegilops (GenBank accession no. XM_020292626.1), and WRKY transcription factor 11 in Aegilops (GenBank accession no. EU665440.1), but a predicted protein in H. vulgare (AK360029.1). The nucleotide sequence consistency of these was 94.83%, 94.83%, 94.83%, and 87.96%, respectively. HORVU4Hr1G011590 was annotated to encode four predicted proteins (NIASHv2089N03, NIASHv1011A24, NIASHv3105H03, and FLbaf169m16) in H. vulgare (GenBank accession no. AK369358.1, AK354776.1, AK375802.1, and AK248640.1) and WRKY transcription factor 1 in Aegilops (GenBank accession no. XM_020293740.1), the nucleotide sequence consistency of which was 99.09%, 98.47%, 98.77%, 98.77%, and 89.94%, respectively. HORVU1Hr1G090190 was annotated to encode WRKY transcription factor 1 in H. vulgare (GenBank accession no. DQ863110.1) and WRKY transcription factor 55 in Aegilops (GenBank accession no. XM_020341128.1). Nucleotide sequence consistency of them was 100.00% and 87.077%, respectively. Research showed that miR159 played an important role in disease symptom induction by a severe strain of Cucumber mosaic viru [69]. The inhibition of Fny-CMVinduced symptoms in Arabidopsis containing mutant alleles for the targets MYB33 and MYB65 of miR159ab con rmed the importance of miR159 in pathogenesis [69]. Phukan et al. showed WRKY included the WRKYGQK motif followed by a Cx4 5Cx22 23HxH or Cx7Cx23HxC zinc-nger motif [70]. In this study, the conserved domain results revealed that HORVU7Hr1G083270, HORVU4Hr1G011590, and HORVU1Hr1G090190 possess WRKY DNA -binding domains that belong to WRKY transcription factors. Most of the WRKY genes studied were found to be involved in the disease-resistant of plants. In apple, for example, MdWRKY26 can induce eight pathogenesis-related (PR) genes to enhance host resistance to leaf spot disease [64]; MdWRKY31 regulates plant resistance to Botryosphaeria dothidea through the salicylic acid signaling pathway by interacting with MdHIR4 [71]. In Arabidopsis, the overexpression of AtWRKY33 accelerated and intensi ed resistance to necrotrophic fungi [72]. These ndings strongly suggested that WRKY functions critically in plant resistance to fungal disease. In our study, the expression of hvu-novel-91 was downregulated, while the expression of WRKYs was upregulated in the leaves after infection with BLS. The expression of miRNAs and their target genes had opposite trends in the leaves after infection with BLS, which suggested that HORVU7Hr1G083270, HORVU4Hr1G011590, and HORVU1Hr1G090190 might be regulated by hvu-novel-91, hvu-novel-46m, and hvu-miR159b during BLS-resistant of Tibetan hulless barley.
Additionally, the conserved domain results showed that HORVU5Hr1G096940 and HORVU5Hr1G096950, possessing cytochrome P450 domains, belong to the CYP450 family. In rice, a CYP-like protein encoded by Os08g01480 helps plants to resist abiotic stresses [73]. In Arabidopsis, the enzymes encoded by CYP94B3 and CYP94C1 could promote the decomposition and inactivation of the plant hormone jasmonyl isoleucine complex, which affects the jasmonic acid metabolism pathway and thus affects the expression of stress-resistant genes [74]. In this study, the expression of hvu-novel-52 was downregulated, while the expression of HORVU5Hr1G096940 and HORVU5Hr1G096950 was upregulated in the leaves after infection with BLS, which suggested that HORVU5Hr1G096940 and HORVU5Hr1G096950 might be involved in the BLS disease process via the activity of hvu-novel-52 during BLS-resistant of Tibetan hulless barley.
Hvu-novel-11 targeted HORVU7Hr1G073050, which was located on chromosome 7H and annotated to encode the external alternative NAD(P)H-ubiquinone oxidoreductase B3 (NDB3) protein of in A. thaliana. In NCBI, HORVU7Hr1G073050 was annotated to a predicted protein (NIASHv2009M23) of H. vulgare (GenBank accession no. AK362739.1) and an external alternative NAD(P)H-ubiquinone oxidoreductase B2 (NDB2) of Aegilops (GenBank accession no. XM_020292911.1). Nucleotide sequence consistency of them was 99.80% and 93.42%, respectively. A series of treatment groups representing speci c AOX and NDH genes were identi ed in response to mitochondrial inhibition, plastid inhibition and other abiotic stresses [75]. In this study, the expression of hvu-novel-11 was downregulated, while that of HORVU7Hr1G073050 was upregulated in the leaves after infection with BLS, and the expression of KL14S was signi cantly higher than that of ZS. This result suggested that HORVU7Hr1G073050 might be involved in the BLS disease process via the activity of hvu-novel-11 during BLS-resistant of Tibetan hulless barley.

Conclusion
In summary, 4 known miRNAs and 4 novel miRNAs were revealed in this empirical study, the potential roles of which in BLS resistance might involve the regulation of downstream target genes. The expression levels of 8 miRNAs and 10 target genes were tested, suggesting that these miRNAs and their target genes jointly played a vital role in the defense response of Tibetan hulless barley to BLS, although their relationship in the process of resistance to BLS requires further investigation and elucidation. Fourteen laccase genes from GO analysis and four plant-pathogen interaction genes from the KEGG analysis were deleted because the expression of these genes did not change signi cantly in the result of RNA-Seq; their function requires further study in the future. This study may improve our understanding of miRNA