Tillering characterization of tall fescue
The two tall fescue genotypes (Ch-3 and Ch-5) showed a remarkable difference in tillering development (Fig. 1). Ch-3 with high tiller production rate had 122 tillers, but Ch-5 had only 9 tillers after 2 months of the establishment (Fig. 1A, B). Furthermore, we found that Ch-3 plants were still undergoing vegetative growth while Ch-5 plants have begun heading growth (Fig. 1A). Therefore, the difference in tiller number between the two tall fescue genotypes reached its maximum during this growth stage, which was the most appropriate period to collect for sRNAs sequencing. But, the tall fescue ‘Houndog V’ plants showed different growth stage (Fig. 1C and D.) At the 2.5-leaf stage, plants were at the vegetative growth and the tiller number was zero, which named Pre-tillering. However, at the 4.5-leaf stage, plants began to tiller and the tiller number attained two or three, which named Tillering. Hence, in order to identify the miRNAs involved in tillering more accurately, samples were collected both at vegetative growth and tillering stages.
Deep sequencing of sRNAs
To identify and characterize the role of miRNAs in the tillering development of tall fescue, the samples from two genotypes Ch-3 and Ch-5 with different tillering, and the Pre-tillering and Tillering samples from ‘Houndog V’ at the different tillering development stage were used for sRNAs libraries construction and sequencing. Using the Illumina HiSeq 2500 system, a total of 222 million high-quality clean reads were obtained and each sample produced more than 18 million clean reads averagely. By choosing 18-30 nt clean reads, about 55.97% (5 075 789), 43.4% (5 814 258), 47.9% (7 516 722) and 47.7% (7 501 722) of the total reads were perfectly mapped to the full-length transcriptome data of tall fescue for Ch-3, Ch-5, Pre-tillering and Tillering, respectively (Additional file 1). Through scanning the size distribution based on total unique reads of the four samples, 24-nt sRNAs were the most abundant in all samples, accounting 36.8%, 36.2%, 35.33%, and 35.94 in Ch-3, Ch-5, Pre-tillering and Tillering, respectively (Additional file 2). The total mapped reads were grouped into various noncoding sRNA categories, including known miRNA (0.16-0.43%), rRNAs (12.87–24.14%), tRNAs (0%), snRNAs (0.13–0.49%), snoRNAs (0.13–0.25%) and unmatched sRNAs (75.28–84.4%) as shown in Table 1. The total reads that could be annotated to known miRNAs were 9 079, 28 516, 23 705, and 28 673 for Ch-3, Ch-5, Tillering, Pre-tillering, respectively. Interestingly, Ch-3 with high tiller production rate showed less number of known miRNAs than that of Ch-5. Consistent with this result, the Tillering samples had fewer known miRNAs compared with the Pre-tillering samples.
Table 1
Distribution of unique reads among different categories in tall fescue.
Category
|
Ch-3
|
Ch-5
|
Tillering
|
Pre-tillering
|
Total
|
5075789 (100%)
|
5814258 (100%)
|
7516722 (100%)
|
7501722 (100%)
|
KnownmiRNA
|
9079 (0.16%)
|
28516 (0.43%)
|
23705 (0.34%)
|
28673 (0.40%)
|
rRNA
|
1236690 (24.14%)
|
1034765 (17.95%)
|
1136985 (14.87%)
|
1210625 (16.10%)
|
tRNA
|
1 (0.00%)
|
1 (0.00%)
|
3 (0.00%)
|
2 (0.00%)
|
snRNA
|
12537 (0.27%)
|
27204 (0.49%)
|
10006 (0.13%)
|
10312 (0.14%)
|
snoRNA
|
7674 (0.15%)
|
7637 (0.13%)
|
18656 (0.25%)
|
16876 (0.22%)
|
Other
|
3809809 (75.28%)
|
4716135(81.00%)
|
6327366 (84.40%)
|
6235234 (83.14%)
|
Identification of known and novel miRNAs
To identify known and novel miRNAs in tall fescue, all the unannotated unique reads that perfectly mapped to tall fescue full-length transcriptome data were aligned to plant miRNAs in miRBase21 database. A total of 3 051 miRNA precursors were identified and then compared with mature miRNAs from Brachypodium distachyon, Oryza sativa, Zea mays, Hordeum vulgare, Aegilops tauschii and Festuca arundinacea. Accordingly, 148 known miRNAs belonging to 70 families were identified (Additional file 3, 4). Among them, 60 families were well conserved and found in more than two plant species. Especially, miR395, miR399, miR169, miR166, miR160, miR172 and miR395 were highly conserved and identified in nearly 30 plant species (Additional file 4). Besides, we analyzed the first base preference for known mature miRNA (Fig. 2). Among these four groups: Ch-3, Ch-5, Pre-tillering and Tillering showed preference towards U for the first base in 18~23-nt sRNAs, but the first base preference in 24~30-nt sRNAs was different between the four groups. Ch-3 and Tillering had more U preference than Ch-5 and Pre-tillering in 24~35-nt sRNAs. On the other hand, 10 families such as miR5048, miR5185, miR9863, miR7708, etc. were less conserved and found only in one plant species. In addition to these known miRNAs, 60 novel miRNAs were identified by predicting the hairpin structures of their precursor sequences (Additional file 5), and the abundance of U was decreased for the first base preference in 18~30-nt sRNAs for the identified novel miRNAs (Additional file 6).
Differential expression of miRNAs during tillering development
To identify the genome-wide small RNAs involved in tillering development, the different expression levels of known and novel miRNAs were clustered using log10 (TPM+1). The miRNAs detected in the two genotypes or different tillering development stages were separated (Fig. 3; Additional file 7), indicating that the tiller development was regulated by miRNAs. Based on the read count value, a total of 34 miRNAs from 29 miRNA families exhibited differential accumulation between Ch-3 and Ch-5 genotypes (Additional file 8). We found 24 significantly up-regulated and 10 down-regulated miRNAs in Ch-3 when compared to Ch-5. Furthermore, the abundance of various miRNAs at the tiller development stage was altered. In tillering plants, the expression levels of 16 miRNAs were increased, while 13 miRNAs were decreased compared with the non-tillering plants (Additional file 9).
To identify the key miRNAs controlling tillering in tall fescue, we performed the Venn diagram analysis and measured 14 miRNAs co-up-regulated and 4 miRNAs co-down-regulated at both Ch-3/Ch-5 and Tillering/Pre-tillering groups (Additional file 10). The differentially co-expressed miRNAs were 5 novel and 13 known miRNAs belonging to 15 families (Table 2). Furthermore, in order to check whether the 18 miRNAs is relative to tall fescue tillering, we measured the expression level of six randomly picked miRNAs from them using Stem-Loop qRT-PCR. As shown in Additional file 11, bdi-miR160f, novel_22, novel_23, osa-miR156a, osa-miR408-3p and osa-miR394 showed the same change pattern in Ch3 and Ch5 plants compared with the sequencing data. In addition, novel_22, novel_23, osa-miR156a, osa-miR408-3p and osa-miR394 showed higher transcriptional level in axillary bud than that in out growth bud samples, and they also displayed more expression level in Ch3 than that in Ch5. The bdi-miR160f showed lower transcriptional level in Ch3 than that in Ch5, and it was lower in axillary bud than that in out growth bud samples. The finding indicates that the 18 miRNAs may play key roles in tall fescue tillering development.
Table 2
Co-up-regulated or co-down-regulated miRNAs at both Ch-3/Ch-5 and Tillering/Pre-tillering groups.
miRNA
|
CH3
|
CH5
|
Tillering
|
Pre-tillering
|
padj
|
Target genes number
|
Mature sequence
|
ata-miR5168-5p
|
38.87
|
0.00
|
16.10
|
4.45
|
***
|
26
|
GGGUUGUUGUCUGGUUCAAGG
|
bdi-miR160f
|
5.57
|
27.78
|
3.11
|
66.33
|
***
|
17
|
UGCCUGGCUCCCUGUAUGCC
|
bdi-miR167e-3p
|
407.37
|
207.92
|
11.23
|
5.00
|
***
|
9
|
AGGUCAUGCUGGAGUUUCAUC
|
bdi-miR397b-5p
|
106.33
|
12.46
|
81.01
|
56.33
|
***
|
13
|
AUUGAGUGCAGCGUUGAUGAA
|
hvu-miR397a
|
117.00
|
18.92
|
157.90
|
87.78
|
***
|
4
|
CCGUUGAGUGCAGCGUUGAUG
|
novel_13
|
11.59
|
68.34
|
17.12
|
306.91
|
***
|
2
|
AUUGUGACUUACAAACUAGGACGG
|
novel_18
|
0.91
|
13.89
|
4.89
|
26.46
|
***
|
6
|
AUUGUGCACUAACCCUCUAAUGAG
|
novel_2
|
432.73
|
43.29
|
1869.22
|
762.40
|
***
|
16
|
UGAAGUGUUUGGGGGAACUC
|
novel_22
|
16.67
|
0.00
|
127.31
|
61.78
|
***
|
31
|
UUUUUCUGUAGCUCAAAAUAC
|
novel_23
|
14.08
|
0.73
|
133.74
|
62.30
|
***
|
2
|
UAGAGAUGUGGCAUUACUUUGAGU
|
osa-miR156a
|
309.58
|
86.59
|
1795.78
|
844.18
|
***
|
17
|
UGACAGAAGAGAGUGAGCAC
|
osa-miR394
|
154.64
|
20.11
|
1442.09
|
678.90
|
***
|
10
|
UUGGCAUUCUGUCCACCUCC
|
osa-miR397a
|
103.28
|
18.10
|
139.14
|
83.35
|
***
|
6
|
UCAUUGAGUGCAGCGUUGAUG
|
osa-miR397b
|
117.00
|
18.92
|
157.39
|
87.37
|
***
|
8
|
UUAUUGAGUGCAGCGUUGAUG
|
osa-miR408-3p
|
98.37
|
25.53
|
274.00
|
152.32
|
***
|
6
|
CUGCACUGCCUCUUCCCUGGC
|
osa-miR444b.2
|
33.42
|
0.60
|
433.59
|
177.03
|
***
|
24
|
UGCAGUUGUUGUCUCAAGCUU
|
zma-miR1432-5p
|
1.87
|
38.38
|
2.09
|
105.41
|
***
|
9
|
CUCAGGAGAGAUGACACCGAC
|
zma-miR528a-3p
|
7.36
|
0.00
|
38.00
|
0.00
|
***
|
6
|
CCUGUGCCUGCCUCUUCCAUU
|
Note:Thered indicates the miRNA is positive with tillering, and the green indicates the miRNA is negative with tillering. Transcript abundance of miRNAs in Ch-3, Ch-5, Tillering and Pre-tillering was present with readcount value. *** indicates the padj |
Prediction and enrichment of miRNA target genes
To gain additional insights into the miRNA-genes pathway that related to tiller development in tall fescue, we obtained a total of 28 927 potential target genes for all known and novel miRNAs from Ch-3, Ch-5, Tillering and Pre-tillering samples (Additional file 12). The number of target genes for each miRNA showed a remarkable difference, such as bdi-miR5049-3p and bdi-miR5067 had only one target gene but bdi-miR845 and osa-miR414 targeted more than 30 genes. Here, miRNA regulated targeted genes through two miRNA actions: target messenger RNA cleavage and translational repression. GO enrichment analysis was performed to evaluate the potential functions of target genes of differential expression miRNA. A total of 2 093 target genes of miRNA with differential accumulation between Ch-3 and Ch-5 were identified and categorized into 52 GO functional subcategories (Fig. 4; Additional file 13).
Results showed that the oxidation-reduction process in the biological process (BP) and membrane protein complex in the cellular component (CC) were the most enriched GO terms with 212 (10.4%) and 112 (5.5%) genes, respectively. Between Tillering and Pre-tillering, 1 140 target genes of miRNA were identified and categorized into 42 GO functional subcategories which were fewer than between Ch-3 and Ch-5. Molecular function (MF) showed the highest enrichment in Tillering and Pre-tillering, and oxidoreductase activity (32, 9.58%), coenzyme binding (19, 5.69%) and cofactor binding (21, 6.29%) were the most abundant GO terms, which were different with the dominant GO terms enriched in BP and CC in Ch-3 and Ch-5. However, we found some co-enriched target genes in both Ch-3/Ch-5 and Tillering/Pre-tillering groups using KEGG analysis (Fig. 5; Additional file 14). About 689 target genes were mapped into 20 KEGG pathways in the Ch-3/Ch-5 group, while 100 genes were mapped into 20 KEGG pathways in the Tillering/Pre-tillering group. Fortunately, five KEGG pathways such as “ubiquitin-mediated proteolysis”, “phagosome”, “fatty acid biosynthesis”, “oxidative phosphorylation” and “biosynthesis of unsaturated fatty acids” showed co-enrichment in both Ch-3/Ch-5 and Tillering/Pre-tillering groups. The most co-enrichment of miRNA target genes were detected in the “ubiquitin-mediated proteolysis” and “oxidative phosphorylation”, suggesting that these two pathways may play a great role in controlling tillering development in the grass.
Validation of miRNA expression patterns by qRT-PCR
To validate the small RNA sequencing results involved in tillering development, we analyzed the expression level of twelve miRNAs using Stem-Loop qRT-PCR. Twelve miRNAs were randomly selected including two novel miRNAs, seven conserved miRNAs and three non-conserved miRNAs. As shown in Fig. 6, miRNA novel-2, osa-miR444b.2 and bdi-miR397b-5p were co-up-regulated, while bdi-miR5067 showed co-down-regulated in the Ch-3 and Tillering grass compared with the Ch-5 and Pre-tillering grass, which was in line with miRNA sequencing results. The abundance of miRNA novel-2, osa-miR160a-5p and osa-miR171b were significantly increased in the Pre-tillering grass compared to Tillering ones, but the expression did not change between Ch-3 and Ch-5 grasses. However, bdi-miR167e-3p exhibited higher expression level but ata-miR156d-3p and osa-miR168a-5p showed lower expression levels in Ch-3 than that in Ch-5, while their expression level was not changed between Tillering and Pre-tillering grass. The transcript abundance of bdi-miR5181b had no difference between both tall fescue genotypes Ch-3 and Ch-5 and between Tillering and Pre-tillering grass. Interestingly, ata-miR9772a-5p showed increased transcript abundance in Ch-3 than that in Ch-5, but had a lower expression level in the Tillering grass compared with the Pre-tillering grass. Additionally, the correlation between small RNA sequencing and qRT-PCR results were evaluated in terms of transcript abundance. As shown in Additional file 15, the qRT-PCR measurements were moderately correlated with miRNA profiles revealed by Illumina sequencing (y = 0.0085x + 2.4609; R2= 0.5431; P < 0.001), indicating that the small RNA sequencing data was accurate and effective.