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 field area (per square meter) diseased by BLS of Kunlun14 was 2.67%, while it was 92.67% for Z1141, based on field experiments conducted from 2015 to 2017 (Fig. 1C and D). The percentages differed significantly between the two barley cultivars (P < 0.01), confirming 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).
Identification 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 five miRNAs. The expression levels of different miRNAs differed significantly, 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).
After prediction, a total of 56 novel miRNAs were obtained (Table 3). The 10 most highly expressed novel miRNAs were hvu-novel-9, novel-1, hvu-novel-2, hvu-novel-16, hvu-novel-11, hvu-novel-46, hvu-novel-52, hvu-novel-7, hvu-novel-91, and hvu- novel-18, and their accumulative reads accounted for 96.21% of the total conserved reads. Among them, hvu-novel-2 had the highest number of reads, with an average of 68.81 reads in each library, which accounted for 34.31% of the total conserved reads (Fig. 3B).
Identification of miRNAs expressed in two Tibetan hulless barley genotypes infected with BLS.
Expression profiling analysis revealed that the expression of many miRNAs was modulated by BLS. To identify the differentially expressed miRNAs, Venn diagram analysis was used to predict both conserved and novel miRNAs. There were 28, 27, 19, and 22 conserved miRNAs detected in the libraries of KL14N, KL14S, ZN, and ZS, respectively; 12 of these miRNAs were common to the four groups (Fig. 4A; Additional file 1). There were 54, 28, 50, and 45 novel miRNAs detected in the libraries of KL14N, KL14S, ZN, and ZS, respectively. Of these, 25 were common to the four groups (Fig. 4B; Additional file 2). A total of 24 miRNAs differed significantly between KL14N vs. KL14S and ZN vs. ZS (P < 0.05) (Fig. 4C; Additional file 3), of which there were 11 conserved miRNAs and 13 novel miRNAs. The expression levels of 24 miRNAs are displayed in a heatmap (Fig. 5; Additional file 4). Based on their fold-change values, 20 miRNAs (including 10 up-regulated miRNAs and 10 down-regulated miRNAs) were found between KL14S and KL14N (Additional file 5), and eight miRNAs were detected between ZS and ZN (six up-regulated miRNAs and two down-regulated miRNAs; Additional file 5). Ultimately, 24 differentially expressed miRNAs were found, of which 11 were conserved miRNAs (hvu-miR6199, hvu-miR168-5p, hvu-miR6198, hvu-miR5049c, hvu-miR171-5p, hvu-miR156a, hvu-miR6181, hvu-miR397a, hvu-miR6195, hvu-miR168-3p, and hvu-miR159b) and 13 were novel miRNAs (hvu-novel-43, hvu-novel-119, hvu-novel-1, hvu-novel-45, hvu-novel-20, hvu-novel-111, hvu-novel-52, hvu-novel-46, hvu-novel-79, hvu-novel-94, hvu-novel-9, hvu-novel-11, and hvu-novel-91) (Fig. 4C). The differing miRNA expression patterns emphasized the uniqueness of the various miRNA responses to BLS; hence, this suggests that BLS differentially affects miRNA expression in Tibetan hulless barley cultivars with differing disease resistance.
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 identifiable for the 24 miRNAs (Additional file 6). Predicted target genes were classified 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 significantly enriched GO terms. GO classification based on sequence homology revealed that 427 out of all the assembled genes were categorized (Additional file 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. Specifically, 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 significantly enriched in 14 genes, most of which were laccase gene family members (Additional file 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 file 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.
Transcriptome analysis showed that a total of 131 genes were significantly changed in Kunlun14 and Z1141 under BLS infection (Additional file 10). Among these genes, we concentrated on the functional genes and transcription factors possibly related to disease-resistant traits according to conserved domain analysis (http://pfam.xfam.org/search) and gene annotation at Swiss-Prot (http://www.gpmaw.com/html/swiss-prot.html), Uniprot (https://www.uniprot.org/), and NCBI (https://www.uniprot.org/). Additionally, considering the negative regulatory relationship between miRNA and target genes and that the relative transcript levels were greater than 20 in one library of KL14N, KL14S, ZN, or ZS (according to the results of Additional file 10), we found seven functional genes and three transcription factors associated with disease resistance. The reference genome sequence information of the candidate genes is shown in Additional file 11. Two members of the P450 genes (Gene ID: HORVU5Hr1G096940 and HORVU5Hr1G096950) were found. Three members of the WRKY transcription factor genes (Gene ID: HORVU1Hr1G090190, HORVU4Hr1G011590, and HORVU7Hr1G083270) were found. RGA gene (Gene ID: HORVU4Hr1G027270), LIN (Gene ID: HORVU2Hr1G104640), SAM (Gene ID: HORVU3Hr1G067760), PSD (Gene ID: HORVU4Hr1G088470), and NDB (Gene ID: HORVU7Hr1G073050) were found. The miRNAs of these target genes were hvu-miR168-3p, hvu-miR156a, hvu-miR171-5p, hvu-miR159b, hvu-novel-91, hvu-novel-46, hvu-novel-52, and hvu-novel-11.
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 real-time 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 file 4). We also validated the expression of 10 target genes identified 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 file 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).