Seeds of the rice cultivar (Oryza sativa cv. Nanhui 511 (NH511)) were germinated in Murashige and Skoog (MS) media for 3 d and then transferred to tap water for 3 d before being transferred into a hydroponic solution containing KNO₃ as the only N source. The seedlings were then divided into two groups and transferred to HN or LN conditions for six days of cultivation. The seedlings were cultivated in a growth chamber under a photoperiod of 12 h/12 h (light/dark) (~230 μmol m^-2 s^-1) at 28°C/25°C, and the solution was renewed every 3 d. For the different nitrogen treatment assays and RNA-seq, the seedlings under normal nitrogen (1 mM KNO₃) were transferred to LN (0.2 mM KNO₃) and HN solutions (5 mM KNO₃) for 3 hours, respectively. Seedlings from different N treatments were collected and frozen in liquid nitrogen for RNA extraction. For time-course expression analysis under HN treatment, the seedlings were transferred to HN solutions (5 mM KNO₃) and RNA samples were collected at different time points under high nitrogen conditions (5 mM nitrate) by TRIzol reagent (Invitrogen, America) according to standard procedures.

Total small RNA was isolated from the samples treated for 3 h and grown under HN and LN conditions using the mirVana™ RNA Isolation Kits (Thermo Fisher, Vilnius, Lithuania), with three biological replications used for this assay. Subsequently, the RNA samples from these replicates were delivered to Biomarker Technologies (Beijing, China) for library construction and small RNA sequencing, after which transcriptome data were obtained and analyzed. The data content has been submitted to GEO database on 2023.4.18, and been publicly released on 2023.4.21. To review GEO accession GSE230023, the following links allow you to view the uploaded data: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE230023.

A total of 2 μg of RNA was used for reverse transcription according to standard procedures. The Vazyme miRNA 1st Strand cDNA Synthesis Kit (by stem-loop) (Vazyme, Nanjing, China) was used to synthesize miRNA cDNA. qPCR analysis and experimental data were obtained using a Roche Light Cycler 480 (Roche, Switzerland), after which the expression variation in the N-related genes was calculated. The primers used in this study are listed in Supplementary Table S5 .

Differential expression analyses of the two conditions were performed using the DESeq R package (version 1.10.1) (Love et al., 2014). DESeq provides statistical routines for determining differential expression in digital miRNA expression data, using a model based on a negative binomial distribution. The resulting P values were adjusted using Benjamini and Hochberg’s approach to control for the false discovery rate. miRNAs with an adjusted p<0.05 found by DESeq were assigned as differentially expressed.

