The OsNAR2.1 overexpression transgenic line (Ov199) used in this study, which was named as pUbi-OsNAR2.1, have been described in detailed in previous studies (Chi et al., 2011; Chen et al., 2017). Briefly, we amplified the OsNAR2.1 ORF sequence from cDNA isolated from the Oryza sativa L. ssp. Japonica and ligated it into the expression vector containing the Ubi promoter. The expression construct was later transferred into Agrobacterium tumefaciens strain EHA105 by electro-poration, followed by transformation into the rice as described by Upadhyaya et al. (Fan et al., 2020). A previously described procedure was adopted to generate the OsNAR2.1 RNAi line (RNAi) (Yan et al., 2011).

In field experiments, the seedlings of OsNAR2.1 overexpression line (Ov199), OsNAR2.1 RNAi line (RNAi) and wild type (NP) were grown in plots at the Nanjing Agricultural University in Nanjing, Jiangsu. Soil chemical properties and fertilizer application details are given in (Upadhyaya et al., 2000).

In hydroponic experiments, the seeds were surface-sterilized and germinated in the International Rice Re-search Institute (IRRI) solution. The detailed nutrient composition and pH of IRRI have been reported previously (Chen et al., 2016). The lines of Ov199, RNAi and NP were grown in a greenhouse with IRRI solution for one week under 14 h light (30^°C)/10 h dark (22^°C) photoperiod and 60% relative humidity conditions. For different nitrogen conditions treatment experiments, the seedlings of NP were grown under nitrogen deficient conditions for three days, followed by treating with 0.125 mM NH₄NO₃ (LN), 1.25 mM NH₄NO₃ (NN) and 2.5 mM NH₄NO₃ (HN) for one week. The RNA samples of NP were extracted after the treatment.

