The wheat landrace WTB with high dormancy and PHS resistance was subjected to HT and NT treatments, as well as whole-transcriptome sequencing. The WTB plants were potted in a natural field environment (31°58′N, 117°240′E, Hefei, China, Supplemental Table 1 ) during 2017–2018 cropping seasons before anthesis. All fields were maintained as disease- and weed-free. After selecting the representative 60 main stem spikes to mark in anthesis, the plants were transferred to two illumination incubators for HT (35°C day/25°C night) or NT (25°C day/20°C night) treatment (16 h day/8 h night, 70% relative humidity) from 21 DPA to 35 DPA. Seeds were collected at 21 DPA, 28 DPA, and 35 DPA, named 21DPA-NT, 28DPA-NT, 28DPA-HT, 35DPA-NT, and 35DPA-HT; immediately frozen in liquid nitrogen; and stored at −80°C for RNA extraction.

The TaGA20ox1 mutant (designated ga20ox1, jm_chr3D_508613303) and the TaCDPK21 mutant (designated cdpk21, jm_chr7A_270043204) of wheat cultivar JM22, created using EMS provided by Yantai Jien Biotechnology Company (http://jm.ytjebc.com/), were used to determine the seed germination percentage (GP). The cdpk21 mutant plants were grown in growth chambers (16 h day/8 h night, 70% relative humidity). The ga20ox1 mutant plants were also grown in growth chambers (16 h day/8 h night, 70% relative humidity) and treated at HT (NT as control) from 21 DPA to 35 DPA, which is consistent with the treatment methods used for the WTB plants.

All Arabidopsis plants used had the Arabidopsis thaliana Columbia (Col-0) accession background. The atcdpk24 mutant (SALK_015986C) was ordered from The Arabidopsis Information Resource (TAIR, https://www.Arabidopsis.org/). Homozygous mutant plants were identified in the F₂ generation through PCR by using primers obtained from the Salk Institute Genomic Analysis Laboratory (SIGnAL) database ( Supplemental Table 2 ). For planting, seeds were imbibed at 4°C for 3 days for stratification and were harvested from plants cultivated in soil in a greenhouse at 21°C ± 1°C (16 h day/8 h night, 60% relative humidity).

Rice (Oryza sativa L. japonica Nipponbare) plants were grown in a walk-in growth chamber under a 12-h light (28°C)/12-h dark (22°C) photoperiod regime. Homozygous lines (T₂ generation) were selected for the subsequent experiments.

GP was used to evaluate SD. Wheat seed GP was determined according to the method of Zuo et al. (2019). Germinated seeds were selected on the basis of seed coat rupture at the embryo site, and GP was calculated as the number of germinated seeds by the 7th day of the germination test divided by the total number of seeds.

For Arabidopsis seeds, germination phenotypes were determined using the method of Cao et al. (2020). Briefly, approximately 50 seeds of individual plants were sown on filter paper saturated with sterile distilled water in Petri dishes and cultivated in a growth chamber (16 h light/8 h dark, 21°C ± 1°C). The seeds were collected at 25 DPA, and the GP was measured on the 7th day of the germination assay. Seeds with a protruded radicle were considered to be germinated.

For rice seeds, germination phenotypes were determined using the method of Li et al. (2020). Rice seeds with a ≥ 1-mm-long radicle were considered to be germinated.

To determine ABA, GA, and α-amylase activities, freshly harvested rice seeds were imbibed in the growth chamber (16 h light/8 h dark, 21°C ± 1°C) for 24 h before sampling. Approximately 100 mg samples of the imbibed seeds were used. ABA and GA levels were quantified using the ELISA Kit (Meimian, BY-JZF0146 and BY-JZF0147) according to the manufacturer’s instructions. Alpha-amylase activity was measured using a kit (Solarbio, BC0615) according to the manufacturer’s instructions.

Total RNA was isolated from 15 samples (21DPA-NT, 28DPA-NT, 28DPA-HT, 35DPA-NT, and 35DPA-HT, 3 repeats for each sample) by using the MiniBEST Plant RNA Extraction Kit (TaKaRa) in the same manner as for transcriptome sequencing. The complementary DNA (cDNA) library was constructed using the Ribo-Zero Magnetic Kit (Madison, USA). The procedure was performed in accordance with that described in a study by Zhang et al. (2016). Libraries for small RNA sequencing were constructed using the NEB NextR UltraTM small RNA sample library prep kit (NEB, USA). Strand-specific sequencing for mRNAs and small RNAs was performed on an Illumina HiSeq4000 system by Biomarker technologies (Beijing, China); sequencing was performed using the standard Illumina protocol according to Ma et al. (2015).

Reference genome sequences and annotations were downloaded from the wheat IWGSC v1.1 genome website (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.1/). Reads containing the adaptor sequence, duplicate sequence, bases with a quality score of < Q10 more than 50%, and with uncertain bases (reads with > 10% unknown sequences “N”) were removed. For the RNA-seq data, the clean reads were filtered from the raw reads mapped to the wheat reference genome by using Bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2) (Langmead and Salzberg, 2012) and HISAT2 (http://www.ccb.jhu.edu/software/hisat) (Kim et al., 2015) with default parameters. The Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) method was used to measure transcript expression levels by using Cufflinks software (http://cole-trapnell-lab.github.io/cufflinks) (Trapnell et al., 2012). Reads uniquely mapped to the reference sequences (with ≤ 1 mismatch) and FPKM ≥ 1 in at least one sample were used for identifying differentially expressed genes (DEGs). DEG analysis was performed using DEGseq (https://bioconductor.org/biocLite.R) with criteria of fold change ≥ 2 and P < 0.01 (Wang et al., 2010). Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to functionally categorize DEGs by using the KEGG (http://www.genome.jp/kegg/genes.html) databases, respectively (Kanehisa and Goto 2020). TBtools software (https://github.com/CJ-Chen/TBtools) (Chen CJ et al., 2020) was used to identify the overlapping differentially expressed RNAs in different samples and tissues. Transcription factor (TF)-related sequence information was obtained from the Plant Transcription Factor Database (http://planttfdb.gao-lab.org) (Jin et al., 2015).

Following the assembly, cuffcompare (a tool of Cufflinks, http://cole-trapnell-lab.github.io/cufflinks) was used to investigate the location between transcripts and known mRNAs and lncRNAs (Trapnell et al., 2012). The non-coding transcripts were identified using four software, including Coding Potential Calculator2 (http://cpc2.gao-lab.org/) (Kang et al., 2017), Coding-Non-Coding Index (https://github.com/www-bioinfo-org/CNCI) based on the coding potentials that scored < 0 (Sun et al., 2013), Coding Potential Assessment Tool (CPAT, http://lilab.research.bcm.edu) (Wang et al., 2013) based on the length, and coverage of the Pfam protein database (http://pfam.xfam.org/) (Mistry et al., 2020) based on no transcripts was aligned. The lncRNAs were identified when they satisfied all methods and were classified into several categories according to their genomic location and previous description of the Rfam by Inference of RNA Alignments (INFERNAL, http://eddylab.org/infernal/) (Nawrocki and Eddy, 2013). Coding genes located within 100-Kb upstream and downstream of lncRNAs were identified as cis-targets.

To obtain a set of predicted circRNA loci from the RNA-seq data, Find_circ v1.2 (https://github.com/marvin-jens/find_circ) (Memczak et al., 2013) with output filtered for a minimum of 6 reads was used. The raw counts were first normalized using transcripts per million (TPM). CircRNAs with fold changes ≥ 2 and P ≤ 0.05 were classified as circRNAs with significant differential expression.

Following the trimming of adaptor sequences in the 5’ and 3’ ends of sequenced reads, we filtered out low-quality reads (reads with unknown sequences “N”), adaptor sequence fragments, and empty reads. By using the SOAP2 (http://soap.genomics.org.cn) (Li et al., 2009) with parameter settings for perfect matches, reads were aligned to the genome. Then, the unique sequences of length 18–25 nucleotides (nt) were blasted to rice precursors in miRBase v22 (https://www.mirbase.org/) (Kozomara et al., 2019) to detect known miRNAs. Mireap v0.1 (https://help.rc.ufl.edu/doc/Mireap) (Li et al., 2011) was used for predicting novel miRNAs. TPM was used to normalize the counts, and miRNAs with fold changes ≥ 2 and P ≤ 0.05 were categorized as miRNAs with significant differential expression. Target genes of differentially expressed miRNAs were predicted using psRNA Target (http://plantgrn.noble.org/psRNATarget) (Dai and Zhao, 2011).

A set of selected DEGs from differentially expressed transcripts were validated through qRT-PCR by using the same RNA samples as those used for transcriptome sequencing. For transgenic materials, seeds imbibed for 24 h were used for quantitative analysis. The target gene sequences were used to design gene-specific primer pairs by using Primer Premier v5.0 ( Supplemental Table 2 ). The SYBR green Premix Ex TaqTM II (TaKaRa, RR081A) quantitative PCR system was used for qRT-PCR analysis. Approximately 0.5 μg of isolated RNA was used for first-strand cDNA synthesis by using the Prime Script RT reagent kit (Takara, RR037A) following the manufacturer’s instructions for mRNAs and lncRNAs. For cirRNAs, RNase R (Jisai, R0301) was used to remove genomic linear DNA before reverse transcription. Furthermore, instead of oligo (dT) primers, random primers were added to synthesize first-strand cDNA for circRNAs. The method of poly (A) tail was used for miRNA quantitative confirmation, and the miRNA First-Strand Synthesis Kit (Takara, No.638315 and No.638316) was used based on the manufacturer’s instructions. The following amplification program was used: denaturation at 95°C for 2 min and 40 cycles of amplification (95°C for 10 s, followed by 59°C for 30 s). Melting curve analysis was performed from 65°C to 95°C, with 5 s increments of 0.5°C. The wheat Actin gene was used as an internal control for mRNAs, lncRNAs, and circRNAs. The wheat U6 gene was used as an internal control for miRNAs. All samples were analyzed in triplicate. The 2^-△△Ct (Livak and Schmittgen, 2001) relative quantification method was used to calculate the relative expression level.

The TaCDPK21-7A (TraesCS7A02G267000) cDNA sequence from ‘Chinese Spring’ was amplified and cloned into the EcoR I and BamH I sites of the binary vector pCAMBIA1301a under the control of the 35S promoter. The vector was transformed into the Agrobacterium tumefaciens strain CV3101. Subsequently, it was transformed into Arabidopsis Col-0 and atcdpk24 plants by using previously described methods (Clough and Bent, 1998).

The miR27319 full-length sequence (approximately 300 nt before and after the miRNA site) from ‘Chinese Spring’ was amplified and cloned into the Eco31 I sites of the binary vector pCAMBIA1300 under the control of the 35S promoter. To create a vector of short tandem target mimic (STTM), a 48-nt-specific sequence of miR27319 (with both ends including miRNA target binding sites and Eco31 I sites) was recombined into the expression vector pCAMBIA1300. These vectors were transformed into the A. tumefaciens strain CV3101 and subsequently transformed into mature rice seed-derived calli (Nipponbare, Nip) (Hiei et al., 1994) and into the Arabidopsis Col-0 plants (Clough and Bent, 1998) by using previously described methods. The transgenic plants were selected using hygromycin, and transgene insertion was validated through DNA extraction and PCR ( Supplemental Table 2 ).

Clean Illumina sequences were deposited in the National Center for Biotechnology Information Databank (NCBI) (accession number: PRJNA895954).

SPSS v22.0 (http://www.spss.com) software was used for analyzing phenotypic data, and correlation coefficients were estimated using the Pearson method. The Mann–Whitney U-test was performed through stepwise, step-down comparisons.

The GP values of the seeds collected from 28DPA-NT, 35DPA-NT, 28DPA-HT, and 35DPA-HT were 0%, 0%, 9%, and 63%, respectively ( Supplemental Table 3 ), indicating that HT obviously reduces the SD level and promotes germination compared with NT, which is consistent with the previous study results.

As shown in Supplemental Table 4 , 2070.31 million reads in chain-specific libraries were obtained. Of them, 2063.48 million reads were clean reads, accounting for 99.67% of the total reads. Average reads of 85.18% were aligned to the wheat reference genome ( Supplemental Table 4 ). Pearson correlation analysis showed a good relationship for biological replicates (0.952-0.993**, Supplemental Table 5 ). In total, 130312 mRNAs were obtained, 41316 of which were identified as new genes ( Supplemental Table 6 ). Additionally, we obtained 40970 lncRNAs and 4414 circRNAs ( Supplemental Tables 6 , 7 ). Overall, most lncRNAs were approximately 300-500 nt long ( Figure S1A ), and the expression levels of lncRNAs were generally lower than those of mRNAs ( Figure S1B ). The percentages of lncRNA classes suggested that most lncRNAs were long intergenic-lncRNAs (lincRNAs) (86.5%), followed by antisense-lncRNAs (8.2%), intronic-lncRNAs (2.7%), and sense-lncRNAs (2.6%) ( Figure S1C ). CircRNAs were mainly produced from genes and most of them were approximately 200-400 nt long ( Figure S2A ). In total, 2439 circRNAs were derived from the gene ID in IWGSC v1.1, and 1166 circRNAs were derived from new genes predicted in this study ( Figure S2B ).

In total, 344.62 million reads were obtained from miRNA libraries, and 191.41 million clean reads were obtained ( Supplemental Table 4 ). Because primary miRNA transcripts contain one or more miRNA stem loops that could be processed into mature miRNAs, 785 miRNAs were identified from these libraries, including 79 known and 706 novel miRNAs ( Supplemental Table 6 ). The majority of the small RNAs were 21- and 24 nt long ( Figure S2C ). The proportions of first base distribution are shown in Figure S2D ; the proportion of the U base (47%) was higher than those of A (33%), C (14%), and G (6%).

Overall, 3605 and 2372 mRNAs were differentially expressed in 28DPA-NTvsHT and 35DPA-NTvsHT, respectively ( Figure 1A ; Supplemental Table 8 ). Of them, 373 down-regulated mRNAs and 380 up-regulated mRNAs were common to the two groups ( Figure 1A ). In 28DPA-NTvsHT, DEGs were involved in benzoxazinoid biosynthesis; valine, leucine, and isoleucine biosynthesis; starch and sucrose metabolism; and carbon fixation in photosynthetic organisms, as determined through KEGG enrichment analysis ( Figure 1B ). In 35DPA-NTvsHT, DEGs were involved in protein processing in the endoplasmic reticulum and in amino sugar and nucleotide sugar metabolism ( Figure 1C ).

According to the Venn diagram, 94 and 66 miRNAs were differentially expressed in 28DPA-NTvsHT and 35DPA-NTvsHT, respectively. Of them, 14 down-regulated miRNAs and 10 up-regulated miRNAs were common to the two groups ( Figure S3A ; Supplemental Table 9 ). In total, 174 and 137 lncRNAs were differentially expressed in 28DPA-NTvsHT and 35DPA-NTvsHT, respectively. Of them, 28 down-regulated lncRNAs and 3 up-regulated lncRNAs were common to the two groups ( Figure S3B ; Supplemental Table 10 ). In addition, 13 and 11 circRNAs were differentially expressed in 28DPA-NTvsHT and 35DPA-NTvsHT, respectively. Only 3 up-regulated circRNAs were common to the two groups ( Figure S3C ; Supplemental Table 11 ).

To determine the functional relationships between ncRNAs and the target genes involved in HT-regulated dormancy, we constructed a regulatory network based on differentially expressed mRNAs, miRNAs, and lncRNAs ( Figure S4 ; Supplemental Table 12 ). The network map revealed that the relationships between coding and non-coding RNAs were one-to-one, one-to-many, or many-to-one; otherwise, lncRNAs could act as a common target for miRNAs and mRNAs. Among the 849 differentially expressed mRNAs identified in both 28DPA-NTvsHT and 35DPA-NTvsHT, 374 contained miRNA identical regions, which implied that these mRNAs may contain miRNA-binding sites and act as targets for miRNAs ( Supplemental Table 8 ). In the regulatory network, unconservative_chr5A_27319 (designated miR27319) was particularly located at the key node, suggesting its key role in HT-mediated dormancy; however, this result needs to be further validated experimentally ( Figure S4 ; Supplemental Table 12 ). Moreover, 6 lncRNAs were found to be co-expressed with 7 mRNAs, including 5 involved in cis-regulatory mechanisms and 2 trans-regulatory mechanisms ( Supplemental Table 12 ).

To validate transcriptome sequencing results, qRT-PCR was performed to determine the expression patterns of randomly selected 64 mRNAs, 11 miRNAs, 50 lncRNAs, and 5 circRNAs from the 35DPA-NT and 35DPA-HT samples. As shown in Figures S5A, B , and Supplemental Table 13 , the expression of all these selected mRNAs exhibited the same change trends determined through RNA-seq analysis, indicating that the transcriptome sequencing results were reliable. The aforementioned selected candidates were involved in various biological pathways, including 25 genes encoding TFs (such as C3H, MYB, bHLH, MADS, WRKY, GATA, and B-box), 17 genes involved in plant hormone pathways (such as ABA, GA, ethylene, and auxin pathways), 4 peroxidase genes, 2 heat shock protein family genes, 3 heat stress TF genes, 5 cytochrome P450 family genes, and 8 genes involved in other pathways ( Supplemental Table 13 ).

The expression patterns of 9 of the 10 miRNAs were consistent with the sequencing data ( Figure S5C ; Supplemental Table 13 ). According to both the sequencing data and qRT-PCR results, the expression of the known tae-miR408 and tae-miR395b was downregulated in 35DPA-NTvsHT, whereas that of the three reported miRNAs (tae-miR9779, tae-miR9772, and tae-miR9774) was upregulated, indicating that these miRNAs are responsive to HT ( Figure S5C ; Supplemental Table13 ). Additionally, 4 novel miRNAs (unconservative_chr1A_1868, unconservative_chr5A_27319, unconservative_chr5A_26651, and unconservative_chr3A_15287) were validated to be differentially expressed after HT treatment ( Supplemental Table 13 ).

The expression patterns of 50 lncRNAs with a high fold change in 35DPA-NTvsHT were consistent with those obtained from the sequencing data ( Figure S5D and Supplemental Table 13 ). In particular, five lncRNAs (MSTRG.1841654.4, MSTRG.3326929.1, MSTRG.3111190.1, MSTRG.918018.1, and MSTRG.364559.2) exhibited more than 3-fold differences between 35DPA-NT and 35DPA-HT in both RNA-seq analysis and qRT-PCR, implying their critical roles in HT-mediated dormancy ( Supplemental Table 13 ).

The primers for circRNAs were designed for qRT-PCR validation. Five selected circRNAs (chr2D:624323990|624324666, chr2D:624324010|624324666, chr3B:127321663|127322231, chr3D:866248|867183, and chr3D:477885083|477889020) were observed to be processed through back-splicing, confirming the reliability of high-throughput transcriptome sequencing data ( Figure S5E ; Supplemental Table 13 ).

According to Tai et al. (2021), 81 SD-associated genes have been summarized in wheat, with eight of them being differentially expressed in WTB seeds after HT treatment. Expression of the four known dormancy genes (TaDOG1-3A, TaDOG1-3D, TaMFT-3A, and TaMFT-3D) was downregulated in the WTB seeds after HT treatment, which was further validated through qRT-PCR ( Supplemental Table 13 ). The other four known genes (TaSnRK2-2D, TaAFP-2D, GID1-1A, and GID1-1D), which are involved in the ABA and GA signaling pathways, were also responsive to HT in the WTB seeds. This suggested that ABA and GA signals may participate in HT-mediated dormancy ( Supplemental Table 14 ).

Notably, TaDOG1-3A, GID1-1D, TaMFT-3A, and TaMFT-3D were predicted to be targeted by the novel lncRNAs MSTRG.4103020.1, MSTRG.861889.1, MSTRG.729470.1, and MSTRG.891662.2, respectively. This will contribute to broaden our understanding of the molecular mechanisms of TaDOG1 and TaMFT in SD regulation ( Supplemental Table 14 ).

HT significantly upregulated the expression of gibberellin 20 oxidase 1-B-like gene (designated TaGA20ox1, TraesCS3D02G393900) (4.87-fold in 28DPA-NTvsHT and 2.84-fold in 35DPA-NTvsHT), and its expression pattern was verified through qRT-PCR ( Figure S6A ). To validate the role of TaGA20ox1 in HT-mediated SD, we obtained the EMS mutant ga20ox1 with a SNP mutation (G/A) at the 128th amino acid in the functional domain (a non-heme dioxygenase N-terminal domain) of TaGA20ox1. This resulted in the conversion from glutamate (Glu) to lysine (Lys) in the wheat cultivar JM22 ( Figure 2A ). To avoid the interference of the genetic background, the ga20ox1 mutant was backcrossed with JM22 for two generations. In the BC₂F₂ (ga20ox1) mutant seeds, the relative expression of GA20ox1 was significantly decreased in mutant lines compared with JM22 ( Figure 2B ). And the endogenous GA content of TaGA20ox1 was significantly decreased in ga20ox1, revealing that the TaGA20ox1 mutation affects endogenous GA biosynthesis ( Figure 2C ). The ga20ox1 seeds exhibited obviously lower GPs (average GP 27%) than JM22 seeds (average GP 68%) ( Figures 2D, E ). Moreover, compared with JM22 seeds (average GP 91%), ga20ox1 seeds were less sensitive to HT (average GP 37%) ( Figures 2D, E ). These results confirm that TaGA20ox1, as a crucial component in the GA biosynthesis pathway, participates in HT-mediated dormancy.

The calcium signaling pathway is known to be related to SD and seed germination (Zhu et al., 2007; Zhao et al., 2011). Being a calcium sensor, the expression of the gene calcium-dependent protein kinase 21-like (designated TaCDPK21, TraesCS7A02G267000) was significantly upregulated (1.69-fold) in 35DPA-NTvsHT. The TaCDPK21 expression pattern in 35DPA-NTvsHT was further verified through qRT-PCR ( Figure S7A ). Combining increased GP and decreased SD phenotypes, we suggest that TaCDPK21 negatively regulates SD and positively regulates germination.

To preliminarily investigate the role of TaCDPK21 in SD, we searched its Arabidopsis homologous AT2G31500 (AtCDPK24) on the basis of the protein sequence blast (http://plants.ensembl.org/index.html) and obtained the mutant atcdpk24 (SALK_015986C). This mutant had a transfer DNA (T-DNA) inserted in the third exon of AtCDPK24 and lacking AtCDPK24 expression ( Figure S7B ). According to germination test results, the GP (average GP 41%) of atcdpk24 mutant lines were lower than those of Col-0 (average GP 82%) ( Figure 3A ), demonstrating that AtCDPK24 has a negative role in SD. Considering the homologous relationship between them, we speculate that TaCDPK21 may have a function similar to that of AtCDPK24. Subsequently, we constructed a fusion with the full-length coding sequence of TaCDPK21 under the control of the constitutive cauliflower mosaic virus 35S promoter for overexpression in the atcdpk24 mutant background (35S:TaCDPK21/atcdpk24). Two independent transgenic lines (35S:TaCDPK21/atcdpk24-#3 and 35S:TaCDPK21/atcdpk24-#6) were selected for further studies ( Figures 3B , S7C ). The phenotypic data revealed that TaCDPK21 overexpression increased the GPs (average GP 79%) of the atcdpk24 mutant lines (average GP 41%), suggesting that TaCDPK21 negatively regulates SD and positively mediates germination ( Figure 3A ). Simultaneously, transgenic Arabidopsis plants ectopically expressing TaCDPK21 were also generated (35S:TaCDPK21-#2 and 35S:TaCDPK21-#12) ( Figures 3B , S7C ). According to the results, the GP of the 35S:TaCDPK21 plants (average GP 98%) were higher than those of Col-0 (average GP 82%). This finding supports the negative role of TaCDPK21 in SD and the positive role in germination ( Figure 3A ).

To authenticate the roles of TaCDPK21 in wheat SD, we identified a cdpk21 EMS mutant (JM22) with a point mutation (G/A) in the protein kinase domain site of TaCDPK21, resulting in premature termination of translation ( Figure 3C ). To eliminate the interference of the genetic background, the cdpk21 mutant was backcrossed with JM22 for two generations. The BC₂F₂ (cdpk21) mutant seeds showed lower GPs (average GP 37%) than JM22 seeds (average GP 75%) ( Figures 3D, E ). Moreover, compared with JM22 seeds (average GP 95%), cdpk21 seeds were less sensitive to HT (average GP 48%) ( Figures 3D, E ). The results confirm that TaCDPK21 participates in HT-mediated dormancy.

To preliminarily validate the role of miR27319 in SD, miR27319 expression patterns were first verified through qRT-PCR ( Figures 4A , S8A ). MiR27319 shared the high similar sequence between wheat and rice, and had three base differences between wheat and Arabidopsis ( Figures S9A, B ). Then the miR27319 precursor was heterologously overexpressed in Arabidopsis and rice, the mature transcripts of miR27319 in Arabidopsis and rice were simultaneously silenced using the STTM method.

Two overexpression lines (35S:pre-miR27319-#9 and 35S:pre-miR27319-#11) and two STTM lines of miR27319 (STTM-AtmiR27319-#1 and STTM-AtmiR27319-#2) in Arabidopsis were generated. The pre-miR27319-OE transcript levels were higher in overexpression lines than in Col-0, and the mature transcripts of miR27319 were lower in the STTM lines than in Col-0 ( Figure 4B ). Seed GPs of miR27319 overexpression lines and STTM lines of Arabidopsis were determined. According to the results, the GPs of 35S:pre-miR27319-#9 (average GP 33%) and 35S:pre-miR27319-#11 (average GP 34%) were lower than those of Col-0 (average GP 80%), whereas the GPs of STTM-AtmiR27319-#1 (average GP 98%) and STTM-AtmiR27319-#2 (average GP 94%) were higher than those of Col-0 (average GP 80%) ( Figures 4C, D ). This suggests the positive role of miR27319 in SD.

Similarly, 15 pre-miR27319 overexpression rice (Nip, WT) lines (35S:pre-miR27319) and 15 miR27319 STTM rice lines (STTM-OsmiR27319) were generated ( Figures 5A–C ). qRT-PCR analysis revealed that the primary transcripts of miR27319 were higher in the overexpression plants than in the WT lines ( Figure 5B ), and the mature transcripts of miR27319 were lower in the STTM lines than in the WT lines ( Figure 5C ). Compared with the WT lines (average GP 76%), 15 pre-miR27319 overexpression lines exhibited lower GPs (39%–63%, average GP 53%) (P < 0.05), whereas 15 STTM-OsmiR27319 lines exhibited higher GPs (80%–98%, average GP 87%) (P < 0.05) ( Figures 5D, E ). The linear regression analysis revealed that the relative expression of miR27319 in 35S:pre-miR27319 and STTM-OsmiR27319 lines had a significant negative correlation with GPs (R² = 0.78** and R² = 0.77**, respectively ( Figures S8B, C ), indicating that miR27319 positively regulates SD in rice.

Phytohormones GA and ABA are the core factors determining SD and seed germination in diverse plant species. GA promotes seed germination by elevating α-amylase activity, whereas ABA inhibits germination and induces dormancy (Ali-Rachedi et al., 2004; Kucera et al., 2005). To preliminarily explore the molecular mechanism of miR27319 in SD regulation in rice, we first measured the ABA and GA contents as well as α-amylase activity in the seeds of the representative 35S:pre-miR27319-#6, 35S:pre-miR27319-#11, STTM-OsmiR27319-#17, and STTM-OsmiR27319-#30 lines as well as in the WT seeds. The seeds of 35S:pre-miR27319-#6 and 35S:pre-miR27319-#11 lines had higher ABA levels and lower GA levels and α-amylase activity than the WT seeds ( Figures 6A–C ). The seeds of STTM-OsmiR27319-#17 and STTM-OsmiR27319-#30 lines showed the opposite trends ( Figures 6A–C ). Subsequently, we investigated the expression of several key genes involved in ABA and GA biosynthesis/catabolism/signaling pathways in the seeds of 35S:pre-miR27319 and STTM-OsmiR27319 as well as the WT plants. Three key genes involved in the ABA biosynthesis/signaling pathways (OsNCED1 and OsABI5) and the GA catabolism pathway (OsGA2ox1) exhibited upregulated expression in the seeds of 35S:pre-miR27319 and downregulated expression in the seeds of STTM-OsmiR27319 compared with the WT seeds ( Figure 6D ). By contrast, the expression of the GA biosynthesis gene OsGA20ox1 and the α-amylase gene OsAMY1B exhibited opposite trends ( Figure 6D ). The results suggest that miR27319 regulates SD through the ABA and GA biosynthesis/catabolism/signaling pathways.