Known and novel miRNAs were identified using the miRDeep2 software (Friedlander et al., 2012), which is a comprehensive software package for miRNA identification that matches reads to the genome based on the distribution information of reads on the precursor sequence. Target genes of the total miRNAs were predicted using the TargetFinder software (Allen et al., 2005). Possible targets of miR166k-3p and miR444b.1 were predicted using psRNATarget (https://www.zhaolab.org/psRNATarget/) under the condition of an expectation value of less than 3.5, respectively.

Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was performed using the GOseq R package based on the Wallenius non-central hyper-geometric distribution. The structure of miRNAs were predicted by RNAfold WebServer (http://rna.tbi.univie.ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi) and the targets degradome data were collected and analyzed from TarDB (http://www.biosequencing.cn/TarDB/).

Our previous study showed that rice cultivars exhibit different morphological and physiological variations when grown under different nitrogen conditions. In the present study, a typical indica rice restorer cultivar, namely NH511, was selected for further analysis. To fully understand its different growth responses to different nitrogen supplies, the plants were first germinated and cultivated under normal nitrogen conditions for five days. The seedlings were then divided into two groups and transferred to HN or LN conditions for six days of cultivation. We found that NH511 was sensitive to HN supply conditions, exhibiting promoted growth and increased lateral root numbers and lengths under HN conditions than under LN conditions ( Figures 1A, B ). Moreover, the lateral root tips of NH511 grown under HN and LN conditions were stained with KI-I₂ and observed under an optical microscope, and our observations further suggested that lateral root development was impaired by influencing root tip maturation under LN conditions ( Figures S1A, B ). Subsequently, transcriptomic read data were obtained from seedlings under HN and LN conditions by small RNA sequencing ( Table S1 ), and the statistical data of total reads and unique reads were analyzed. The results showed that the mutually unique reads accounted for 10.28% of unique reads under HN and LN conditions, 54.92% under HN, and 34.26% under LN ( Figure 1C ). For the total reads, the mutual reads accounted for 62.62% of total reads under HN and LN conditions, 23.04% under HN, and 14.34% under LN ( Figure 1D ). Taken together, these transcriptomic data demonstrated that NH511 exhibits favorable variations in miRNA profiles in response to different nitrogen supplies.

Furthermore, we investigated the miRNA profile and characteristics of NH511 in response to different nitrogen conditions. First, the RNA samples extracted from seedlings grown under different nitrogen conditions were sequenced, and the raw sequencing data were processed to obtain unique reads. Furthermore, the filtered data were further processed to facilitate miRNA alignment, identification, and prediction based on the rice genome. A total of 586 miRNAs were identified under HN conditions, including 458 known and 128 novel miRNAs, whereas 572 miRNAs were identified under LN conditions, including 444 known and 128 novel miRNAs ( Figure 2A ). Statistical analysis of the length distribution of clean data was performed based on the analysis and statistics of the original sequencing data. The results showed that in response to nitrogen, most of the miRNAs were distributed within 20–24 nt, with miRNAs with the length of 21 nt and 24 nt occupying a large proportion, which is consistent with the typical length range of miRNAs in plants ( Figure 2B ). Moreover, owing to post-transcriptional base editing of miRNAs, which leads to seed sequence changes and thus changes of target genes, the sequence variations of miRNA were analyzed using isomiRID software (Data sheet2). The target genes of the total miRNAs were predicted and analyzed using the TargetFinder software ( Table S2 ). Finally, the miRNA nucleotide bias at each position and the first nucleotide bias for all miRNAs were analyzed, and the results showed that most miRNAs started with A or U at the first nucleotide; A and U also occupied a large proportion of nucleotides at each position ( Figures S2A, B ).

The observed differences in the performance of NH511 under HN and LN conditions prompted us to identify differentially expressed miRNAs under different nitrogen conditions. Therefore, the miRNA differential expression profiles were obtained from the transcriptome data based on DESeq, which provides statistical routines for determining differential expression in digital miRNA expression data using a model based on negative binomial distribution. Our results showed that a total of 100 miRNAs were differentially expressed between HN and LN conditions, including 75 upregulated and 25 downregulated miRNAs ( Figure 3A ). Hierarchical cluster analysis was also performed for the selected differentially expressed miRNAs under HN and LN conditions ( Figure S3 ). Moreover, the Volcano Plot also showed differences in miRNA expression levels and showed that the difference in miRNA expression between the HN and LN conditions was statistically significant ( Figure 3B ). Additionally, small RNA sequencing was used to accurately quantify miRNA abundance in the rice restorer cultivar NH511, and the differentially expressed miRNAs with a relative change greater than 2-fold were identified and defined as N-responsive miRNAs ( Figures 3C and S4 ).

Based on the above results, a gene ontology (GO) enrichment analysis was performed to further characterize the main biological functions of target genes of differentially expressed miRNAs under high nitrogen supply. All differentially expressed miRNAs were divided into three categories: biological processes, cellular components, and molecular functions. These three categories were further divided into 42 enriched subcategories ( Figure 4A ). To fully understand the cellular pathways of target genes differentially expressed miRNAs under HN and LN conditions, an enrichment analysis based on KEGG was performed. Three KEGG pathway categories were identified: cellular processes, metabolism, and genetic information processing ( Figure 4B ). Specifically, ‘protein processing in endoplasmatic reticulum’ was the most enriched in cellular processes, which could be a consequence of different nitrogen supplies. Regarding metabolism, ‘phenylpropanoid biosynthesis’ was the most overrepresented because phenylpropanoids include a variety of organic compounds synthesized from amino acids such as phenylalanine and tyrosine, which are related to nitrogen metabolism. For genetic information processing, spliceosome, ribosome, and ribosome biogenesis are the top three pathways, which probably influences the gene expression and protein coding involved in nitrogen uptake and assimilation.

To further validate the differential miRNA expression under HN and LN conditions, a qPCR expression analysis was performed, and the nine selected typical miRNAs (miR1441, miR443, miR812q, miR5801, miR1861b, miR166k-3p, miR530-5p, miR444b.1, and miR395v) were validated using qPCR primers under HN and LN conditions ( Table S5 ). The results showed that most miRNAs had an expression pattern similar to that of their small RNA sequencing profiles ( Figures 5A–I ). Furthermore, three miRNAs (miR443, miR1861b, and miR166k-3p) exhibited a significant upregulation pattern in response to HN treatment, but were downregulated expressed in the LN treatment ( Figures 5B, E, F ). The expression of miR444b.1 and miR395v was downregulated under HN conditions but upregulated under LN conditions ( Figure 5H, I ). These data provide a foundation for the functional elucidation of candidate miRNAs and the possibility of selecting elite rice restorers for high-NUE rice breeding.

In the present study, differentially expressed miRNAs were identified from transcriptomic data under HN and LN conditions, and partially representative miRNAs were validated by qRT-PCR analysis. Based on qRT-PCR results and previous studies on identification and characterization of targets genes of these validated miRNAs, we selected miR166k-3p and miR444b.1 for the further investigation. The two candidates, miR166k-3p and miR444b.1, were characterized and shown to be upregulated and downregulated under HN conditions, respectively. First, possible targets of miR166k-3p and miR444b.1 were predicted using psRNATarget (https://www.zhaolab.org/psRNATarget/) under the condition of an expectation value of less than 3.5, respectively ( Tables S3 and S4 ). As shown by the results, miR166k-3p mainly regulated homeodomain containing proteins and other unknown expressed proteins. Homeodomain containing proteins are a large family of transcription factors (TFs) that contain a highly conserved DNA-binding domain of 60 amino acids known as the homeodomain. They mediate protein-protein or protein-DNA interactions during plant growth and development. Furthermore miR444b.1 regulated MADS-box family proteins, which have been reported to be involved in many important pathways of environmental adaptation in plants. Furthermore, the secondary structures of miR166k-3p and miR444b.1 were predicted and analyzed for the first time by RNAfold WebServer (http://rna.tbi.univie.ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi) ( Figure S5 ), and the degradomes of crucial targets were predicted and identified suing the TarDB database (http://www.biosequencing.cn/TarDB/) for miR166k-3p ( Figure 6 ) and miR444b.1 ( Figure 7 ) respectively. These results provided detailed information on the target gene processing by the two candidate miRNAs as well as a foundation for further studies.

To further investigate the expression variations in miR166k-3p and miR444b.1 under high nitrogen supply, a qRT-PCR assay was carried out. The expression levels of miR166k-3p with its target gene OsHOX10 (LOC_03g01890) and miR444b.1 (LOC_02g36924) with its target gene OsMADS27 under the HN treatment were validated using qPCR primers at different time points. The results showed that miR166k-3p expression was upregulated in response to HN condition, which was reversed in the expression pattern of the target gene OsHOX10 ( Figure 8A ). Similarly, the expression of miR444b.1 was suppressed, and the expression of its target gene OsMADS27 was induced under the HN treatment ( Figure 8B ). Taken together, these results further elucidated the response of the two candidate miRNAs (miR166k-3p and miR444b.1) to nitrogen supplies as well as their different expression patterns, which provides possibilities for further studies of these two candidate miRNAs and new insights for high-NUE rice breeding using elite restorer cultivars.