The pooled RNA samples extracted from the leaves, flowers and stems of rice obtained during the blooming stage were sent for sequencing services by the Kangchen Bio-tech Inc., China. The quality and quantity of RNA samples were analyzed by NanoDrop ND-1000 model. The quantified RNA was used to synthesize cDNA, followed by the ligation of 5’- and 3’-adapters to create the sequencing library. The cDNA samples were diluted to a final concentration of 8 pM prior to generate clusters on Illumina cBlot using TruSeq SR Cluster kit (#GD-402-4001, Illu-mina). Sequencing was performed on Illimina HiSeq 2000 using TruSeq Rapid SBS Kit (#FC-402-4001, Illumina), followed by data acquisition and processing sequence. The sequence data have been deposited in the NCBI GEO under accession number GSE224933.

All unmatched sequences from the full small RNA sequencing data to the known miRNA annotations were analyzed with miRDeep2 package (http://www.mdcberlin.de/en/research/research_teams/systems_biology_of_gene_regulatory_elements/projects/miRDeep/) to construct potential models of miRNA precursors by exploring the Dicer cleavage sites in order to predict the novel miRNAs (Zhang et al., 2020).

The previously described protocols for total RNA isolation and concentration measurement were used with modifications (Friedländer et al., 2008; Yan et al., 2011). The synthesis of miRNA first-strand cDNA was performed with miRNA First-Strand cDNA Synthesis SuperMix (Vazyme, Co. R323-01, Nanjing, China), prior to the amplification of qRT-PCR products with the AceQ qPCR SYBR Green Master Mix kit (Vazyme Biotech Co. Q311-02, Nanjing China.) using a Step One Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The details of primers for the PCR are shown in Supplementary Table 5 .

The collected data were tabulated and analyzed for significant differences using the IBM SPSS Statistics 20 program and One-Way ANOVA, followed by Tukey’s test (P < 0.05).

In rice, OsNAR2.1 is a high affinity nitrate transporter partner protein that plays a critical role in the absorption of nitrate (Yan et al., 2011; Liu X et al., 2014). Previously, the transgenic rice plants p35S:OsNAR2.1 and pOsNAR2.1:OsNAR2.1 were found to promote plant growth and increase the N concentration, N content and yield of rice (Chen et al., 2017; Chen et al., 2020). Also, overexpression of OsNAR2.1 was found to enhance drought tolerance and grain yield under drought stress conditions in rice (Chen et al., 2019). In order to explore the regulatory mechanism of OsNAR2.1 in rice growth, the agronomic traits of the pUbi : OsNAR2.1 overexpression line (Ov199), OsNAR2.1 RNAi suppression line (RNAi) and the non-transgenic line (NP) were investigated ( Figure 1 ). The plant height was substantial increased in Ov199 lines, but was reduced in RNAi lines as compared to NP ( Figures 1A, B ). In addition, the study reported that the seeds weight per panicle, the number of seeds per panicle, the seed setting rate and the yield per plant in Ov199 were significantly increased by 57.3%, 63.8%, 20.2% and 35.0%, respectively than the NP line. Similar agronomic traits were significantly reduced by 24.8%, 27.0%, 32.1% and 51.1%, respectively in RNAi line as compared to the NP (Liu X et al., 2014). And the dry weight of Ov199 was increased compared with NP, and reduced in RNAi lines ( Figure 1C ). Total N content of Ov199 was substantial increased in Ov199 lines compared with NP, but opposite in RNAi lines ( Figure 1D ). Hence, our findings strongly suggest that OsNAR2.1 plays a vital role in the rice growth and the output in rice production.

In order to explore whether the overexpression or silencing of OsNAR2.1 affect the expression of small RNAs, we performed small RNA sequencing by using Ov199, RNAi and NP lines. After excluding the low-quality reads and adaptors, the align clean reads of the small RNAs for the Ov199 and RNAi lines were found to be higher in quantity than those of NP ( Supplementary Table 1 ). Further analysis of the retained small RNA from all lines have found that the prominent sizes of the small RNAs left were 21-25 in length ( Figure 2A ). The 24 nt small RNAs accounted for 25.21%, 24.05%, 26.42% of total small RNAs in NP, Ov199 and RNAi line, respectively. The proportion of 24 nt small RNA in Ov199 was lower than that of NP and RNAi ( Figure 2A ). Meanwhile the 21 nt small RNAs accounted for 7.31%, 8.02%, and 7.41% of total small RNAs in NP, Ov199 and RNAi, respectively. The proportion of 21 nt small RNAs in Ov199 was higher than that of NP and RNAi ( Figure 2A ).

The microarrays analysis of miRNAs from the three experimental lines revealed that one hundred sixty-five known miRNAs belong to the seventy-six miRNA families were differentially expressed ( Figure 2B ). The number of known miRNA family members varies greatly on the basis of the differences in endogenous nucleotides. The miR166 family has the highest number of members, which can be classified into seventeen different miRNA families ( Figure 2B ). Besides, the miR156, miR167, miRNA444 and miR1862 families are groups that contain 5-15 members. The rest of the known miRNAs fall into twenty-three miRNA families each containing 2-4 members and forty single-membered miRNA families ( Figure 2B ).

The expression abundance of miRNAs in NP, Ov199 and RNAi was different based on the normalized transcript per million (TPM) analysis of miRNA data. Through hierarchical clustering and expression analysis, we found that the expression level of 45% of miRNAs in Ov199, 52% of miRNAs in RNAi and 35% of miRNAs in NP from was upregulated ( Figure 3 ; Supplementary Table 2 ). Among the known miRNAs, miR1876 was the only miRNA upregulated/expressed in all experimental lines. The number of transcription units for NP, Ov199 and RNAi were 36,075, 32,893 and 43,939, respectively. The expression level was highest in the miR166 family, followed by miR167, miR444, miR168 and the remaining miRNA families ( Supplementary Table 2 ). Some miRNAs grouped under the same family were found to be differentially expressed in different experimental lines ( Supplementary Table 2 ). For example, miR156a has the lowest expression in NP, but the miRNA with the lowest expression in Ov199 and RNAi was miR156b-3p and miR156j-3p, respectively ( Supplementary Table 2 ).

In order to identify the N-responsive miRNAs, we comparatively analyzed the expression of miRNAs from all the experimental lines. The comparative results from three comparison groups (i.e., Ov199/RNAi, Ov199/NP, RNAi/NP) showed that there were 62, 64 and 59 differentially expressed miRNAs in each respective group with P <0.05 significant difference ( Figure 4A ; Table S2 ). Nine out of all differentially expressed miRNA were found to be significantly expressed in all comparison groups, while majority of them only altered in two comparison groups. There were 3, 4 and 5 miRNAs with significant different levels exclusively expressed in group Ov199/RNAi, Ov199/NP and RNAi/NP, respectively ( Figure 4A ).

Differential expression analysis of miRNA between transgenic lines and the control (NP) showed that only the peak expression of miR169f in Ov199 and RNAi was consistent with the changes in OsNAR2.1 expression and total nitrogen concentration ( Figure 4B ; Supplementary Table 2 ). The expression peak of miR529a, miR1882e-3p and miR5150-5p in Ov199 and RNAi were opposite to the expression level of OsNAR2.1 and total N concentration ( Figure 4B ; Supplementary Table 2 ). The expression patterns of miR169f, miR529a and miR1882e-3p were consistent with the sequencing results as verified by RT-PCR ( Figure 5 ). The expression level of miR169f was significantly upregulated in Ov199 and downregulated in RNAi as compared to the NP ( Figure 5A ). In contrary, the expression of miR529a and miR188e-3p were significantly downregulated in Ov199 and upregulated in RNAi ( Figures 5B, C ). Therefore, OsNAR2.1 changes the nitrogen concentration and content in rice and maybe affect the expression of multiple miRNAs.

Identification of novel miRNAs from the transgenic OsNAR2.1 lines was performed with miRDeep analysis package. A total of 150 novel miRNAs were identified from all experimental lines with an approximate length of 17-23 nt, of which the most common length falls in 24 nt ( Supplementary Table 3 ). There were 13, 32 and 10 specific novel miRNAs that only expressed in NP, Ov199 and RNAi, respectively; while 54 novel miRNAs were found in all the three lines ( Figure 6A ). Besides, there were 64 novel miRNAs found to be overlapped in group NP/Ov199, whereas 66 and 73 novel miRNAs overlapped in group NP/RNAi and Ov199/RNAi, respectively ( Figure 6A ). Compared with NP, 22 novel miRNAs were upregulated in Ov199 but 21 of them were downregulated in RNAi ( Figure 6B ). By referring to the known miRNA families, the abundance of novel miRNAs was very low, and more than 75% of the novel miRNAs have less than 100 transcript units ( Supplementary Table 3 ). There were only five novel miRNAs with high abundance of more than 1000 transcript units that is miR-novel-chr1-37617, miR-novel-chr1-38364, miR-novel-chr1-43165, miR-novel-chr4-27017 and miR-novel-chr8-11525 ( Supplementary Table 3 ). And the prediction target genes of miR-novel-chr1-38364, miR-novel-chr1-43165 and miR-novel-chr4-27017 are LOC_Os09g3939.2, LOC_Os09g10820.1, LOC_Os08g9080.2, respectively, and others no corresponding target gene ( Supplementary Table 4 ). All these five highly abundance novel miRNAs expressed differentially in all the three lines. In general, the expression level of novel miRNAs was relatively downregulated in RNAi compared to the NP and Ov199 ( Supplementary Table 3 ). The Ov199 line has the highest number of novel miRNAs (76.7%) from its total number of miRNAs, which was much higher than the NP (59.3%) and RNAi (63.3%) lines. The results implied that the overexpression of OsNAR2.1 gene in Ov199 line might help to promote the expression of novel miRNA in the rice plants.

In order to examine the effects of different nitrogen concentrations in the regulation of miRNA expression, we performed miRNA expression verification on rice under high nitrogen (HN), low nitrogen (LN) and normal nitrogen (NN) treatments. A total of five miRNAs were selected as the candidates for the verification purposes that is miR1874, miR5150, chr3-36174, chr4-27017 and chr5-21745 ( Figure 7 ). The results of qRT-PCR showed that the expression of miR1874 was inhibited in both HN and LN treatments as compared to the NN. The inhibition on the chr3-36174 was also observed in HN, but was not in LN ( Figure 7 ). The expression of miR5150, chr4-27017 and chr5-21745 was up-regulated under HN treatment in relative to the NN, but was not significantly different under LN condition ( Figure 7 ).

The nitrogen concentration of rice can be altered with the expression of OsNAR2.1, which can result in the epigenetic changes in rice, including the expression pattern and intensity of miRNAs. It has been reported that miR169 in maize, common bean (Pharsalus vulgaris) and Arabidopsis responded to N starvation (Valdés-López et al., 2010; Xu et al., 2011; Zhao et al., 2011; Liang et al., 2012; Trevisan et al., 2012). Besides, the miR393/AFB3 has been linked to the regulation of root system architecture in response to the internal and external N availability in Arabidopsis (Vidal et al., 2010). The lateral root growth is usually associated with the nitrate uptake, but the overexpression of miR393 was found to alter the primary and lateral root growth after nitrate application (Vidal et al., 2010). Recently, an auxin-related regulatory miR167/ARF8 module has been reported to regulate the ratio between emerging and initiating lateral roots in Arabidopsis. The cell-specific regulation of a transcriptional circuit mediates lateral root growth in response to the N via miRNA167 (Gifford et al., 2008). Furthermore, the expression of recently described miR529 was found to render some agronomic traits of rice, such as rice tiller number, panicle type and hydrogen peroxide resistance. The miR529 targets five SPL genes (OsSPL2, OsSPL14, OsSPL16, OsSPL17 and OsSPL18), where the OsSPL14 expression showed to be inversely correlated to the miR529a (Yue et al., 2017). The OsSPL14 has been shown to improve the nitrogen uptake in rice with an expression pattern oppositely responses to the level of nitrogen supply (Srikanth et al., 2016; Wu et al., 2020). Studies shown that the endogenous N content and concentration of rice via transgenic OsNAR2.1 lines to determine miRNAs that could be responding to the changes of endogenous N ( Figures 1 , 2 , 4 ). Therefore, it is speculated that the expression peaks of these miRNAs might be correlated to the expression of OsNAR2.1 and the nitrogen concentration in rice.

High-throughput sequencing approach has been using in many studies to analysis genome-wide miRNA expression in plants (Zhao et al., 2007; Chi et al., 2011; Liang et al., 2012). Apart from that, transgenic technology provides a permissible alternative for in vivo study to examine the effect of endogenous nitrogen concentration of rice towards the miRNA expression, which allows further attempts to improve the understanding between miRNA and the rice NUE. Nitrogen-deficient conditions cause plants to deplete N internally. As an emergency response, they regulate the ability to absorb and transport N by regulating the genes and small RNA expression. The report shown that during nitrification and absorption, miR393 is activated by the N signal transmitted, and effected the function of the phytohormone sensing organ auxin signaling F-box proteins (AFB3) (Bao et al., 2014; Sajjad et al., 2021). In maize plants deficient in N, six miRNAs were downregulated, including miR408, miR169, miR166, miR528*, and miR169* (Yang et al., 2019).

In our study, the number of miRNAs being expressed was highest in the nitrogen-efficient line (Ov199), followed by the nitrogen-inefficient line (RNAi) and the wild-type (NP), regardless of the known or novel miRNAs ( Supplementary Table 1 ). Among the 159 known miRNAs, both transgenic lines have about the same number of the highest expressing miRNAs, i.e., 45% and 52% of known miRNAs in Ov199 and RNAi, respectively. Meanwhile, there was only 35% of the known miRNAs with the highest expression in NP. The results indicate that OsNAR2.1 may play a role in the transgenic lines to alter the rice endogenous nitrogen environment, resulting in the differential expression pattern of a large number of miRNAs in rice.

Previously, many known miRNAs have been reported to be involved in nitrogen related metabolisms in plants. The nitrogen concentration which was alter in plants may affect the expression of miRNA. More than 40 miRNAs, such as miR164, miR167, and miR399, are responsive towards nitrogen starvation in different parts of maize (Xu et al., 2011; Trevisan et al., 2012). The miRNAs from 7 families, including miR156, miR157, miR399, participate in the regulation of NUE in different parts of rice (Cai et al., 2012). A total of 19 miRNAs families, such as miR156, miR157, miR160, are involved in the regulation of NUE in Arabidopsis (Zhao et al., 2011; Liang et al., 2012). Also, there are 25 miRNAs in soybeans involved in the regulation of soybean NUE, such as miR396, miR164, miR168 (Valdés-López et al., 2010).

By comparing the miRNA sequencing data from all the experimental lines, we have identified four miRNAs (miR169f, miR529a, miR1882e-3p and miR5150-5p) expressed in a way that significantly inversed between the transgenic lines and the wild-type ( Figure 5 ). Among them, miR169 has been reported to respond to the nitrogen deficient conditions in corn (Xu et al., 2011; Trevisan et al., 2012), soybean (Liang et al., 2012), and Arabidopsis (Zhao et al., 2011; Liang et al., 2012). The protein-coding gene NFYA, which is the target gene of miR169, encodes for the nitrate transporter proteins NRT1.1 and NRT2.1 as reported before (Zhao et al., 2011). MicroRNA529 (miR529) targets the OsSPL14 gene by reducing its expression level with the increase of nitrogen supply, showing its important roles in nitrogen metabolism of rice (Srikanth et al., 2016; Wu et al., 2020). The results showed the involvement of miR169 and miR529 in rice NUE regulation. To date, there is no direct evidence to support the expression of miR1882e-3p and miR5150-5p being related to the NUE in rice. In addition, some miRNAs can respond to changes in endogenous nitrogen as well as changes in exogenous nitrogen, such as miR1874, miR5150, chr3-36147, chr4-27017 and chr5-21745 ( Figures 6 , 7 ). By changing the NUE in rice, we can find novel mature miRNAs related to nitrogen, which will be useful for future research on NUE-related miRNAs.

In summary, the manipulation of rice endogenous nitrogen concentration through genetic modification can be an effective approach for the discovery of novel N-responsive miRNAs. Thus, the current research provides new insights for the detection of nitrogen related novel miRNAs.