Abstract Background Pre-harvest sprouting (PHS) is one of the most important problems associated with the severe decrease of yield and quality under disaster weather of continuous rain in wheat harvesting stage. At present, the functions and mechanisms related to the involvement of post-transcriptional regulation has not been studied very clearly in PHS resistance. Results This study compared the differences of germinated seeds in miRNAome between the PHS-tolerant and PHS-susceptible white wheat varieties. A total of 1879 miRNAs were identified from three different stages during seed germination. In order to further obtain candidate miRNAs, the different datasets of differentially expressed miRNAs were excavated by using differential-expression and time-series analysis. Combined with degradome data, the miRNA-mRNA networks analysis was performed after genome-wide screening of target genes, and then KEGG enrichment highlighted that the starch and sucrose metabolism pathway related to PHS was specifically enriched in an especial target-gene dataset derived from R12R18-HE miRNAs. Based on transcriptome data, a network associated with starch metabolism was systematically and completely reconstructed in wheat. Then, the starch degradation pathway controlled by seven miRNA-RNA pairs were supposed to be the essential regulation center for seed germination in wheat, which also could play a critical role on the PHS resistance. Conclusion Our findings revealed the complex impact of the miRNA-mediated mechanism for forming intrinsic and inherent differences, which resulting in significant difference on PHS performance between white wheat varieties. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-024-06039-8. Keywords: Wheat, Pre-harvest sprouting, MiRNA, Seed germination, Starch Introduction The pre-harvest sprouting is a worldwide natural disaster. Periods of rain during pre-harvest could cause unharvested cereal crops to sprout, such as wheat, barley, sorghum, and rice. Resistance to PHS is a complex quantitative trait that is controlled by both genetic factors and external environmental factors [[34]1], including relative humidity and temperature, awn and spike structure, germination-inhibitory compounds, α-amylase (α-Amy) activity, grain structure and color, phytohormones, seed dormancy [[35]2]. Among them, the molecular mechanisms about seed dormancy and grain color have been studied more and more thoroughly. Seed dormancy is the main genetic component for PHS resistance and determines the resistance level of wheat. Therefore, the genetic regulation mechanism of seed dormancy is often taken as the main content when studying the mechanism of PHS resistance in previous reports [[36]3]. Several resistant genes such as MKK3, Vp1, PM19, MFT, PHS1, PHS-3D, ABI5, FUS3 and DOG1, were characterized from wheat for grain dormancy [[37]2]. Among them, there was controversy about the candidate genes for Phs-A1. Based on Transcriptomic analysis of wheat near-isogenic lines, the PM19-A1 and PM19-A2 were identified as the candidate genes located at Phs-A1 locus, the major dormancy QTL for regulating seed dormancy [[38]4]. Afterwards, TaMKK3-A were isolated as the real candidate for Phs-A1 by map-based cloning and found that it encoded a MAPKinase Kinase [[39]5]. White-grained wheat has been reported to be more susceptible to PHS than red-grained wheat over a hundred years ago. After the efforts of Lang et al, the dual role Myb10-D/PHS-3D which can regulate the concentrations of flavonoids and ABA, ultimately affecting both the grain color and seed germination, were revealed as the molecular basis of relationship between grain color and PHS resistance in wheat [[40]6]. Subsequent research has confirmed this result, from which CRISPR/Cas9-mediated restoration of Tamyb10 could convert a white wheat variety into a red one and improve its PHS tolerance [[41]7]. Nevertheless, there are still a few excellent gene resources which confer the resistance of white-grained wheat to pre-harvest sprouting in some special germplasms of wheat landraces. For example, a novel allele of viviparous-1 (Vp-1Bf) associated with high seed dormancy were reported to make the resistance of white wheat to PHS in a Chinese wheat landrace, Wanxianbaimaizi [[42]8]. Except the major QTL TaPHS1, two stable minor QTLs for resistance to pre-harvest sprouting (PHS) were identified from a white wheat variety “Danby”, which could greatly improve PHS resistance in white wheat [[43]9]. Plants make decisions throughout their lifetime based on complex networks, in which miRNAs play indispensable roles in various biological processes of plant growth and development, such as tissue differentiation and organ development, signal transduction, biological and abiotic stress, etc [[44]10]. miRNA-mediated regulation during seed dormancy, seed germination and seed vigor has been studied in different plants. Previous research revealed that in Arabidopsis, higher miRNA156 levels enhanced seed dormancy [[45]11]. However, miRNA156 levels increased in embryos reduced seed dormancy of Ginkgo biloba [[46]12], and in rice the mutation of miRNA156 lead to enhance seed dormancy [[47]13]. Sequestering miR165/166 enhances seed germination in Arabidopsis thaliana under normal condition and ABA treatment [[48]14]. Deep sequencing of small RNA reveals that miR9736-z, miR5059-z, miR399-y, and miR163-z etc. are the key miRNAs regulating seed germination-related genes in Arabidopsis [[49]15]. The integrated analyses of transcriptome and miRNAome data showed that the seed imbibition process is regulated by a variety of DEGs and microRNAs in mung bean, including 23 putative genes targeted by miR156, miR171b-3p, miR166e-3p, and miR169-1, etc [[50]16]. The research about miRNA-target pairs participating in hormone interaction in barley revealed that miR393-mediated auxin response regulation affects grain development and influences gibberellic acid and abscisic acid homeostasis during germination [[51]17]. By STRING database assay, a miRNA-mediated gene interaction network regulating seed vigor in rice was revealed, which comprised at least two miRNA-target pathways: the miR5075-mediated oxidoreductase related pathway, and the miR164e related pathway [[52]18]. 791 mature miRNAs, associated with improved maize seed vigor induced by gibberellin, were obtained with different expressions by comparing the miRNAs sequences in miRbase database [[53]19]. Further analysis also demonstrated that miRNAs play critical roles in transcriptional regulation of critical pathways during PHS by interacting with target genes [[54]20]. For instance, overexpression of miR9678 delayed germination and improved resistance to PHS in wheat through modulating the GA/ABA signaling [[55]21]. However, the other miRNA-mediated regulation mechanisms for seed germination or dormancy conferred resistance to PHS is not very well understood. The white wheat cultivars, which were more popular for millers, had higher flour yield than red-grained wheat cultivars at the same grade of whiteness of flour. The improvement of white-grained wheat with strong resistance to PHS is important for expanding its production in order to meet the market demands [[56]1]. Meanwhile, the molecular mechanism of small RNA related to PHS resistance in wheat also needs to be explored. Hence, dissecting the regulatory function of miRNAs in modulating seed dormancy and germination is very important for understanding the molecular mechanism of resistance to PHS in white wheat. Here, identification and comparison of important miRNAs associated with seed germination were performed between two wheat cultivars with different germination speeds to resolve the mechanisms of PHS resistance regulated by small RNA. Ultimately, we found out several key miRNA-mRNA interaction pairs participated in starch degradation affecting pre-harvest sprouting resistance of wheat, which will provide an important theoretical basis for pre-harvest sprouting resistance. In our future work, it will employ the key miRNAs through molecular breeding or transgenic breeding techniques for the genetic improvement and cultivation of PHS-resistant varieties, which is of great significance for ensuring human food security under the extreme climate environment. Materials and methods Plant material and analysis of pre-harvest sprouting resistance Two white wheat varieties HZ1005 (PHS-susceptible) and Huamai1168 (PHS-tolerant) were used in this study. The wheat variety "Huamai1168" was bred from offspring of "Chuan8910/Huaai01//Zhoumai12/Emai12" through pedigree selection for several years. "HZ1005" is a new wheat line developed through pedigree selection of "Zhoumai18/Zhengmai9023". They were planted in experimental fields of Huazhong Agricultural University. The identification of pre-harvest sprouting about HZ1005 and Huamai1168 was performed in two consecutive years from 2021 to 2022. In each year, when wheat went into the full bloom stage, the hang tags were used to mark the flowering spikes. After 35 days, thirty spikes of HZ1005 or Huamai1168 were selected and cut, then soaked in sterilized water for 1 hour after sterilization. Afterwards, the spikes were covered with absorbent papers, placed in a seed germination box, and incubated at 18℃. Spray water 3 times a day to ensure that the absorbent paper is moist. After 3 days, the spikes were threshed to check the germination of the seeds. The root reaching the length of the seed and the bud reaching half of the length of the seed were used as the criterion for seed germination. The whole spike germination rate (SGR) was used to evaluate the pre-harvest sprouting. Preharvest sprouting rate = number of germinated seeds/total number of seeds × 100%. T test was used for statistical analysis (∗P < 0.05; ∗∗P < 0.01; ∗∗P < 0.001). In order to further explore the difference between the speed of seed germination of two wheat varieties, the seed germination experiment was carried out to observe the difference of seed germination between HZ1005 and Huamai1168. The seed germination rates were counted at 24, 48, 72 hours. Small RNA library construction and sequencing The germinated seeds from HZ1005 and Huamai1168 were collected at different time points (12h, 18h and 36h) during germination. Each time point had three biological replicate samples. Fifteen consistently germinated seeds were selected from each sample, immediately frozen in liquid nitrogen and stored at − 80°C for RNA extraction. Total RNA was extracted by using TRIzol reagent (Invitrogen, Grand Island, NY, USA). The RNA amount and purity of each sample was quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed by Bioanalyzer 2100 (Agilent, CA, USA) with RIN number > 7.0, and confirmed by electrophoresis with denaturing agarose gel. The small RNA (sRNA) libraries were constructed with TruSeq small RNA sample prep kits (Illumina, San Diego, United States). Here, total 18 small RNA libraries (RG12-1, RG12-2, RG12-3, RG18-1, RG18-2, RG18-3, RG36-1, RG36-2, RG36-3, SG12-1, SG12-2, SG12-3, SG18-1, SG18-2, SG18-3, SG36-1, SG36-2 and SG36-3) were built. Then, all libraries were sequenced using Illumina HiSeq2500 (LC Sciences, Houston, Texas, United States). Bioinformatics analysis of sRNA and identification of miRNA Raw reads were subjected to an in-house program, ACGT101-miR (LC Sciences, Houston, Texas, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. Subsequently, unique sequences with length in 18 ~ 25 nucleotide were mapped to specific species precursors in miRBase 22.0 by BLAST search to identify known miRNAs and novel 3p- and 5p- derived miRNAs. Length variation at both 3’ and 5’ ends and one mismatch inside of the sequence were allowed in the alignment. The unique sequences mapping to specific species mature miRNAs in hairpin arms were identified as known miRNAs. The unique sequences mapping to the other arm of known specific species precursor hairpin opposite to the annotated mature miRNA-containing arm were considered as novel 5p- or 3p derived miRNA candidates. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in miRBase 22.0 by BLAST search, and the mapped pre-miRNAs were further blasted against the specific species genomes to determine their genomic locations. The above two we defined as known miRNAs. In addition, sRNA precursors containing classic hairpin structure, which were unannotated in the former steps, were used to predict novel miRNAs by miREvo and mirdeep2. The expression levels of miRNAs in each library were normalized. Data normalization followed the procedures as described in a previous study [[57]22]. Then, differential expressed miRNAs were screened using Student’s T-test or ANOVA based on the experiments design. The significance threshold was set to be 0.01 and 0.05 in each comparison groups. mRNA library construction and sequencing Total RNA was extracted by using TRIzol reagent (Invitrogen, Grand Island, NY, USA). High purity 1.5 µl of total RNA from each sample was reverse-transcribed to construct cDNA libraries according to the instructions of the Illumina TruSeq RNA Sample Preparation Kit (Illumina, United States). Total 18 mRNA libraries (RG12-1, RG12-2, RG12-3, RG18-1, RG18-2, RG18-3, RG36-1, RG36-2, RG36-3, SG12-1, SG12-2, SG12-3, SG18-1, SG18-2, SG18-3, SG36-1, SG36-2 and SG36-3) were constructed. Finally, paired-end sequencing was performed on an Illumina Novaseq™ 6,000 (LC Sciences, United States). Mapping and transcriptome assembly For transcriptome sequencing, the raw data were processed by SOAP for mapping reads to the ribosome database. Then the reads that had been mapped to the ribosome database were removed. Subsequently, the clean reads were obtained by removing missing bases and low-quality bases (base quality ≤ 10). High-quality clean data were mapped to the reference genome (Chinese Spring IWGSC RefSeq v.1.0) using hisat2. After the transcriptome was generated, StringTie and ballgown ([58]http://www.bioconductor.org/packages/release/bioc/html/ballgown.ht ml) were used to estimate the expression levels of all transcripts and perform expression levels for mRNAs by calculating FPKM (Fragments Per Kilobase of exon model per Million mapped fragments). Degradome library construction and target identification For the degradome sequencing, equal amounts of all RG (RG12-1, RG12-2, RG12-3, RG18-1, RG18-2, RG18-3, RG36-1, RG36-2, RG36-3) or SG (SG12-1, SG12-2, SG12-3, SG18-1, SG18-2, SG18-3, SG36-1, SG36-2, SG36-3) RNA samples used in small RNA sequencing were mixed together to generate two degradome libraries for searching the potential target of wheat miRNAs. The cDNA library was constructed and sequenced on Illumina Hiseq 2500 (LC Sciences, United States). Raw sequencing reads were obtained using Illumina’s Pipeline v1.5 software. After that, reads were mapped to the wheat genome to obtain cDNA sense and antisense tags. The tags mapped to cDNA or mRNA sequences were then used to predict cleavage sites. The pure reads with 20 and 21 nucleotides were used to identify potentially cleaved targets by PAREsnip ([59]http://srna-workbench.cmp.uea.ac.uk/tools/paresnip/) and CleaveLand 3.0 ([60]http://sites.psu.edu/axtell). The targets predicted were classified into five categories (0, 1, 2, 3 and 4). Based on the expression characteristics of the wheat transcriptome data, t-plots were built for the high-efficiency analysis of the potential miRNA targets. Finally, all candidate targets were used for function annotation and gene ontology (GO) analysis. GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis To systematically analyze and identify potential biological pathways, the OmicShare cloud platform were used for GO and KEGG pathway enrichment analysis with the target genes of different miRNAs. GO enrichment analysis include three categories: biological process (BP), molecular functions (MF), and cell components. The GO enriched terms were identified with P values < 0.05. The KEGG pathway terms with P values < 0.05 were also identified and the Top 25 of KEGG enrichment were selected for further analysis. Additionally, we selected the KEGG pathway analysis to execute functional annotation on target gene datasets TGs1, TGs2, and TGs3 respectively. Quantitative real-time PCR The total RNA was reverse-transcribed to miRNAs using the Mir-X™ miRNA First-Strand Synthesis and SYBR Kit (Clontech, CA, USA). Mir-X miRNA qRT-PCR SYBR Kits use a single-step, single-tube reaction to produce first-strand cDNA, which is then specifically and quantitatively amplified using a miRNA-specific primer. The miRNA-specific 5’ primers of the miRNAs were designed according to the manufacturer’s instructions [[61]23]. U6 for miRNAs were used as internal control. The RNA was reverse transcribed to cDNA using the SuperScript III reverse transcriptase (Invitrogen). The reverse-transcription reaction was performed at 16°C for 30 min, followed by 60 cycles at 30°C for 30 sec, 42°C for 30 sec and 50°C for 1 sec. The reaction mixture was incubated at 85°C for 5 min to inactivate the reverse transcriptase. Then, quantitative real-time PCR (qRT-PCR) was performed using CFX96 Real Time System (Bio-Rad) to validate the sRNA-seq and RNA-seq data. Actin was used as an endogenous control. The reaction procedure was as follows: 10 min at 95°C, followed by 40 cycles of 15s at 95°C and 60s at 60°C. The reactions were performed with three biological replicates, and a melting curve analysis was carried out to ensure specific amplification of the products. The primers were listed in Table S1. Result Differences of preharvest sprouting resistance between HZ1005 and Huamai1168 HZ1005 and Huamai1168 are two white varieties with similar grain morphology of seeds harvested in fine weather (Fig. [62]1A). Under the continuous rainy condition before harvest, it was found that few seeds of wheat variety Huamai1168 appeared the sprouted phenomenon in comparison with the HZ1005 in the field (Fig. [63]1B, Fig. S1A). The preharvest germination tests of the single spike from the years of 2021 and 2022 also showed that the preharvest germination rates of HZ1005 (48.2% and 56.7%) were significantly higher than that of Huamai1168(1.8% and 7.2%) (Fig. [64]1C). The comparative test of germination rates was performed by using healthy seeds from both varieties. On the first day, the germination rates revealed the 4-fold and 2-fold changes between HZ1005 and Huamai1168 from 1 to 2 days, and it was up to 97% for HZ1005. On the third day, the germination rates of Huamai1168 were at above 96%, which was almost the same as HZ1005(Fig. [65]1D). It indicated that the seeds of these two varieties had the same germination vigor, but they had difference in the germination speed. To show the difference between the varieties of HZ1005 and Huamai1168 in germination speed, the method of visual observation and record temporal was carried out and listed in chronological order. Seeds were soaked in sterile water for 6 hours after full absorption and then arranged on filter papers. It was obvious that the seeds of HZ1005 began to break through the seed coats at 24h, but the seeds of Huamai1168 began to sprouting until nearly at 48h (Fig. [66]1E). Finally, all seeds can germinate and grow normally, and there was no significant difference in seedling growth after 144h (Fig. S1B). Therefore, the difference of germinated speed rather than the germination rate was the most important factor that caused the difference in resistance to the pre-harvest sprouting (PHS) between the two varieties. The germinated growth was slow for Huamai1168, which was resulting in its tolerance to PHS under the field condition. Therefore, the difference in germination speed, rather than the difference in germination rate, is an important factor contributing to the difference in PHS resistance between the two varieties. Accordingly, HZ1005 was defined as RG (rapid germination) variety, and Huamai1168 was defined as SG (slow germination) variety in this study. Fig. 1. [67]Fig. 1 [68]Open in a new tab Phenotype identification of pre-harvest sprouting (PHS) and seed germination for HZ1005 and Huamai1168. A Grain morphology of seeds harvested in fine weather. B PHS phenotypes of HZ1005 and Huamai1168. C Using whole spike germination test to evaluate PHS rate of HZ1005 and Huamai1168 in years of 2021 and 2022. Significance is determined by Student’s t-test, ***P < 0.001. D Seed germination rate of HZ1005 and Huamai1168. E The method of visual observation and record temporal order to show the difference between the varieties of HZ1005 and Huamai1168 in germination speed Overview of miRNAome data and identification of miRNAs during seed germination To investigate the germination-mediated miRNAs that may be involved in germination growth and PHS tolerance regulations in wheat. We performed miRNAome by using the germinated seeds at different time points (12h, 18h and 36h) for identifying the number, type, and expression difference of miRNAs from RG and SG varieties. Here, total 18 sRNA libraries were constructed and sequenced. The preliminary data obtained from each library are presented in Table S2. Overall, 17,149,204 to 25,727,433 raw reads were obtained from these 18 sRNA libraries. After filtering out the low-quality sequences, repeats, and adapters sequences, 2,343,971 to 5,829,406 clean reads were obtained. Further, the remaining sequences were aligned to mRNA, Rfam and Repbase database to further discard mRNA, other ncRNAs (rRNA, tRNA, snRNA, snoRNA), and repeat sequences, then the unmapped sequences were used for miRNAs identification. Finally, all these reads were classified into ten groups. By comparing the number and percentage of different reads among different libraries, we found that the average ratios of valid reads in RG libraries were higher than those in SG libraries (Fig.S2A-F, Table S2). A total of 1680 pre-miRNAs and 1879 (1862 unique) miRNAs were identified in all samples. All these miRNAs were divided into six groups, then 1331 known miRNAs (401, 150, 131, 607 and 42 in gp1a, gp1b, gp2a, gp2b and gp3, respectively) and 548 novel miRNAs (548 in gp4) were obtained (Table S3). Among the different groups, the number of pre-miRNAs and unique miRNAs in gp2b for each sample was the largest, and the number in gp3 for each sample was the least (Fig. S2G and I, Table S3). All mapped sRNAs within the length range of 18–25 nucleotides (nt) were counted for the total and unique reads, most of samples showed the highest abundance at 24-nt (Fig. S2 J and K). Interestingly, the SG12 samples had higher abundance in 24 nt sequences than that of RG12 samples (Fig. S2K). The length distribution of miRNAs showed that those sequences of 21 nts were the most abundant size class of the known miRNAs followed by 18, 19, 20 and 22 nts, while two most abundant novel miRNAs were 21 and 24 nts (Fig. S2L and M). Analysis of nucleotide bias at each position of miRNAs showed that the first nucleotide tended to be uracil (U) (Fig. [69]2A). However, we found that in the novel miRNAs of 24nt, more miRNAs had adenine (A) at the the first base rather than uracil (U) (Fig. [70]2B). Overall, about 70% of miRNAs were expressed at low levels in each sample (Norm value < 10), and the Norm values of 500–1000 in RG samples were higher than those in SG samples at any time point (Fig. [71]2C). The miRNAs are highly evolutionarily conserved in different species. To reveal the conservation of identified miRNAs with other species in this study, we compared pre-miRNAs with other species in miRbase. Then, 104 miRNAs were found to be highly conservative with their Homologue in soybean, which was the largest number among all tested species, but only 99 homologous miRNAs were found from wheat (Fig. [72]2D). Therefore, most of miRNAs identified in our study were different from those already existed in the miRBase of wheat database. According to the family analysis, 825 known miRNAs were assigned to 85 known families. The TOP 10 miRNA families contained about 55% (450) miRNAs, in which MiR1122 the largest family had the most members of up to 103, followed by the miR156 family and miR5067 family with 58 and 46 members respectively (Fig. [73]2E). In order to find miRNAs distributed in clusters, we co-localized all known and novel miRNAs on wheat genome. Based on the screening standard of inter-distance < = 50000 nts, total 121 pre-miRNA clusters containing 327 miRNAs were identified. Moreover, the distribution of these pre-miRNA clusters was uneven in different chromosomes, and the chromosome 5D had the largest number with 10 clusters, followed by 6A, 6B, 1D and 3D with 8 clusters (Table S4). Interestingly, 7 MiR6300 members from 3 consecutively clusters on 1B and 10 MiR7757 members from 2 consecutively clusters on 3A, which constituted two hot spots of MiR6300 and MiR7757 families respectively (Fig. [74]2F), may play very important roles in seed germination. Fig. 2. [75]Fig. 2 [76]Open in a new tab Comparative analysis of miRNAs in this study. A The distribution of the first base of known miRNAs with the length of 18–25 bp. B The distribution of the first base of novel miRNAs with the length of 18–25 bp. C Distribution of the miRNA expression levels in different samples. D The conservation analysis of miRNAs investigated in this study. E The number of miRNAs in different miRNA families. F Analysis of the miRNA cluster in wheat chromosomes based on the genome-wide location of pre-miRNAs To enhance the reliability and stability of miRNAs identified here, we filtered out the miRNAs expressed only in a library of three biological repeats, then 1103 miRNAs obtained. Afterward, 609, 551 and 562 miRNAs were detected in SG12, SG18 and SG36 respectively, of which 424 miRNAs were co-expressed in the three time points (Fig. S3A). 549, 736 and 821 miRNAs were detected in RG12, RG18 and RG36 respectively, of which 431 miRNAs were overlapped among three time points (Fig. S3B). In addition, 354 miRNAs were co-expressed in all samples, and the 270 miRNAs were special-expressed in RG samples but only 46 miRNAs were special-expressed in SG samples (Fig. S3C). Based on p-value < 0.05, there were 408 differentially expressed miRNAs between RG12 and SG12, of which 206 were up-regulated and 202 were down-regulated. There were 372 differentially expressed miRNAs between RG18 and SG18, of which 119 were up-regulated and 253 were down-regulated. There were 295 differentially expressed miRNAs between RG36 and SG36, of which 96 were up-regulated and 199 were down-regulated (Fig. S3D-G, Table S5). To verify the results of miRNA-Seq data, six miRNAs (zma-MIR164h-p5_2ss4CG17CG_1, tae-miR159a, tae-miR9652-5p, bdi-MIR169c-p3_1, osa-miR444a-3p.2_R + 1, and PC-5p-8258_1835) were selected randomly for the validation of their expression patterns by qRT-PCR. It could be seen that the expression patterns of these miRNAs from qRT-PCR data showed the same trend as those from miRNAome data (Fig. S4). Therefore, the RNA-Seq results were reliable in this study. Comparative miRNA expression profiles by differential-expression analysis and time-series analysis In order to further analyze the roles of miRNAs in seed germination, we used different strategies to search for candidate miRNAs. Firstly, the differentially expressed miRNAs in common among all three comparison groups were excavated from the two varieties of wheat. Based on the overlap analysis, only 94 shared miRNAs were obtained, but 154, 121 and 122 differential miRNAs were specifically detected in SR12, SR18, and SR36, respectively (Fig. [77]3A). Because the 94 miRNAs were differentially expressed in all three time points between SG and RG samples, they were defined as SR-oDE (SG and RG-overlapped differential expression) miRNAs and selected to be candidates for subsequent analysis. Moreover, cluster analysis showed that most of SR-oDE miRNAs were downregulated in SG samples (Fig. S5), which was consistent with the result of Upset analysis, 49 downregulated miRNAs accounted for 52.1% of the total (Fig. [78]3B). Fig. 3. [79]Fig. 3 [80]Open in a new tab Differential expression and expression pattern analysis for exploring the important miRNAs based on miRNAome data. A The overlap analysis of differentially expressed miRNAs among three comparison groups. B The upset analysis of differentially expressed miRNAs among different samples. C The overlap analysis of differentially expressed miRNAs between the comparison groups of SSS and RRR. D Identification of the temporal expression trends for SSS & RRR overlapped miRNAs in SG or RG samples. E Identification of the temporal expression trends for SSS-special & RRR-special miRNAs in SG or RG samples. F The expression patterns of SSS and RRR overlapped miRNAs from Profile 1 and Profile 8 were shown by cluster diagram. G The expression patterns of SSS-special & RRR-special miRNAs from Profile 1 and Profile 8 were shown by cluster diagram Secondly, the stage-differential miRNAs were excavated among different time points of one variety. Based on the time series analysis, 449 and 626 differentially expressed miRNAs were identified from the multi-group comparisons of SG12vsSG18vsSG36 (SSS) and RG12vsRG18vsRG36 (RRR), respectively (Table S6). 364 shared miRNAs were obtained between the SSS and RRR comparative groups, specially 85 miRNAs presented in SSS and 262 presented in RRR (Fig. [81]3C). In order to further hunted for more important candidates, the expression patterns of 364, 85 and 262 miRNAs were performed the trend analysis under the same criterion (fold-change > 1.5) in SG or RG samples. Then, SSS and RRR overlapped miRNAs, SSS-special and RRR-special miRNAs were all divided into eight subclasses of profiles (Fig. [82]3D and E). In particular, the most dramatically changed miRNAs in Profile 1 and Profile 8 were selected in each group. After aggregating all miRNAs in Profile 1 and Profile 8 subclasses, the expression patterns of 91 SSS and RRR overlapped miRNAs were shown by cluster diagram (Fig. [83]3D and F). Obviously, we found that the expression levels of miRNAs in R12 and R18 in cluster II was relatively low but the expression levels of miRNAs in R12 and R18 in cluster IV was relatively high (Fig. [84]3F). The expression patterns of 53 SSS-special and RRR-special miRNAs were also presented by cluster diagram together (Fig. [85]3E and G), then we found that the expression levels of miRNAs in R12 and R18 in cluster I was relatively low but the expression levels of miRNAs in R12 and R18 in cluster II was relatively high (Fig. [86]3G). Subsequently, the miRNAs with high expression levels in R12 and R18 in these two heatmaps were summarized into a group as R12R18-HE miRNAs, and the miRNAs with low expression levels in R12 and R18 in the two heatmaps were summarized into a group as R12R18-LE miRNAs (Fig. [87]3F and G, Table S7). Finally, we got the three most important candidate-miRNA groups SR-oDE miRNAs, R12R18-HE miRNAs, and R12R18-LE miRNAs, contributing to the differences of germinated speed between RG and SG wheat. Degradome analysis Degradome analysis was used to identify mRNAs that were post-transcriptionally silenced by miRNAs. To analyze the miRNA-induced mRNA degradation during seed germination, degradome sequencing was performed for high throughput screening the target genes of identified miRNAs in this study. Detailly, 16,002,713 and 11,780,460 raw reads were detected from RG and SG degradome libraries respectively. After filtering out the reads < 15 nts, 11,748,726 mapped reads were got in RG library, accounting for 73.42% of the raw reads, which were covered for 94474 transcripts in wheat genome, In SG library 8,565,069 mapped reads were obtained, accounting for 72.71% of the raw reads, which were matched to 86847 transcripts (Table S8). The expression levels of transcripts from degradome were analyzed for the preliminary understanding of differential degradation pattern in response to germination between RG and SG libraries. The results showed that 42203 transcripts were differentially expressed, including 17542 up-regulated transcripts and 24661 down-regulated transcripts in the comparison group of RG vs SG (Fig. S6A). We performed GO and KEGG enrichment analysis on these differentially expressed genes. Through GO enrichment, we found that TCA cycle, translation initiation factor activity and glycolysis process were the three most significant enriched pathways (Fig. S6B). Through KEGG classification and enrichment analysis, we found that translation and carbohydrate metabolism were the two pathways with the highest number of genes (Fig. S6C), and the TCA cycle was the most significantly enriched pathway (Fig. S6D). GO and KEGG pathway enrichment analysis for the identification of key miRNAs The miRNA–mRNA modules contributing to the contributing to seed germination were discovered based on the combination of sRNA-seq and degradome-seq data. subsequently, 1636 miRNAs and 22308 mRNAs were constructed into 42203 miRNA–mRNA pairs (Table S9). To explore the specific functions of candidate miRNAs mentioned above, their target genes were screened and obtained from degradome data. Afterwards, three target gene datasets TGs1, TGs2, and TGs3 were constructed for SR-oDE miRNAs, R12R18-HE miRNAs, and R12R18-LE miRNAs, respectively. Then, 1660 miRNA–mRNA pairs were constructed for SR-oDE miRNAs and their targets (Fig. [88]4A, Table S10), 1614 miRNA–mRNA pairs were constructed for R12R18-HE miRNAs and their targets (Fig. [89]4B, Table S10), and 701 miRNA–mRNA pairs were constructed for R12R18-LE miRNAs and their targets (Fig. [90]4C, Table S10). To reveal miRNA-regulated mRNA targets with regulatory functions contributing to seed germination, the target genes from TGs1, TGs2, TGs3 datasets were further analyzed by GO annotation and KEGG enrichment, respectively. In general, functional classification of target genes was performed by GO analysis under biological processes (BP), cell components (CC) and molecular functions (MF) categories. It was shown that in all three target genes datasets, the "metabolic processes", "cells" and "binding" are the most abundant items in BP, CC and MF categories, respectively (Fig. S7A-C). Moreover, the "Protein folding" (GO:0006457) and "double-stranded DNA binding" (GO:0003690) are the two most significantly enriched items in TGs1 (Fig. S7D). The "heat shock protein binding" category (GO:0031072) is the most significantly enriched item in TGs2 (Fig. S7E). "5-phosphoribose 1-diphosphate Biosynthesis process" (GO:0006015), "5-phosphoribose 1-diphosphate metabolic process" (GO:0046391), and "ribose phosphate diphosphokinase complex" (GO:0002189) are the top three significantly enriched items in TGs3 (Fig. S7F). Obviously, the GO analysis results showed the complex and significant differences among different target gene datasets, resulting in hard to identify and check effectively which dataset is more important. Fig. 4. [91]Fig. 4 [92]Open in a new tab Identification and expression analysis of candidate miRNAs and their targets based on miRNA-mRNA network and KEGG enrichment analysis. A The miRNA-mRNA network and KEGG enrichment analysis based on SR-oDE miRNAs and their targets. B The miRNA-mRNA network and KEGG enrichment analysis based on R12R18-HE miRNAs and their targets. C The miRNA-mRNA network and KEGG enrichment analysis based on R12R18-LE miRNAs and their targets. D The expression profiles of candidate miRNAs and their targets involved in starch and sucrose metabolism were verified by qRT-PCR Subsequently, KEGG analysis was used to explore the metabolic pathways regulated by target genes. 94, 105 and 78 pathways were identified and classified for the three target gene datasets TGs1, TGs2, and TGs3 respectively (Table S11). The top 25 of KEGG enrichment analysis are shown in Fig. [93]4 (Fig. [94]4A-C). Specially, it was interesting that the KEGG term “starch and sucrose metabolism” was only enriched in TGs2 dataset (Fig. [95]4B), possibly close association with seed germination. Definitely, a great number of studies have suggested that the starch, hydrolyzed to small oligosaccharides and glucose, which are transported to the embryo to support the growth of the developing seedling, plays a very important role in germinated grain [[96]24]. In order to comprehensively and accurately grasp gene expression information in our study, the mRNA profile was employed to identify the gene expression features during seed germination via high-throughput sequencing. Total 52454 genes were expressed during seed germination in wheat (Table S12). Afterwards, the expression information of the relevant genes involved in starch and sucrose metabolism of TGs2 dataset were screened and filtered from the transcriptome data. Under a |fold-change| > 2 threshold for selecting candidates, it was shown that only three genes TraesCS7B02G034600, TraesCS5A02G265800 and TraesCS1A02G422900 displayed significantly different expression between RG and SG varieties (Table S13). Meanwhile, qRT-PCR were employed to verify the expression levels of miRNAs and target genes. We found that the expression level of the gene (TraesCS7B02G034600) encoding Pullulanase (PULA) in RG36 was three times higher than that in SG36. TraesCS7B02G034600 was the target of tae-miR9662b-p5_1ss9CG whose expression level in RG36 was also nearly three times higher than that in SG36. TraesCS5A02G265800, encoding Beta-glucosidase (BG), the expression level of which in RG12 and RG18 was obviously higher than that in SG12 and SG18, was the target of gma-miR6300_L-1R + 5 whose expression level in RG18 and RG36 was also nearly three times higher than that in SG18 and SG36. TraesCS1A02G422900, encoding Glucan endo-1,3-beta-glucosidase (GEBG), the expression level of which in RG12 and RG18 was obviously higher than that in SG12 and SG18, was the target of ata-miR9674c-3p whose expression level in RG12 and RG18 was also higher than that in SG12 and SG18 (Fig. [97]4D). Constructing starch and sucrose metabolic network for the exploration of miRNA–Target pairs in starch metabolism The starch and sucrose metabolism (SSM), as one of most important pathways in which the degradation of starch and the increase of soluble sugar provide energy for seed germination, has always been the focus in the studies of PHS. In order to systematically and completely explore the starch metabolic pathway in our study, according to KEGG database the related genes involved in starch and sucrose metabolism were scanned and selected from the mRNA-seq data. After that, total 510 expressed genes were obtained (Fig. [98]5A, Table S14). On this basis, the metabolic network from starch to glucose was reconstructed and drown for wheat, including 30 kinds of enzymes encoded by these 510 structural genes which were assigned to 26 metabolic nodes (Fig. [99]5B). Fig. 5. [100]Fig. 5 [101]Open in a new tab the metabolic network from starch to glucose systematically were reconstructed in wheat based on RNA-seq data and the key miRNA-RNA pairs involved in starch and sucrose metabolism were revealed. A Hierarchical cluster analysis of all expressed genes involved in starch and sucrose metabolism from different samples. B The metabolic network from starch to glucose was reconstructed and drown for wheat in seed germination. The expression trends and expression levels of structural genes between RG and SG varieties were shown in different metabolic nodes for representing possible enzyme active intensities. The expression of each gene showed in the histogram was the total expression levels of all transcripts of all homologous genes detected in the samples. AMY: Alpha-amylase; BAM: Beta-amylase; ISA: Isoamylase; GBE: 1,4-alpha-glucan-branching enzyme; GLGSY: Glycogen synthase; SSY: Starch synthase; GBSS: Granule-bound starch synthase; AGP: Alpha-1,4 glucan phosphorylase; GPA: Glucose-1-phosphate adenylyltransferase; PGM: Phosphoglucomutase; UGPA: UTP–glucose-1-phosphate uridylyltransferase; GPG: glucose-1-phosphate guanylyltransferase; SPSA: Sucrose-phosphate synthase; SPP: Sucrose-phosphatase; G6PC: Glucose 6-phosphate phosphatase; CESA1: Cellulose synthase (GDP-forming); CESA2: Cellulose synthase (UDP-forming); SUS: Sucrose synthase; AG: Alpha-glucosidase; DEP: 4-alpha-glucanotransferase; GUN: Endoglucanase; BG: Beta-glucosidase; CALS: Callose synthase; GEBG: Glucan endo-1,3-beta-glucosidase; TPS: alpha-trehalose-phosphate synthase; TPP: trehalose-phosphate phosphatase; TRE: trehalose; INV: invertase; SST: sucrose:sucrose 1-fructosyltransferase; PULA: Pullulanase. C The putative mode for understanding the differentiation of starch degradation and glucose accumulation caused the differentiation of germination and growth process between RG and SG varieties. D Construction of the miRNA-mRNA network involved in starch and sucrose metabolism. E The two largest miRNA families MIR6300 and MIR5054 involved in starch metabolism could regulate AMY1 and AMY2 expression levels. F The interaction of SPP and SUS genes with their relevant miRNAs was identified and classified from the miRNA-mRNA network Here, it was shown that the expression trends and expression levels of most genes were similar in the network between these two varieties, and only a few genes changed in the expression levels were obvious, such as BG, GEBG, PULA identified and mentioned above. In addition, we also found that the expression levels of AMY, SUS, SPP were extremely different among SG and RG samples. The expression level of AMY2 in RG samples was significantly higher than that in SG samples. The expression level of AMY1 in RG18 was more than double that in SG18.The expression level of SUS in RG 18 and RG36 was significantly higher than that in SG 18 and SG36, respectively. The expression level of SPP in RG samples was more than double that in SG samples at each time point (Fig. [102]5B). However, the expression levels of TREH in SG36 was more than twice that found in RG36, but the FPKM expression values of TREH was very low in all samples and not more than 25. Overall, the seven important candidate enzymes, AMY2 and AMY1 associated with starch degradation, SUS and SPP involved in sucrose generation, BG, GEBG and PULA related to D-glucose production, were notably up-regulated in RG samples at almost all the timepoints (Fig. [103]5B). The fold changes of these gene expression levels between different comparisons confirmed this result (Fig. S8A). Due to the gene-dosage effect of high expression of related genes, it can provide more enzymes for catalysis in RG variety. Therefore, our results suggested that starch degradation and glucose accumulation are extremely fast in RG compared to SG, which probably caused the germination and growth process of RG variety faster than SG variety. That is, the shorter time may be required for seed germination in RG variety (Fig. [104]5C). Moreover, we also found that AMY1, AMY2, PULA, and SUS had the greatest expression differences during germination, and the AMY1, BG and GEBG were the predominantly expressed genes in germinated seeds (Fig. S8B). In addition, the metabolic pathway from α-D-Glucose 1-phosphate to GDP-glucose then to cellulose was absent because none of GPG and CESA1 expressed in wheat. Because there were no detectable G6PC expressed during seed germination, the metabolic reaction is virtually absent from D-Glucose 6-phosphate to D-glucose (Fig. [105]5B). In order to fully understand the role of miRNAs in starch metabolism, we scanned the corresponding miRNAs for all relevant genes identified in starch pathway by using the degradation database (Fig. S9), and then widely constructed a miRNA-target network, including 261 miRNAs and 230 target genes formed 574 miRNA-mRNA pairs (Fig. [106]5D, Table S15). Among them, we found that the MIR6300 and MIR5054 were the two largest families involved in starch metabolism, 17 members of MIR6300 and 7 members of MIR5054 regulated 36 and 23 target genes respectively (Fig. [107]5D). Especially, the important miRNA-mRNA pairs, such as AMY1 gene members targeted by bdi-miR5054 and AMY2 gene members targeted by gma-miR6300_1ss18GC, were revealed in this study (Fig. [108]5E). Moreover, the expression levels of bdi-miR5054 were significantly higher in RG samples than that in SG samples, thus it showed significant synergistic regulatory patterns between bdi-miR5054 and AMY1 members. However, the expression levels of miR6300_1ss18GC were higher in SG samples than that in RG samples, thus it showed significant antagonistic regulatory patterns between miR6300_1ss18GC and AMY2 members (Fig. [109]5E). Eventually, the lower AMY1 and AMY2 expression levels may lead to the slower degradation of starch in SG variety. Of all the SPP gene members, only TraesCS5A02G002700 (SPP2) was identified to be targeted by hbr-MIR166a-p5_2ss6AG17GA and osa-MIR1863b-p3_1ss19AT, but they never showed the significantly differences in expression patterns between the different varieties. 5, 2, 3 and 1 genes encoding SUS1, SUS2, SUS4 and SUS7 respectively were found in wheat from this study. Especially, TraesCS7A02G158900, TraesCS7B02G063400 and TraesCS7D02G159800, which are the three of five SUS1 members, were dominant genes with the extremely high expression in SUS family during seed germination. Meanwhile, the expression levels of them were significantly higher in RG samples than that in SG samples at every timepoint. Interestingly, these three genes could be targeted by five miRNAs in which only tae-MIR9672a-p5 and PC-5p-1563_8312 showed higher expression levels and greater expression differences in different samples (Fig. [110]5F). Therefore, three miRNA-mRNA pairs, PC-5p-1563_8312–TraesCS7A02G158900, PC-5p-1563_8312–TraesCS7D02G159800, and tae-MIR9672a-p5–TraesCS7B02G063400, have been constructed and their antagonistic regulatory patterns were revealed. It indicated that the higher expression levels of tae-MIR9672a-p5 and PC-5p-1563_8312 contributed to the lower expression levels of SUS1 in SG samples (Fig. S10). The other two SUS4 genes with low expression levels had no significant differences between SG and RG samples (Fig. [111]5F), although which could be targeted by relevant miRNAs. Therefore, lower SUS1 expression levels may eventually lead to a slower rate of sucrose production in SG variety. Regulatory network of candidate miRNAs and target genes correlated with seed germination in PHS transition Starch is broken down to sugars which could supply the developing embryo. Based on the comprehensive analysis described above, we proposed a possible model which control starch and sucrose metabolism pathway involved in PHS transition by miRNAs. It was found that the difference in starch degradation rate in seeds leads to difference in seed germination speed. Then, we discovered and constructed the miRNA-mRNA regulatory network based on the 7 key miRNAs and their targets associated with starch and sucrose metabolism pathway. In the regulatory network, two major α-amylase isozyme families AMY1 and AMY2 were regulated by bdi-miR5054 and miR6300_1ss18GC, respectively. PULA was the target gene of tae-miR9662b-p5_1ss9CG. BG was the target gene of gma-miR6300_L-1R + 5. GEBG was the target gene of ata-miR9674c-3p. In addition, SUS1 members were regulated by tae-MIR9672a-p5 and PC-5p-1563_8312. In the pathway, the α-amylase catalyes the breakdown of starch for producing maltose and dextrin, then dextrin is decomposed and conversed into glucose by pullulanase. The content of sucrose was affected by sucrose synthetase. The glucan endo-1,3-beta-glucosidase converts 1,3-β-D-Glucan to D-Glucose and the beta-glucosidase converts cellodextrin to D-Glucose. Due to the roles of miRNAs in the regulation of these genes encoding the enzymes involved in starch and sucrose metabolism, the starch degradation and sugars production were much more slowly in SG variety compared with that in RG variety during seed germination. Eventually, when encountering disaster weather of continuous rain in harvesting stage, the germination speed and seedling growth of SG variety may be slower than that of RG variety, so SG variety showed PHS resistance (Fig. [112]6). Fig. 6. Fig. 6 [113]Open in a new tab Putative miRNA regulatory network of PHS transition in wheat. Seven key miRNAs regulate the accumulative variation of starch, maltose, dextrin, sucrose and glucose, which were the important components in starch metabolism, to achieve the regulation of PHS. The miR6300^1 represents miR6300_1ss18GC, and miR6300^2 represents gma-miR6300_L-1R + 5 Discussion PHS is actually a bidirectional regulation of seed dormancy and germination by the environment [[114]25]. PHS is closely associated with the degree of seed dormancy, which inhibits seed germination under optimal environmental conditions [[115]26]. In fact, the relationship between seed dormancy and seed germination is inseparable in the study of PHS. Identification of pre-harvest sprouting resistance mainly include whole-spike germination test and grain germination test [[116]27]. One other study has demonstrated an increase in seed germination with the extension of post-ripening time, indicating a gradual release of seed dormancy [[117]28]. Previous studies have revealed the roles in regulating seed germination in model plants such as Arabidopsis thaliana, rice (Oryza sativa), maize and other plants (Zea mays). After deep sequencing of small RNA, 590 differentially expressed miRNAs were uncovered from the early seed germination mutant (eno2^−) of Arabidopsis thaliana [[118]15]. A total of 247 miRNAs (104 known and 143 novel) were identified as the key factors for responding to ABA and GA in maize embryos during seed germination [[119]29]. A total of 68 differentially expressed miRNAs were identified and used to construct the interaction network between miRNAs and their target genes, which resulted in 571 miRNA–mRNA pairs during Magnolia sieboldii seed germination [[120]30]. A total of 145 known miRNAs belonging to 47 families and 78 novel miRNAs were identified from the germinated seeds of ancient eudicot Nelumbo nucifera by using small RNA sequencing [[121]31]. 88 conserved miRNAs representing 25 defined families and 13 novel miRNAs were specifically implicated in the germination of salt-stressed R. soongorica seeds [[122]32]. Multi-omics was used to reveal the molecular mechanism of seed germination in maize, and then 437 known and 354 novel miRNAs with different expression levels were identified for regulating germination induced by gibberellin (GA) [[123]33]. Through the analysis of small RNA sequencing to reveal the molecular mechanism of seed germination in Mung Beana, there were 284 miRNAs found in the imbibed and dried seeds of mung bean containing 213 known and 71 novel miRNAs [[124]16]. Here, we performed our investigation on miRNAs-involved regulation for seed germination by a systematical multi-omics analysis in wheat. Through small RNA sequencing analysis, a total of 1879 (1862 unique) miRNAs distributed on all chromosomes of wheat were identified in this study, of which 1331 were known miRNAs and 548 were novel miRNA (Table S2, Fig. S1 and Fig. [125]1F). After combined miRNAome data with degradome data, 1636 miRNAs and 22308 mRNAs were constructed into 42203 miRNA–mRNA pairs (Table S8). Therefore, this work provides sufficient data for further studying on the role of miRNAs in wheat germination. Of course, the availability of so many miRNAs is also closely related to the large genome size of wheat. It also indicated that the molecular mechanism involved in wheat seed germination may be more complex than that of other plants. Through a comprehensive analysis, it was found that the miR167–ARF–sugar metabolism genes (GmHXK, GmPFK, GmTPS, and GmFRK) pathway facilitated carbon allocation after starch degradation, potentially contributing to low seed starch content in soybean [[126]34]. The sugarcane miRNAs, miR172, miR164, miR396 and miR169, regulated AP2/ERF, NAC, GRF and bZIP TF members associated with sugar metabolism to especially control of sucrose accumulation in sugarcane [[127]35]. Zma-miR159i-3p and Zma-miR159k-3p were predicted to post-transcriptionally regulate the expression of ZmMYB138 and ZmMYB115 which directly influenced the activity of SSG promoters, indicating that miRNA–TF–SSG regulatory network has an essential role in maize starch biosynthesis [[128]36]. According to the analysis of transcriptome data and small RNA-seq data, it is showed that NnumiR396a-NnSS2 could prevent the synthesis of amylopectin and NnumiR396b-NnPGM2 could affect the synthesis of total starch in Lotus [[129]37]. However, few studies have systematically and comprehensively depicted miRNA-mediated sugar metabolism pathway. Here, by integrating the miRNAome data with the transcriptome data and degradome data, we found that miRNAs play important roles in starch and sucrose metabolism (Fig. [130]4B). Subsequently, 510 genes involved in starch and sucrose metabolism were identified through genome-wide scanning, and then a regulatory network of miRNAs-targets containing 574 miRNA-mRNA pairs was constructed for studying on the regulation of the biological process from starch to glucose (Fig. [131]5A-D, Table S13). Therefore, our study provided the favorable basis for further exploring miRNA-mediated molecular mechanism involved in the starch metabolism. The starch and sucrose metabolism, as one of most important pathways in which the degradation of starch and the increase of soluble sugar provide energy for seed germination, has always been the focus in the studies of PHS. During the transition of seeds from a dormant state to germination, carbohydrate metabolism and plant hormone signal transduction pathways are activated [[132]38]. Dormant buds require carbon sources to promote growth, and the plant adjusts sugar metabolism between storage and soluble sugars to maintain dormancy and sprouting [[133]39]. Though the systematical determination of sucrose, starch, and glucose levels, it was found that sucrose synthase plays an important role in regulating carbohydrate metabolism during seed dormancy release in P. calleryana [[134]40]. Compared with dark, light could promote the germination of celery seed by inducing the degradation of starch and increasing the content of soluble sugar [[135]41]. In wheat, the full suite of amylolytic enzymes including TaAMY1 and TaAMY2, and proteolytic enzymes are produced to mobilize both C and N reserves for the germinating embryo during PHS [[136]42]. TaAMY1 and TaAMY2 appeared to share common specific structural features with the barley HvAMY2 and HvAMY1, respectively, but only TaAMY2 have indications of the presence of an extra SBS sugar tongue site suggesting a different level of efficiency or a different target [[137]43]. Abnormal accumulation of soluble sugar (α-gluco-oligosaccharide and sucrose) by TaAMY1 over-expression reduced the grain dormancy and enhanced abscisic acid (ABA) resistance. However, the impact of the overexpressed TaAMY1 on germination was very limited [[138]44]. The overexpression of TaAMY2 caused ABA insensitivity and resulted in grains with almost total absence of dormancy before any after-ripening, but the role of α-amylase in starch degradation during germination was really confused, because the TaAMY2 is not required to break down starch in the first stage of germination [[139]45]. Germination studies have shown that HvAMY1 is accumulated during the first 24 h of imbibition, followed by HvAMY4 from 24–72 h of imbibition [[140]46]. Similar to previous studies in barley, we found that TaAMY1 were expressed during the first 12 h of imbibition, followed by TaAMY2 from 18–36 h of imbibition in wheat. Additionally, our results showed that the expression levels of AMY1 and AMY2 were significantly higher in RG samples than that in SG samples (Fig. [141]5B), which may lead to the slower degradation of starch in SG than that in RG. However, the role and impact of TaAMY1 on the early stage of germination remains unknown. The previous study also suggested that sucrose synthase plays an important role in regulating carbohydrate metabolism during seed dormancy release in P. calleryana [[142]40]. In our study, 5, 2, 3 and 1 structure genes encoding SUS1, SUS2, SUS4 and SUS7 respectively were found in wheat from this study. Especially, three SUS1 members TraesCS7A02G158900, TraesCS7B02G063400 and TraesCS7D02G159800 were dominant genes with the extremely high expression in SUS family during seed germination (Fig. [143]5F). Meanwhile, the expression levels of them were significantly higher in RG samples than that in SG samples at every timepoint (Fig. S9). Glucan endo-1,3-beta-glucosidases have been implicated in cell division, pollen development, regulation of plasmodesmata signaling, cold response, seed germination, and maturation [[144]47]. In rice, The GWAS analysis of seed germination indicated that the candidate gene LOC_Os07g07340 (glucan endo-1,3-beta-glucosidase) affecting root dry mass was associated with germination rate [[145]48]. Beta-glucosidase is an enzyme involved in the degradation of cellulose and plays an essential part in many biological processes (Purification and characterization of a recombinant beta-glucosidase in Escherichia coli). In rice, beta-glucosidases displayed dramatic expression levels at 24 h compared with that of 6 h during early seed germination phase [[146]49]. However, there are few reports about glucan endo-1,3-beta-glucosidase and beta-glucosidase contributing to PHS in previous studies. Here, it was implied that the activities of these two enzymes were higher in the rapidly germinated variety (RG) than that of slowly germinated variety (SG) (Fig. [147]5B). Therefore, maybe they were relevant to the PHS phenotype of RG seeds in field under excessive rainfall conditions. Under the influence from the high expressions of AMY1, AMY2, SUS1, PULA, BG, and GEBG enzymes, it is quite possible that the increase of soluble sugar have been accelerated in RG variety during different germination stages, which leads to the higher germination speed of RG than that of SG (Fig. [148]5C). Presumably, it was also relevant to the PHS phenotype of RG variety. Our results also confirmed that the starch degradation plays a very important role in PHS resistance (Fig. [149]6). The previous studies have shown that miRNAs were involved in PHS regulation of wheat. A wheat-specific microRNA (miR9678) could improve the PHS resistance of wheat by targeting a long noncoding RNA (WSGAR) and triggering the generation of phased small interfering RNAs which play a role in the delay of seed germination [[150]21]. The wheat miR27319 played the positive role in regulation of seed dormancy through phytohormone abscisic acid and GA biosynthesis, catabolism, and signaling pathways [[151]50]. The integrative analysis of transcriptome and microRNAs revealed the underlying transcriptional co-regulation mechanism of miR1120, miR1130, miR1133 and miR408 for controlling GA biosynthesis during seed germination of wheat [[152]20]. Here, we also found several new miRNAs which could play important roles for PHS resistance in wheat though mediating the starch metabolism pathway initiated during seed germination (Fig. [153]6). Conclusions Although PHS transition in wheat is a highly complex process, it is certain that miRNAs must play very important roles in the regulation of seed germination and dormancy. By integrating the miRNAome data with the transcriptome data and degradome data, the miRNA–mRNA pairs who were related to seed germination speed between two white wheat varieties were screened. Through mining and analyzing inherent mechanism, the network effectively integrating by miRNA–mRNA regulatory models of bdi-miR5054–AMY1, miR6300_1ss18GC–AMY2, tae-miR9672a-p5–SUS1, PC-5p-1563_8312–SUS1, tae-miR9662b-p5_1ss9CG–PULA, gma-miR6300_L-1R + 5–BG, and ata-miR9674c-3p–GEBG during the process ofstarch metabolism were supposed to be the essential regulation center for seed germination of wheat, which also could play a critical role on the PHS resistance. Hence, this study not only provides more clues for understanding the control of gene expression during seed germination, but also reveals the miRNA-mediated regulation mechanism for forming intrinsic and inherent differences, resulting in the significant difference of PHS phenotype between the white wheat varieties. Overall, our study will provide an important theoretical basis for pre-harvest sprouting resistance improvement in wheat breeding. Supplementary Information [154]Additional file 1.^ (2.6MB, docx) [155]Additional file 2.^ (15.1MB, zip) Acknowledgements