Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Identification of microRNAs in developing wheat grain that are potentially involved in regulating grain characteristics and the response to nitrogen levels

BMC Plant Biology volume 20, Article number: 87 (2020) Cite this article

Abstract

Background

MicroRNAs (miRNAs) play crucial roles in the regulation of plant development and growth, but little information is available concerning their roles during grain development under different nitrogen (N) application levels. Our objective was to identify miRNAs related to the regulation of grain characteristics and the response to different N fertilizer conditions.

Results

A total of 79 miRNAs (46 known and 33 novel miRNAs) were identified that showed significant differential expression during grain development under both high nitrogen (HN) and low nitrogen (LN) treatments. The miRNAs that were significantly upregulated early in grain development target genes involved mainly in cell differentiation, auxin-activated signaling, and transcription, which may be associated with grain size; miRNAs abundant in the middle and later stages target genes mainly involved in carbohydrate and nitrogen metabolism, transport, and kinase activity and may be associated with grain filling. Additionally, we identified 50 miRNAs (22 known and 28 novel miRNAs), of which 11, 9, and 39 were differentially expressed between the HN and LN libraries at 7, 17, and 27 days after anthesis (DAA). The miRNAs that were differentially expressed in response to nitrogen conditions target genes involved mainly in carbohydrate and nitrogen metabolism, the defense response, and transport as well as genes that encode ubiquitin ligase. Only one novel miRNA (PC-5p-2614_215) was significantly upregulated in response to LN treatment at all three stages, and 21 miRNAs showed significant differential expression between HN and LN conditions only at 27 DAA. We therefore propose a model for target gene regulation by miRNAs during grain development with N-responsive patterns.

Conclusions

The potential targets of the identified miRNAs are related to various biological processes, such as carbohydrate/nitrogen metabolism, transcription, cellular differentiation, transport, and defense. Our results indicate that miRNA-mediated networks, via posttranscriptional regulation, play crucial roles in grain development and the N response, which determine wheat grain weight and quality. Our study provides useful information for future research of regulatory mechanisms that focus on improving grain yield and quality.

Background

As one of the largest and most widely planted food crop species in the world, wheat provides approximately 20% of the required human calorie intake [1]. With the rapid increase in the global population, increasing wheat yields is a major goal to meet future food demands [2]. Additionally, with increasing economic development, improving wheat quality, including increasing the protein content, is another important objective for wheat production. Conventional breeding approaches, which depend on the manipulation of genetic variation, have been successful at producing wheat cultivars with improved grain yield and quality; however, understanding the molecular mechanisms that regulate wheat development and grain quality would benefit further wheat improvement programs via modern genetic engineering technology [3].

Small RNAs (sRNAs) are ubiquitous components of plant transcriptomes; they have been found to regulate cellular metabolism, differentiation and growth and are involved in defense against viruses and mobile genetic elements in plants [4]. MicroRNAs (miRNAs) are single-stranded, noncoding sRNAs produced from endogenous transcripts with an imperfect self-complementary stem–loop structure [5, 6]. miRNAs play crucial regulatory roles in organisms by targeting specific mRNAs for cleavage or translational repression [7,8,9]. It has been reported that miRNAs affect various biological processes, such as the development of leaves, roots, and flowers, and affect grain filling [7, 10,11,12]. Grain filling is a critical phase in crop grain development. In rice, osa-miR167 regulates the auxin-miR167-ARF8-OsGH3.2 pathway, which functions during grain filling [13]. High-throughput sequencing revealed that the differentially expressed miRNAs between superior and inferior wheat grain are likely related to cell division, carbohydrate metabolism and hormone biosynthesis [14]. Han et al. [15] demonstrated that miR164, miR160, and miR169 show different expression profiles during grain filling, suggesting that their functions are coordinated during the different stages of wheat seed development. Similarly, Meng et al. [12] identified 873 miRNAs from wheat grain and suggested that many of them potentially regulate grain filling. Using a genome-wide investigation of miRNAs in different wheat tissues, Sun et al. [16] characterized 64 miRNAs that are expressed in developing wheat grain and that might play important functions in regulating grain development.

Furthermore, miRNAs that respond to nitrogen (N) deficiency in higher plants have been identified in previous studies [17,18,19]. In maize, miR169, miR399, miR518 and miR408 respond to N in the roots and leaves, but miR160, miR167, miR168 and miR395 respond to N deficiency only in the roots [18]. Several miRNAs, including miR167, miR171, miR398, miR827, miR408, and miR857, have been identified as being responsive to low-N conditions [20], and the wheat miRNA TamiR444a, which responds to nitrogen deficiency, is involved in plant adaptation to nitrogen starvation stress [21]. Zuluaga et al. [22] identified 84 novel and 161 conserved miRNAs in durum wheat in response to N starvation and found that ttu-miR169c and ttu-novel-61 were strongly downregulated under N-starvation conditions. Some miRNAs, such as miR159a, miR159b and miR399, showed significantly different expression profiles under low N treatment compared with high N treatment [23]. However, little information is available concerning the regulatory role of miRNAs in wheat grain in response to N deficiency, especially at different developmental stages.

Nitrogen is extremely important for plant growth and for agricultural productivity of many crop species [24]. Together with other N metabolites, nitrate can also act as a signal that regulates global gene expression [25]. In China, applying N fertilizer is an important agronomic practices to increase crop yield and economic productivity [26]. In addition to its effect on yield, N fertilizer is one of the key agronomic means of increasing grain protein content and improving wheat grain quality [27, 28]. However, excessive applications of N fertilizer have resulted in environmental pollution [29] and soil acidification [30]. Many efforts have been made to improve N fertilizer-use efficiency to reduce the waste of resources and to solve eco-environmental problems. Understanding the molecular mechanism of plant responses to N fertilizer applications might enable the development of crop varieties with increased nutrient-use efficiency and high grain yield and quality. To gain further insight into the role of wheat miRNAs that respond to different N fertilizer levels wheat grain development, the wheat cultivar Zhengmai 119 (ZM119), which has a high protein content, was grown under high nitrogen (HN) and low nitrogen (LN) conditions, and miRNA expression patterns were evaluated during grain development. The findings in this study will provide new insights into the regulatory function of wheat miRNAs that not only are associated with grain yield and quality but also respond to N supply levels.

Results

Grain characteristics

As shown in Table 1, the wheat cultivar ‘ZM119’ grown under HN conditions had a high grain yield, thousand-kernel weight, total protein content and protein fraction content (albumin, globulin, gliadin, and glutenin), whereas wheat grain from plants grown under LN conditions had low values. Correspondingly, the dough (rheological) characteristics of the wheat flour from the HN treatment had higher values than did the dough from the LN treatment. These results confirm that, compared with no N treatment, N fertilizer applications can increase grain yield and improve wheat grain quality.

Table 1 Grain characteristics of the wheat cultivar ‘ZM119’ grown under HN and LN conditions
Full size table

Deep sequencing of sRNAs in developing wheat grain from plants grown under two N fertilizer levels

To investigate the role of miRNAs in the developing wheat grain under different N fertilizer application levels, we sampled wheat grain at 7, 17, and 27 days after anthesis (DAA) under LN and HN conditions (three biological replicates each) and constructed a total of 18 sRNA libraries (Additional file 1: Table S1). The mean data of the three biological replicates are listed in Table 2.

Table 2 Sequencing data of the small RNA libraries constructed from developing wheat grain produced under two nitrogen fertilizer levels
Full size table

After sequencing, 10,475,211 to 11,752,733 redundant raw reads were produced per library (Table 2). After removing the low-quality reads, 3,955,306 to 6,430,010 clean reads were obtained from each library, representing 36.16–60.05% of the total redundant raw reads. The number of unique raw reads ranged from 1,561,169 to 2,014,214, and the number of unique clean reads ranged from 891,960 to 1,412,072, representing 52.37–70.00% of the total unique raw reads. The majority of the miRNAs in the constructed libraries were between 21 and 24 nt in length (Fig. 1), with the 21-nt class (59.83–74.34%) being the most frequent of the total reads, followed by the 24-nt miRNAs (7.08–14.64%). The 21-nt miRNAs were more frequent in the HN-7 and LN-7 libraries, whereas the 24-nt miRNAs were well represented in the LN-17 and HN-27 libraries.

Fig. 1
figure 1

Length distribution of miRNAs in developing wheat grain of plants grown under two nitrogen application levels and sampled at 7, 17, and 27 DAA. HN-7, HN-17 and HN-27 represent grain at 7, 17, and 27 DAA under high nitrogen application levels, respectively. LN-7, LN-17 and LN-27 represent grain at 7, 17, and 27 DAA under low nitrogen application levels, respectively

Full size image

Identification of known miRNAs and predicted novel miRNAs

To identify the conserved miRNAs in the constructed libraries in this study, the unique reads from each library were subjected to homology searches with known mature miRNAs from plants within miRBase (Additional file 2: Table S2). Additionally, novel miRNAs were predicted from unannotated sequences without similarity to known miRNAs (Additional file 2: Table S2, Additional file 3: Table S3). Since many miRNAs have not been deposited into miRBase, the novel miRNAs identified in this study were also compared with previously published wheat miRNAs [12, 14,15,16, 31,32,33,34,35]. In total, the sequences of 39 miRNAs from among the newly identified miRNAs were the same as those of previously found novel miRNAs (Table S2; the miRNAs are labeled in yellow). Most of the novel miRNAs had a relatively low expression level. The most abundant novel miRNAs were tae-miR9655-p5 and tae-miR2916-p3_2ss2TC17CA, which accounted for 24,656 reads and 27,735 reads in the LN-17 and LN-27 libraries, respectively. A total of 205, 327, and 322 unique miRNAs (of which 92, 190, and 191 were novel, respectively) were found in the HN-7, HN-17, and HN-27 libraries, respectively (Fig. 2). Additionally, a total of 143 miRNAs (55 novel ones) were found in the three HN-7, HN-17, and HN-27 libraries. The HN-7 and HN-17 libraries shared the same 19 miRNAs, of which 11 were novel. There were 245, 332, and 233 miRNAs (of which 121, 199, and 118 were novel, respectively) expressed in the LN-7, LN-17, and LN-27 libraries, respectively. We found that 146 miRNAs (including 53 novel ones) were expressed in all three LN libraries. When the effects of the two N fertilizer application levels were compared, LN-7 and HN-7 were found to share 189 expressed miRNAs, including 83 novel miRNAs. LN-17 and HN-17 shared 283 expressed miRNAs, of which 163 were novel. It is clear that the 17 DAA libraries (HN-17 and LN-17) contained more expressed miRNAs than did the other libraries.

Fig. 2
figure 2

Venn diagram of differentially expressed miRNAs in different libraries at three grain developmental stages under high (HN) and low (LN) nitrogen treatment. HN-7, HN-17 and HN-27 represent grain at 7, 17, and 27 days after anthesis under high nitrogen application levels, respectively. LN-7, LN-17 and LN-27 represent grain at 7, 17, and 27 days after anthesis under low nitrogen application levels, respectively. The number outside (inside) the brackets represents the number of total (novel) differentially expressed miRNAs

Full size image

Differentially expressed miRNAs in the developing wheat grain

To identify the differentially expressed miRNAs in developing wheat grain, we analyzed the relative abundance of miRNAs in the HN-7, HN-17, and HN-27 libraries, and P ≤ 0.05 was used as a threshold to determine significant changes in expression during grain development. Similarly, the abundance of miRNAs in LN-7, LN-17, and LN-27 was also compared. A total of 179 miRNAs (91 known and 88 novel miRNAs) were identified in the comparisons of the three HN treatment libraries (Additional file 4: Tables S4–2 and S4–4). Additionally, 80 known and 76 novel miRNAs were identified among the three LN treatment libraries (Table S4–1, S4–3). All of these miRNAs belong to 51 miRNA families. As shown in Additional file 4 (Table S4–1), the known miRNAs that were differentially expressed between the three LN treatment libraries were divided into four groups. In group A-I, the expression of the identified miRNA showed an increasing trend with grain development; there were 17 miRNAs in this group, including two miR156s and three miR167s, of which tae-miR9674b-5p (Table S4–1) and tae-miR9622a-3p were the most abundant (Table S4–1). The tae-miR9666b-3p miRNA showed an increasing degree of expression in developing grain, with 3, 213, and 674 normalized reads in the LN-7, LN-17, and LN-27 libraries, respectively. This expression profile might indicate that these miRNAs have more important regulatory roles at later stages of grain development. Group A-II miRNAs show a decreasing expression profile with grain development, and there were 30 miRNAs in this group, including two miR159s, three miR171s, and four miR169s. In this group, ata-miR5168-3p was the most abundant, and ata-miR169d-3p_L-2R + 2 showed large fold-changes of − 3.81, − 1.18, and − 4.99 for the LN-17 vs, LN-7, LN-27 vs. LN-7, and LN-27 vs. LN-7 comparisons, respectively. The high levels of expression at the early stages of wheat grain development might be related to grain endosperm differentiation and cell division. The members of group A-III showed a single expression peak, with the highest abundance at 17 DAA. In this group, there were 18 miRNAs, of which tae-miR9655-3p was the most abundant. In addition, there were 15 miRNAs in groups A-IV, including three miR159s, four miR166s, and two miR319s (Table S4–1).

Similarly, there were four miRNA groups (for the known miRNAs) for the HN treatment libraries (Table S4–2). Group B-I contained 15 miRNAs, including two miR156s, one miR408, and three miR9666s. There were 30 miRNAs in group B-II belonging to 17 miRNA families. The 24 miRNAs in group B-III included three miR160s, three miR167s, and two miR531s, of which tae-miR160 (with an increase from 4129 reads to 6157 reads) was the most abundant. There were 22 miRNAs in groups B-IV, including four miR159s and four miR166s. The fold-changes for the HN-17 vs. HN-7 and HN-27 vs. HN-17 comparisons were − 2.23 and − 2.19, respectively, for tae-miR159a (Table S4–2). We noticed that some miRNAs exhibited the same expression profiles during grain development in both the HN and LN treatments. For example, tae-miR408_L-1, osa-miR827, tae-miR9666b-3p, tae-miR9666a-3p, and tae-miR9664–3p_L-1 were in group I (A-I, and B-I) in both the HN and LN treatments. Moreover, miR168, miR396, and miR530 were in group II (A-II, and B-II), and miR-528, miR-1127, and miR9655 were all in group III (A-III, and B-III) (Table S4–1, S4–2). These miRNAs exhibiting the same expression patterns between the two N application levels might be closely related to grain development. Of course, some miRNAs showed different expression profiles with respect to N application levels, which indicates that the regulatory role of these miRNAs might be related to N supply status. For example, ppt-miR894_L-1R + 1 was differentially expressed only in the HN treatment, whereas tae-miR9652-5p showed differential expression only in the LN treatment.

In addition, the expression patterns of novel miRNAs during grain development for the LN treatment and the HN treatment are shown in Table S4–3 and Table S4–4, respectively. Similar to the known miRNAs, the novel miRNAs identified among the LN treatments were also divided into four groups (Table S4–3). Group C-I contains 24 novel miRNAs, of which tae-miR2916-p5_2ss2AG17TG was the most abundant, and 15 miRNAs were not expressed in the LN-7 library. There were 18 and 31 novel miRNAs in groups C-II and C-III, respectively. Only three novel miRNAs were in group C-IV, of which tae-miR9674a-p5 was the most abundant (Table S4–3). Similarly, the differentially expressed novel miRNAs in the HN treatment libraries were also classified into four groups (Table S4–4). In group D-I, there were 30 novel miRNAs, of which 21 showed no expression in the HN-7 library but were highly expressed in the HN-27 library. There were 12 miRNAs in group D-II, of which tae-miR171a-p5 exhibited a large expression fold change (− 4.87) from 7 DAA to 17 DAA. There were 40 and 6 novel miRNAs found in groups D-III and D-IV, respectively (Table S4–4).

Target gene functions of miRNAs differentially expressed during grain development

Based on Gene Ontology (GO) assignments, we identified 3011 potential target genes for 354 differentially expressed miRNAs. The target functions of the miRNAs were analyzed according to the “biological process” (A), “cellular component” (B), and “molecular function” (C) GO categories. In the “biological process” category, the target genes were associated largely with 10 GO terms (Fig. 3), with the three most represented being “growth and development” (24.61%), “regulation of transcription” (22.55%), and “cellular process” (14.11%). Twelve molecular function GO terms were identified, with the three most frequent being “protein binding” (28.28%), “ATP binding” (18.32%), and “DNA/RNA binding” (11.68%). For “cellular component”, the most abundant GO term was “nucleus” (59.77%), followed by “mitochondrion”, with 19.01%.

Fig. 3
figure 3

The major GO categories of “cellular component”, “biological process”, and “molecular function” for the predicted genes of all the differentially expressed miRNAs during grain development

Full size image

To explore the roles of these differentially expressed miRNAs in grain development, we listed the known miRNAs and the novel miRNAs and their target functions in Tables 3 and 4, respectively (target gene IDs in Additional file 5: Table S5). It is clear that the target genes have diverse functions, and the encoded gene products include transcription factors, involve Skp-cullin-F-box (SCF) ubiquitin ligase activity, encode Cytokinin Response Factor 1 (CRF), encode protein kinases and encode heat-shock proteins (HSPs). tae-miR9674b-5p, tae-miR160, ata-miR166a-3p, and tae-miR159a are relatively abundant during all gain development stages (RAGS), and their target genes have multiple functions and encode proteins such as SCF proteins, protein kinase 20-like, serine carboxypeptidase, kinesin-2-like protein pollen semisterility 1 (PSS1) and MYB transcription factors (Table 3), indicating that these miRNAs might have multiple and important roles during all stages of grain development. The miRNAs osa-miR171a, ata-miR172c-3p, three miR396s (ata-miR396a-5p, ata-miR396e-5p, and ata-miR396c-5p), ata-miR393-5p_L-1R + 1 and ata-miR5168-3p were relatively more abundant during early stages (RAES) (HN-7 and LN-7) compared to later stages of development, and the gene products of their predicted targets include the DELLA proteins RGL2 (REPRESSOR OF GA1–3-LIKE 2), CRF1, proteins involved in auxin-activated signaling, serine-type endopeptidase, proteins that respond to gibberellin (GA) and leucine-rich-repeat (LRR) kinase. These miRNAs might play crucial regulatory roles in seed formation. Additionally, several miRNAs, including ata-miR528-5p, bdi-miR531_L-4R + 1_1ss5CT, tae-miR1127b-3p_1ss12TC, tae-miR9652-5p and tae-miR9655-3p, were relatively highly expressed during the middle stages of development (RAMS) (LN-17 and HN-17). The predicted targets of these miRNAs are LRR receptor kinases (involved in starch and sucrose metabolism), NAD(P)H dehydrogenase, UDP-glucuronate decarboxylase, proteins involved in transport, and L-fucosidase, which indicates that they might regulate starch accumulation and processes related to grain filling. The relative abundance of some miRNAs, including two miR156 family members (bdi-miR156h-3p_L + 1 and zma-miR156d-3p_L + 1_1ss9TC), tae-miR167b-5p, tae-miR319, osa-miR827, and tae-miR408_L-1, showed increased abundance at the later stage (RALS); their predicted targets were WRKY transcription factors, SCF proteins, MYB transcription factors, L-ascorbate oxidase, NB-ARC, and SYG1/PHO81/XPR1 (SPX), which suggests that these miRNAs might play important regulatory roles during later stages of grain development.

Table 3 Conserved miRNAs differentially expressed in developing wheat grain and their predicted target gene functions
Full size table
Table 4 Differentially expressed novel miRNAs associated with grain development and their target gene functions
Full size table

For the novel miRNAs (Table 4), two miR2916 members, tae-miR2916-p5_2ss2AG17TG and tae-miR2916-p3_2ss2TC17CA, not only presented the highest abundance in the grain but also increased in abundance during grain filling. These two miRNAs are predicted to target genes encoding serine/threonine kinases, glutathione S-transferase (GSTF) 1, peroxidase, fructose-bisphosphate aldolase, DELLA proteins, and proteins with GTPase activity, indicating that they may have multiple functions during grain development and could play important roles in regulating grain weight and nutrient accumulation. Many miRNAs, including tae-miR5384-p5, PC-5p-3645_2157, PC-3p-4780_1750, PC-5p-24289_247, PC-5p-18073_397, and PC-3p-21502_303, exhibited relatively high expression levels in the middle stage of wheat grain development and were predicted to target genes encoding UDP-glycosyltransferase, metalloendopeptidase, neural precursor cell-expressed developmentally downregulated 8 (NEDD8)-specific protease, and glucosidase.

Cell differentiation plays an important role in determining wheat grain size. To clarify the potential regulatory mechanisms underlying cell differentiation, we constructed a miRNA target gene GO network by selecting miRNAs and their target genes based on their GO annotation (Fig. 4a).

Fig. 4
figure 4

Relationships within a miRNA-gene-GO network for cell differentiation, seed development and response to heat stress as determined by Cytoscape. The red, green, and yellow colors indicate miRNAs, target gene functions and GO terms, respectively

Full size image

There are 27 miRNAs and 20 genes in this network, with five, five, and three miRNAs belonging to the miR171, miR166, and miR319 families, respectively. It is clear that miR166 members target mainly transcription factor-encoding and nucleic acid-related genes, and the miR171 family members target genes mainly involved in the hormone response. Additionally, we constructed a network with miRNAs that had GO annotations related to seed development (Fig. 4b). Fifteen miRNAs and 16 genes are included in this network, including ata-miR167b-3p, which is predicted to target glycine cleavage system H protein 2. The tae-miR815c-p5_2ss6CT21CT miRNA might regulate grain development by targeting cytochrome P450s, which play roles in carbohydrate metabolism. Additionally, a network including eight miRNAs and eight target genes was constructed with GO annotations related to the heat stress response (Fig. 4c). This network shows that the miRNA targets include genes that encode chaperones, transcription factors, HSPs, DnaJ, and components of the abscisic acid (ABA) pathway, indicating that the miRNAs might participate in the heat stress response by regulating various pathways. Additionally, miR172 members might play an important role in the response to heat stress by regulating chaperones.

miRNAs that are differentially expressed between the two N levels and their potential target gene functions

A comparison between HN and LN revealed 40 (25 known and 15 novel), 44 (26 known and 18 novel), and 100 (48 known and 52 novel) differentially expressed miRNAs at 7 DAA, 17 DAA, and 27 DAA, respectively (Additional file 6: Table S6–1, S6–2, and S6–3). To further explore the role of these miRNAs in grain development and the response to N supply levels, the criteria of log2(fold change) ≥1.0 and reads per million of total miRNA reads (RPM) ≥20 were further used to select the differentially expressed miRNAs in the constructed libraries. For the LN-7 vs. HN-7, LN-17 vs. HN-17, and LN-27 vs. HN-27 comparisons, the ratios of upregulated miRNAs were 2/6, 0/6, and 5/10 for the known miRNAs, respectively, and were 5/5, 2/3, and 8/24 for the novel miRNAs, respectively (Table 5).

Table 5 Number of miRNAs that show differential expression between the HN and LN treatments
Full size table

The differentially expressed miRNAs and their target gene functions are listed in Table 6. For the known miRNAs, the abundance of sbi-miR399a was significantly higher in the HN treatments (HN-7 and HN-17) than in the LN treatments (LN-7 and LN-17), and tae-miR9664–3p_L-1 showed significantly higher expression in the HN treatment than in the LN treatment at 17 and 27 DAA. Additionally, ata-miR172b-5p had a higher expression level in the LN-27 treatment than in the HN-27 treatment. For the novel miRNAs, we found more differentially expressed miRNAs at the late developmental stages (LN-27 and HN-27). PC-5p-26413_215 had significantly higher expression in the LN treatment than in the HN treatment at the three grain developmental stages and potentially regulates BRASSINOSTEROID INSENSITIVE 1-associated receptor kinase 1 (BAK1) and LRR receptor-like serine/threonine-protein kinase. We found that 16 novel miRNAs were significantly upregulated in the HN-27 library compared to the LN-27 library, and these include two tae-miR2916 family members (tae-miR2916-p3_2ss2TC17CA and tae-miR2916-p5_2ss2AG17TG) that show fold-changes of 1.74–1.98 between LN-27 and HN-27. Both PC-3p-15988_474 and PC-5p-6533_1330 were more abundant in LN-27 than in HN-27 and potentially function by targeting defense (peroxidase activity)- and transcription factor-encoding genes. PC-5p-36040_123 also showed higher expression in LN-27 compared to HN-27, and its predicted target gene product is disulfide isomerase, which indicates that it might be related to the regulation of grain storage and protein quality.

Table 6 Differentially expressed miRNAs between HN and LN treatments and their predicted target gene functions
Full size table

A network of the miRNAs that are differentially expressed between the HN and LN treatments was constructed on the Cytoscape platform v.3.6.0 (https:cytoscape.org) [36] (Additional file 7: Figure S1). The differentially expressed miRNAs have multiple potential functions, such as targeting carbohydrate metabolism, nitrogen metabolism, signal transduction, hormone responses, and defense. However, it is clear that eight miRNAs, tae-miR2916-p5_2ss2AG17TG, PC-3p-4066_99, tae-miR5386-p5, PC-5p-26413_215, ata-miR164a-5p, tae-miR9655-p5, PC-3p-21502_303, and sbi-miR399a, have the same target—genes encoding transcription factors. Apart from carbohydrate metabolism and signal transduction, five miRNAs, tae-miR2916-p3_2ss2TC17CA, PC-3p-15988_474, osa-miR530-5p_L + 1, tae-miR2592bj-p3_2ss12TC19AT, and ata-miR528-5p, were predicted to target genes that function in transport. We also noticed that tae-miR9664–3p_L-1 not only has multiple potential functions but also is linked to the transcription factor group and the transport group via the RPP13 and RPM1 proteins. The results of these analyses indicate that the differentially expressed miRNAs might respond to different nitrogen fertilizer conditions by regulating multiple metabolic pathways. Additionally, transcription factors and transport proteins may play important roles in the response to N supply status. We also constructed a model to show how the regulation of differentially expressed miRNAs is related to grain development and the response to N supply (Fig. 5).

Fig. 5
figure 5

Model for miRNAs and their target genes associated with grain development and the response to N levels. The black letters in green circles indicate miRNAs, and the red letters in black boxes indicate the target gene function. The black arrows represent the direction of regulation. TF, transcription factor; CRF, cytokinin response factor; SCF, Skp-cullin-F-box; GLU, glucosyltransferase; BRI, brassinosteroid insensitive; ASR, alternative splicing regulator; LRR-RPK, leucine-rich-repeat receptor-like protein kinase; K1P-PSS1, kinesin-1-like protein PSS1; ARF, auxin response factor; HNT, high-affinity nitrate transport; Redox, reduction-oxidation; HSP, heat-shock protein; POD, peroxidase

Full size image

Validation of identified miRNAs and their target genes

To confirm the miRNA expression in this study and to validate the deep sequencing results, we randomly selected four and five miRNAs (including one novel) from the LN and HN libraries, respectively, for quantitative real-time polymerase chain reaction (qRT-PCR) assays. The expression patterns of these miRNAs were similar to those obtained from deep sequencing, suggesting that the sRNA sequencing results obtained in this study are reliable (Fig. 6).

Fig. 6
figure 6

Verification of the expression patterns of nine miRNAs present in developing wheat grain. The different lowercase letters above the columns indicate significant differences (P < 0.05)

Full size image

To verify the expression patterns of the potential miRNA targets, four target genes (the gene IDs and the primer sequences are given in Additional file 8: Table S7) were selected for gene-specific qRT-PCR under the HN and LN treatments (Fig. 7). The expression levels of PC-5p-3645_2157, tae-miR9655-p5, tae-miR9655-3p, and tae-miR319 were found to be upregulated with grain development, and their predicted target genes showed downregulated expression patterns. These results indicate that a negative relationship exists between target gene expression and their corresponding miRNA expression profiles.

Fig. 7
figure 7

Verification of the expression patterns of miRNA target genes in developing wheat grain. The different lowercase letters above the columns indicate significant differences (P < 0.05)

Full size image

Discussion

Grain characteristics under different N treatments

Nitrogen plays an important role in plant growth and agricultural production [24, 26], and N fertilization is an effective agronomic practice to increase wheat grain yield and improve grain processing quality [27, 28]. Here, HN treatment resulted in higher grain yields, protein content and dough development time, which is consistent with the results of previous reports [27, 28, 37]. We also observed that the gluten (gliadin and glutenin) content (27.80 mg g−1) increased more than did the albumin and globulin content (9.74 mg g−1) under HN treatment, indicating different responses of the grain protein component to N fertilizer. According to the national standards of high-quality wheat [38], the grain under HN conditions had a high protein content, and its dough rheological characteristics (water absorption, dough development time and stability) met the standard for strong-gluten wheat, which confirms that N fertilizer applications can improve grain processing quality.

miRNAs involved in wheat grain development

A better understanding of the regulatory roles of miRNAs in grain development and the response to nitrogen application levels will provide new perspectives for improving grain yield and quality. Previous studies have identified several conserved miRNAs in wheat grain, such as miR159, miR160, miR164, and miR156 [12, 15]. In this study, we also identified several highly conserved miRNA families, such as miR156, miR396, miR160, miR172, and miR164, indicating that these miRNAs are abundant in wheat grain and potentially play important regulatory roles in grain development. The reproductive stage of wheat from flowering to maturity includes the grain formation and filling stages; the number of endosperm cells in wheat is directly related to grain yield and quality [39]. A rapid increase in the number of endosperm cells occurs at 6–9 DAA, and the final number is determined by 14–18 DAA [39]. Members of the osa-miR171a, ata-miR172c-3p, miR396 (ata-miR396e-5p, ata-miR396a-5p, and ata-miR396c-5p) and ata-miR5168-3p families were the most abundant miRNAs at 7 DAA, suggesting that these miRNAs might be involved in the early stages of grain development. Higher abundances of miR171a and miR172 members were also found in the early developing grain by Meng et al. [12]. Additionally, by analyzing the miRNA target gene GO terms, we identified 27 miRNAs with target gene products related to cell differentiation (Fig. 4), including osa-miR171a and ata-miR172c-3p. It has been reported that miR172 members control floral organ identity in rice, maize, and barley by regulating the expression of APETALA2 (AP2)-like genes [40,41,42]. The floral and seed development of transgenic rice overexpressing miR172b were defective; the plants also presented relatively low spikelet fertility and reduced seed weight [42]. Here, the predicted target of ata-miR172c-3p is the CRF gene. CRFs compose a small subset of APETALA2/ethylene response factor (AP2/ERF) transcription factors that regulate plant development as part of the cytokinin signal transduction pathway [43], suggesting that ata-miR172c-3p might be involved in wheat seed development via cytokinin regulation. A positive correlation between cytokinin level during grain development and final grain weight has been reported in maize and wheat [44, 45]. Apart from cell differentiation, we also found that the target gene products of ata-miR172c-3p (DnaJ) and ata-miR172b-3p (chaperone binding) are involved in seed development and the response to heat stress (Fig. 4), which might indicate that they play important regulatory roles in determining wheat grain size and weight. Ma et al. [46] found that the miR171-SCL module plays a key role in mediating gibberellin (GA)-DELLA signaling in the coordinated regulation of leaf growth and chlorophyll biosynthesis in the light. We also predict that potential targets for osa-miR171a are genes that encode proteins such as the DELLA protein RG2-like, which is a regulator in the GA signaling pathway. It has been shown that miR396 members affect rice yield by regulating grain size [47], and the base substitution from TT to AA of the target genes can counteract the suppression by miR396, resulting in the activation of the hormone response and an increase in grain size [48, 49]. Schommer et al. [50] also found that the transcription factor TCP4, which is a direct regulator of miR396b, can repress cell proliferation by activating various pathways. Here, three miR396s (ata-miR396e-5p, ata-miR396a-5p, ata-miR396c-5p) were relatively abundant at 7 DAA, and the encoded proteins of the target genes function in response to GA, function as transcription factors, are GTPases, and act as protein kinases, suggesting that the regulation of these miRNAs might be related to grain size. All these results suggest that these miRNAs may play crucial roles in regulating grain endosperm cell development, which is directly related to grain size. However, additional research is needed to explore the target genes and their functions.

Grain filling is an important process in kernel development and ultimately determines kernel weight and flour quality. It has been suggested that 14–18 DAA is a major transition point for the formation of the wheat starchy endosperm and essentially marks the end of endosperm cell division and the deposition of starch and gluten protein [39]. Meng et al. [12] identified some miRNAs that might be involved in the control of grain filling; these miRNAs are involved in various biological processes by targeting many different wheat genes, such as those involved in carbohydrate metabolism, protein metabolism, transport, and transcription. In our study, we identified several conserved miRNAs and novel miRNAs that not only were highly abundant in wheat grain but also were significantly up- or downregulated during grain development; examples include tae-miR160 and ata-miR166a-3p. The gene products of the potential gene targets of these miRNAs include transcription factors, auxin response factors (ARFs), serine carboxypeptidases, kinesin-1-like protein PSS1, and proteins involved in transport. It is well known that miR160 directs the cleavage of several ARF mRNAs to regulate plant growth and seed germination [51, 52]. Expression of a miR160-resistant version of ARF17 in plants causes dramatic developmental defects, and overexpression of miR160 results in fewer starch granules and disorganized root tips [52, 53]. Analysis of the conservation and diversification of the miR166 family in plants showed that miR166 members play a wide and important regulatory role in seed development [54]. Therefore, these data suggest that these miRNAs have multiple regulatory functions in grain filling. In this study, we also identified some significantly upregulated novel miRNAs (such as ae-miR5384-p5, PC-5p-3645_2157, PC-3p-4780_1750, PC-5p-24289_247, PC-5p-18073_397, and PC-3p-21502_303) in the middle stage of grain development (17 DAA). The predicted proteins encoded by the targeted genes include UDP-glycosyltransferase, NEDD8-specific protease, metalloendopeptidase, and glucosidase, which indicate that these upregulated miRNAs might be related to carbohydrate and nitrogen metabolism, which would affect starch deposition and protein accumulation.

Some miRNAs, such as tae-miR319, hvu-miR399, and miR9666, showed little or no expression during the early stage (7 DAA) but then became upregulated as grain development progressed, and these miRNAs were most abundant at 27 DAA. The potential encoded proteins of the genes targeted by these miRNAs include MYBs, DnaJ, NB-ARC domain-containing proteins, and alternative splicing regulators. miR319 is a conserved microRNA that regulates TCP transcription factors involved in various developmental pathways, such as hormone biosynthesis, signaling, and leaf development and senescence [55, 56]. We also predicted that one of the potential targets of miR319 in the response to heat stress is the gene encoding the DnaJ protein. High-temperature stress at grain filling is a major problem in the production of high-quality wheat, and further exploration of the internal mechanism of these miRNAs in response to heat stress would be a potential way to reduce the effects of high-temperature stress.

Grain quality is a complex trait and includes not only protein content but also wheat gluten polymer structure [57]. Chen et al. [58] reported that cysteine and methionine metabolism pathways have an important role in storage protein biosynthesis via the synthesis of sulfur-containing amino acids. Here, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis showed that the target genes of bdi-miR531_L-4R + 1_1ss5CT and ata-miR396e-5p were annotated as involved in cysteine and methionine metabolism. The formation of disulfide (SS) bonds regulated by protein disulfide isomerase is a key factor that determines gluten protein polymerization [57, 59]. We identified one novel miRNA, PC-5p-36040_123, which is predicted to target protein disulfide isomerase, that was significantly downregulated at HN-27. The downregulated expression under high nitrogen treatment is predicted to cause an increase in disulfide isomerase activity and to promote grain protein polymer formation. In addition to SS bonds, hydrogen or electrostatic bonds or other noncovalent bonds also participate in the formation of gluten protein polymers [57]. Additional information concerning the regulatory function of the identified miRNAs and their relationship to the protein polymer structures in grain and flour is needed.

Differentially expressed grain miRNAs in response to nitrogen application levels

Many miRNAs have been identified that respond to nutrients and abiotic stress, such as salt and drought stress [60]. In this study, we identified 50 miRNAs (including 28 novel ones) that showed significant differential expression between HN and LN treatments. The expression of miR399 is upregulated by phosphorus (P) deficiency [61] and is related to P homeostasis by regulating the expression level of the target PHO2 gene, which encodes a ubiquitin-conjugating E2 enzyme, UBC24 [62]. Baek et al. [63] also found that AtMYB2, a transcription factor, can regulate miR399 expression levels in response to P starvation. Similarly, Zhao et al. [23] reported that TamiR399 was downregulated under conditions of N deficiency. In our study, we found that the expression of two miR399s (sbi-miR399a and hvu-miR399) was significantly downregulated in the LN treatment. The two miRNAs are involved in regulating target genes encoding transcription factors, oxidoreductase (defense response), and UDP-glucosyltransferase, suggesting that these miRNAs possibly mediate signal transduction in response to low N by repressing translation.

Arabidopsis miR169 has been shown to be involved in adaptation to LN stress [64]. Zhao et al. [19] also reported that two novel miR169 species (miRC10 and miRC68) critically regulate the adaptation to LN in maize seedlings via their interaction with the mRNA that targets a nuclear transcription factor subunit. The strong downregulation of miR169 in Arabidopsis under LN affects the expression of its targets, NFYA (Nuclear Factor Y, subunit A) family members involved in the ABA-dependent pathway [65]. We found that two miR169s (ata-miR169d-3p_L-2R + 2 and ata-miR169d-5p) had significantly lower expression in LN-27 than in HN-27, and the predicted target gene functions included kinase activity and E3 ubiquitin-protein ligase activity. It was previously reported that Arabidopsis MIEL1 E3 ligase negatively regulates ABA signaling sensitivity by promoting MYB96 turnover [66]. Thus, the regulation of the grain response to N supply levels by miR169 might involve the ABA signaling pathway.

Previous studies have shown that miR827 responds strongly to P starvation by targeting genes encoding SPX-MFS proteins predicted to be implicated in P sensing or transport [67]. There is crosstalk between nitrate and P sensing in Arabidopsis [68], and overexpression of miR827, which targets nitrogen limitation adaption (NLA) mRNA, in Arabidopsis also resulted in similar physiological responses in both NLA and PHO2 (PHOSPHATE2-ubiquitin conjugase) mutants grown under either low nitrate or high P [68]. Fertilizer experiments also confirmed the interaction between N and P fertilizers, and the combination of N and P improved wheat grain yield and quality [69]. Here, we identified one osa-miR827, which is predicted to target a gene encoding an SPX-domain-containing protein that showed significant downregulation in LN-17. Liu et al. [70] found that NLA responded to the remobilization of nitrate between sources and sinks in Arabidopsis by regulating the degradation of NRT1.7, which depended on miR827. This finding suggests that miR827 might influence N accumulation and remobilization by regulating NLA. Additionally, the mode of action occurs via a posttranslational regulatory pathway, not at the level of transcription [71]. For all the identified differentially expressed miRNAs in our study, only one, PC-5p-26413_215, was significantly upregulated in the LN treatment at all three stages. The target gene of PC-5p-26413_215 is predicted to encode BAK1. It is well known that brassinosteroids (BRs) control various aspects of plant development and growth and that BR perception and signal transduction require active BRASSINOSTEROID INSENSITIVE 1 (BRI) and BAK1 [71]. However, additional work is needed to confirm the function of the target genes and the regulatory mechanism of this miRNA.

Nitrogen is an essential inorganic nutrient for wheat growth and plays many important roles in grain development. In this study, we found that 21 novel miRNAs were significantly up- or downregulated only in the LN-27 treatment (examples are PC-5p-6533_1330, PC-3p-15988_474, tae-miR9779-p3, PC-3p-11233_746, and tae-miR5384-p5) in comparison with HN-27, and they are predicted to target multiple genes involved in defense (peroxidase), encoding transcription factors (MYB), involved in nitrogen metabolism (metalloendopeptidase), involved in carbohydrate metabolism (UDP-glycosyltransferase), and involved in transport (hexose transport). Different expression patterns between the LN and HN treatments have also been reported by Zhao et al. [23]. Thus, the differentially expressed miRNAs and their target genes suggest that there are multiple N-responsive modules in wheat grain and that miRNAs may mediate the plant responses and adaptation to N deprivation via posttranscriptional or posttranslational regulation.

Conclusions

We identified 79 miRNAs (46 known and 33 novel miRNAs) that showed significant differential expression during wheat grain development under both HN and LN treatments. We also identified 50 miRNAs (22 known and 28 novel miRNAs) that responded to the level of applied nitrogen. The potential targets of these miRNAs are involved in multiple biological processes, including cell differentiation, carbohydrate metabolism, transcription, signal transduction, transport, and defense. Our results indicate that posttranscriptional regulation by miRNA-mediated networks can play a crucial role in grain development and the N response, which determine wheat grain weight and quality. Our work provides useful information for further exploration of the mechanisms by which miRNAs function to improve grain yield and quality.

Methods

Experimental design and sample preparation

Seeds of the winter wheat cultivar ‘ZM119’, provided by Henan Academy of Agricultural Sciences, were grown during the 2016–2017 growing season at the Scientific and Educational Park of Henan Agricultural University, Yuanyang, China (35°06′ N, 113°56′ E). The cultivar ‘ZM119’ was bred by Henan Academy of Agricultural Sciences, and approved by Henan Province Crop Variety Approval Committee in 2014. The organic matter content in the topsoil was 10.1 g kg− 1, and the available phosphorus and available potassium contents were 35.0 mg kg− 1 and 110.6 mg kg− 1, respectively. Two N fertilizer treatments with three replicates each were applied: LN (0 kg ha− 1) and HN (210 kg ha− 1). The area of each plot was 21 m2 (3 m × 7 m), and each plot received 315 g of P (P2O5) and 315 g of K (K2O) before sowing. Fifty percent of the total N (in the form of urea) was applied before sowing, and 50% was topdressed at the elongation stage. Seeds were planted on 14 October 2016 at a density of 200 seeds m−2. At flowering, wheat spikes having a similar size and flowering on the same day were marked. The wheat spikes were sampled at 7, 17, and 27 DAA. The sampled grain was immediately frozen in liquid nitrogen and then stored in an ultra-low-temperature refrigerator until use. Afterward, mature grain was ground in a Cyclotec Sample Mill (Foss Tecator AB, Höganäs, Sweden) after cleaning and removing the foreign grain.

Small RNA library construction and RNA sequencing

To construct eighteen sRNA libraries, six independent samples, HN-7, HN-17, and HN-27 as well as LN-7, LN-17, and LN-27 (grain sampled at 7, 17, and 27 DAA under HN and LN conditions, respectively) with three biological replicates were separately extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA quantity and quality were assessed using a Bioanalyzer 2100 and an RNA 6000 Nano LabChip Kit (Agilent Technologies, Santa Clara, CA, USA). RNA preparations with a RIN (RNA integrity number) > 7.0 were used for sRNA library construction. Small RNAs were first ligated to 5′ RNA/DNA chimeric oligonucleotide adaptors (Illumina) and then 3′ adaptor to generate ligation products. cDNAs were then produced by reverse transcribing the ligation products. The cDNAs were sequenced using an Illumina HiSeq 2500 instrument (Illumina, San Diego, CA, USA) according to the manufacturer’s recommended protocol.

Bioinformatic analysis and miRNA identification

The raw sequencing reads in the miRNA libraries were subjected to an Illumina pipeline filter (Solexa 0.3), and then an in-house program, ACGT101-miR (LC Sciences, Houston, TX, United States), was used to process the dataset for removing adaptor dimers, common RNA family members, repeats and low-complexity reads. Unique sRNA sequences between 18 and 25 nucleotides in length were subsequently mapped to specific species precursors in miRBase 21.0 (http://www.mirbase.org) using BLAST searches to identify conserved miRNAs and novel 5p- and 3p-derived miRNAs [72]. Internal mismatches in the sequence and length variations at the 3′ and 5′ ends were allowed in the alignment process. The unique sequences that mapped to specific mature miRNA species in the hairpin arms were characterized as known miRNAs. A potential novel 5p- or 3p-derived miRNA was identified when unique sequences were mapped to the other arm of known specific species of precursor hairpins, opposite to the annotated mature miRNA arm. The remaining sequences were further matched to other selected species precursors in miRBase 21.0 by BLAST searches, and the mapped pre-miRNAs were further used as BLAST queries to identify their genomic locations by searching the genomes of specific species. The results of the above two searches were classified as known miRNAs. Finally, unmapped sequences were queried via BLAST against the specific genomes, and hairpin RNA structures containing the sequences from the flanking 120-nt sequences were predicted using RNAfold software (http://rna.tbi.univie.ac. at/cgi-bin/RNAfold.cgi). The key criteria for miRNA prediction are reported in the literature [73]. The secondary structures of precursors were also predicted using RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). The raw data have been submitted to the Sequence Read Archive (SRA) at the NCBI (https://www.ncbi.nlm.nih.gov) under accession number PRJNA563099.

Identification of differentially expressed miRNAs

The expression levels were normalized as RPM values. A t-test was applied to evaluate statistical significance. If the RPM ratios between different treatments were greater than 2 (fold change ≥ 2) and if the P-value was ≤0.05, then the miRNAs were defined as differentially expressed miRNAs.

miRNA target prediction and annotation

The genes targeted by miRNAs were predicted by identifying miRNA-binding sites using computational target prediction algorithms (Target Finder). Target searching was performed using the Triticum aestivum L. IWGSC library (http://plants.ensembl.org/index.html). Potential miRNA targets were functionally annotated by querying the target sequence via BLAST against the NCBI library (https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/). The most abundant miRNA targets were also annotated with GO terms (http://www.geneontology.org/) and KEGG pathways (http://www.genome.jp/kegg/).

miRNA and target mRNA quantification using qRT-PCR

Reverse transcription reactions were carried out using miRNA First-Strand cDNA Synthesis SuperMix (TransScript) according to the manufacturer’s instructions. A SYBR PrimeScript miRNA RT-PCR Kit was used to perform qRT-PCR in a fluorescence detection system (TianGen Biotech, Beijing, China) following the manufacturer’s instructions. All reactions were performed in triplicate for each sample. The 2-∆∆CT method was then used to calculate the relative expression levels of the miRNAs [74]. Target genes of the miRNAs that were differentially expressed during grain development in the HN and LN treatments were selected to validate their expression patterns in developing wheat grain via qRT-PCR. Fisher’s least significant difference test (LSD) was subsequently used to distinguish differences in relative expression levels between different development stages using SPSS (Statistical Program for Social Science) software; P < 0.05 was considered statistically significant. All primer sequences used and the target gene identification results are given in Table S7 (Additional file 8).

Determination of grain total protein content and the individual fraction contents

Four grain protein fractions, albumins, globulins, gliadins, and glutenins, were extracted using the sequential extraction methods of Liu et al. [75]. The total protein content and the individual fraction contents in each sample were analyzed by an automatic Kjeldahl analyzer (Kjeltec 2300, Foss Tecator, Sweden) following the method of ICC (1994) Standard Method 105/2 (International Association for Cereal Chemistry, Vienna, Austria). Differences between HN and LN were evaluated with t-test via SPSS; A value of P < 0.05 was considered to be statistically significant.

Availability of data and materials

All data generated or analyzed during this study are included in this published article (and its supplementary information files). The sequencing raw reads were submitted to the NCBI database with the accession number PRJNA563099 (https://www.ncbi.nlm.nih.gov).

Abbreviations

ABA:

Abscisic acid

ARF:

Auxin response factor

BAK1:

BRASSINOSTEROID INSENSITIVE 1-associated receptor kinase 1

BRI:

BRASSINOSTEROID INSENSITIVE 1

CRF:

Cytokinin response factor

DAA:

Days after anthesis

GA:

Gibberellin

GO:

Gene ontology

GSTF:

Glutathione S-transferase

HN:

High nitrogen treatment

HSP:

Heat-shock protein

LN:

Low nitrogen treatment

LRR:

Leucine-rich repeat

miRNA:

microRNA

N:

Nitrogen

NLA:

Nitrogen limitation adaption

qRT-PCR:

Quantitative reverse transcription-polymerase chain reaction

SCF:

Skp-cullin-F-box; SPXSYK1/PHO81/XPR1

sRNA:

Small RNA

References

  1. Braun HJ, Atlin G, Payne T. Muti-location testing as a tool to identify plant response to global climate change. In: Reynolds MP, editor. Climate change and crop production. Wallingford: CABI; 2010. p. 115–38.

    Chapter  Google Scholar 

  2. Foulkes MJ, Slafer G, Davies WJ, Berry PM, Sylvester-Bradley R, Martre P, et al. Raising yield potential of wheat. III. Optimizing partitioning to grain while maintaining lodging resistance. J Exp Bot. 2011;62:469–86.

    Article  CAS  PubMed  Google Scholar 

  3. Cakir C, Gillespie ME, Scofield SR. Rapid determination of gene function by vius-induced gene silencing in wheat and barley. Crop Sci. 2010;50:S77–84.

    Article  CAS  Google Scholar 

  4. Jinek M, Doudna JA. A three-dimensional view of the molecular machinery of RNA interference. Nature. 2008;457(7228):405–12.

    Article  CAS  Google Scholar 

  5. Xie Z, Allen E, Fahlgren N, Calamar A, Givan SA, Carrington JC. Expression of Arabidopsis MIRNA genes. Plant Physiol. 2005;138(4):2145–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Lee Y, Kim M, Han J, Yeom KH, Lee S, Baek SH, et al. MicroRNA genes are transcribed by RNA polymerase II. EMBO J. 2004;23(20):4051–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.

    Article  CAS  PubMed  Google Scholar 

  8. Mallory AC, Vaucheret H. Functions of microRNAs and related small RNAs in plants. Nat Genet. 2006;38:S31–6.

    Article  CAS  PubMed  Google Scholar 

  9. Voinnet O. Origin, biogenesis, and activity of plant microRNAs. Cell. 2009;136(4):669–87.

    Article  CAS  PubMed  Google Scholar 

  10. Bian H, Xie Y, Guo F, Han N, Ma S, Zeng Z, et al. Distinctive expression patterns and roles of the miRNA393/TIR1 homolog module in regulating flag leaf inclination and primary and crown root growth in rice (Oryza sativa). New Phytol. 2012;196:149–61.

    Article  CAS  PubMed  Google Scholar 

  11. Chen X. A microRNA as a translational repressor of APETALA2 in Arabidopsis flower development. Science. 2004;303:2022–5.

    Article  CAS  PubMed  Google Scholar 

  12. Meng FR, Liu H, Wang KT, Liu LL, Wang SH, Zhao YH, et al. Development associated microRNAs in grains of wheat (Triticum aestivum L.). BMC Plant Biol. 2013;13:140.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Xue LJ, Zhang JJ, Xue HW. Characterization and expression profiles of miRNAs in rice seeds. Nucleic Acids Res. 2009;37(3):916–30.

    Article  CAS  PubMed  Google Scholar 

  14. Wang YY, Shi CN, Yang TX, Zhao L, Chen JH, Zhang N, et al. High-throughput sequencing revealed that microRNAs were involved in the development of superior and inferior grains in bread wheat. Sci Rep. 2018;8:13854.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Han R, Jian C, Lv JY, Yan Y, Chi Q, Li ZJ, et al. Identification and characterization of microRNAs in the flag leaf and developing seed of wheat (Triticum aestivum L.). BMC Genomics. 2014;15:289.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  16. Sun FL, Guo GH, Du JK, Guo WW, Peng HR, Ni ZF, et al. Whole-genome discovery of miRNAs and their targets in wheat (Triticum aestivum L.). BMC Plant Biol. 2014;14:142.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Gifford ML, Dean A, Gutierrez RA, Coruzzi GM, Birnbaum KD. Cell-specific nitrogen responses mediate developmental plasticity. Proc Natl Acad Sci U S A. 2008;105(2):803–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Xu Z, Zhong S, Li X, Li W, Rothstein SJ, Zhang S, et al. Genome-wide identification of microRNAs in response to low nitrate availability in maize leaves and roots. PLoS One. 2011;6(11):e28009.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Zhao Y, Xu Z, Mo Q, Zou C, Li W, Xu Y, et al. Combined small RNA and degradome sequencing reveals novel miRNAs and their targets in response to low nitrate availability in maize. Ann Bot. 2013;112:633–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Liang G, He H, Yu DQ. Identification of nitrogen starvation-responsive MicroRNAs in Arabidopsis thaliana. PLoS One. 2012;7(11):e48951.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Gao S, Guo C, Zhang Y, Zhang F, Du X, Gu J, et al. Wheat microRNA member TaMIR444a is nitrogen deprivation-responsive and involves plant adaptation to the nitrogen-starvation stress. Plant Mol Biol Rep. 2016;34:931–46.

    Article  CAS  Google Scholar 

  22. Zuluaga D, De Paola D, Janni M, Curci PL, Sonnante G. Durum wheat miRNAs in response to nitrogen starvation at the grain filling stage. PLoS One. 2017;12(8):e0183253.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  23. Zhao Y, Guo L, Lu W, Li X, Chen H, Guo C, et al. Expression pattern analysis of microRNAs in root tissue of wheat (Triticum aestivum L.) under normal nitrogen and low nitrogen conditions. J. Plant Biochem. Biot. 2015;24:143–53.

    CAS  Google Scholar 

  24. Vidal EA, Araus V, Lu C, Parry G, Green PJ, Coruzzi GM, et al. Nitrate-responsive miR393/AFB3 regulatory module controls root system architecture in Arabidopsis thaliana. Proc Natl Acad Sci U S A. 2010;107(9):4477–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Gutierrez RA, Lejiay LV, Dean A, Chiaromonte F, Shasha DE, Coruzzi GM. Qualitative network models and genome-wide expression data define carbon/nitrogen-responsive molecular machines in Arabidopsis. Genome Biol. 2007;8:R7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  26. Zhu ZL, Chen DL. Nitrogen fertilizer use in China – contributions to food production, impacts on the environment and best management strategies. Nutr Cycl Agroecosys. 2002;63(2–3):117–27.

    Article  CAS  Google Scholar 

  27. Ma D, Guo T, Wang Z, Wang C, Zhu Y, Wang Y. Influence of nitrogen fertilizer application rate on winter wheat (Triticum aestivum L.) flour quality and Chinese noodle quality. J Sci Food Agric. 2009;89:1213–20.

    Article  CAS  Google Scholar 

  28. Zhang M, Ma D, Ma G, Wang C, Xie X, Kang G. Response of glutamine synthetase activity and gene expression to nitrogen levels in winter wheat cultivars with different grain protein content. J Cereal Sci. 2017;74:187–93.

    Article  CAS  Google Scholar 

  29. Ju XT, Xing GX, Chen XP, Zhang SL, Zhang LJ, Liu XJ, et al. Reducing environmental risk by improving N management in intensive Chinese agricultural systems. Proc Natl Acad Sci U S A. 2009;106(9):3041–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Guo JH, Liu XJ. ZhangY, Shen JL, Han WX, Zhang WF, et al. significant acidification in major Chinese croplands. Science. 2010;327(5968):1008–10.

    Article  CAS  PubMed  Google Scholar 

  31. Chu ZL, Chen JY, Xu HX, Dong ZD, Chen F, Cui DQ. Identification and comparative analysis of microRNA in wheat (Triticum aestivum L.) callus derived from mature and immature embryos during in vitro culture. Front Plant Sci. 2016;7:1302.

    PubMed  PubMed Central  Google Scholar 

  32. Ma X, Xin Z, Wang Z, Yang Q, Guo S, Guo X, et al. Identification and comparative analysis of differentially expressed miRNAs in leaves of two wheat (Triticum aestivum L.) genotypes during dehydration stress. BMC Plant Biol. 2015;15:21.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Xin M, Wang Y, Yao Y, Xie C, Peng H, Ni Z, et al. Diverse set of microRNAs are responsive to powdery mildew infection and heat stress in wheat (Triticum aestivum L.). BMC Plant Biol. 2010;10:123.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Han H, Wang Q, Wei L, Liang Y, Dai J, Xia G, et al. Small RNA and degradome sequencing used to elucidate the basis of tolerance to salinity and alkalinity in wheat. BMC Plant Biol. 2018;18:195.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Achakzai HK, Barozai MYK, Din M, Baloch IA, Achakza AKK. Identification and annotation of newly conserved microRNAs and their targets in wheat (Triticum aestivum L.). PloS ONE. 2018;13(7):20200033.

    Article  CAS  Google Scholar 

  36. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Guttieri MJ, McLean R, Stark JC, Souza E. Managing irrigation and nitrogen fertility of hard spring wheats for optimum bread and noodle quality. Crop Sci. 2005;45:2049–59.

    Article  Google Scholar 

  38. SAC and AQSIQ. Quality classification of wheat varieties: People’s Republic of China National Standard GB/T 17320. Stand. Admin. of the People’s Republic of China, and Gen. Admin. Qual. Supervision, Insp., Quar. of the People’s Republic of China, Beijing, 2013.

    Google Scholar 

  39. Wang WH, Guo WS, Fang MK, Feng CN, Zhu XK, Peng YX. Endosperm cell proliferating and grain filling dynamics in wheat. Acta Agron Sin. 2003;29(5):779–84.

    Google Scholar 

  40. Lauter N, Kampani A, Carlson S, Goebel M, Moose SP. microRNA172 down-regulates glossy15 to promote vegetative phase change in maize. Proc Natl Acad Sci USA. 2005;102(26):9412–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Nair SK, Wang N, Turuspekov Y, Pourkheirandish M, Sinsuwongwat S, Chen G, et al. Cleistogamous flowering in barley arises from the suppression of microRNA-guided HvAP2 mRNA cleavage. Proc Natl Acad Sci U S A. 2010;107(1):490–5.

    Article  CAS  PubMed  Google Scholar 

  42. Zhu QH, Upadhyaya NM, Gubler F, Helliwell CA. Over-expression of miR172 causes loss of spikelet determinacy and floral organ abnormalities in rice (Oryza sativa). BMC Plant Biol. 2009;9:149.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  43. Rashotte AM, Goertzen LR. The CRF domain defines cytokine response factor proteins in plants. BMC Plant Biol. 2010;10:74.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  44. Dietrich JT, Kaminek V, Belvins DG, Reinbett TM, Morris RD. Changes in cytokinins and oxidase activity in developing maize kernel and the effects of exogenous cytokinin on kernel development. Plant Physiol Biochem. 1995;33:327–36.

    CAS  Google Scholar 

  45. Gupta NK, Gupta S, Shukla DS, Deshmukh PS. Differential responses of BA injection on yield and specific grain growth in contrasting genotypes of wheat (Triticum aestivuvm L.). Plant Growth Regul. 2003;40:201–5.

    Article  CAS  Google Scholar 

  46. Ma ZX, Hu XP, Cai WJ, Huang WH, Zhou X, Luo Q, et al. Arabidopsis miR171-targeted scarecrow-like proteins bind to GT cis-elements and mediate gibberellin-regulated chlorophyll biosynthesis under light conditions. PLoS Genet. 2014;10(8):e1004519.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  47. Duan P, Ni S, Wang J, Zhang B, Xu R, Wang Y, et al. Regulation of OsGRF4 by OsmiR396 controls grain size and yield in rice. Nat Plants. 2016;2:15203.

    Article  CAS  Google Scholar 

  48. Hu J, Wang Y, Fang Y, Zeng L, Xu J, Yu H, et al. A rare allele of GS2 enhances grain size and grain yield in rice. Mol Plant. 2015;8:1455–65.

    Article  CAS  PubMed  Google Scholar 

  49. Li S, Gao F, Xie K, Zeng X, Cao Y, Zeng J, et al. The OsmiR396c-OsGRF4-OsGIF1 regulatory module determines grain size and yield in rice. Plant Biotechnol J. 2016;14:2134–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Schommer C, Debernardi JM, Bresso EG, Rodrguez RE, Palatnik JF. Repression of cell proliferation by miR319-regulated TCP4. Mol Plant. 2014;7(10):1533–44.

    Article  CAS  PubMed  Google Scholar 

  51. Liu PP, Montgomery TA, Fahlgren N, Kasschau KD, Nonogaki H, Carrington JC. Repression of AUXIN RESPONSE FACTOR10 by microRNA160 is critical for seed germination and post-germination storage. Plant J. 2007;52:133–46.

    Article  CAS  PubMed  Google Scholar 

  52. Mallory AC, Bartel DP, Bartel B. MicroRNA-directed regulation of Arabidopsis AUXIN RESPONSE FACTOR17 is essential for proper development and modulate expression of early auxin response genes. Plant Cell. 2005;17(5):1360–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Wang JW, Wang LJ, Mao YB, Cai WJ, Xue HW, Chen XY. Control of root cap formation by MicroRNA-targeted auxin response factors in Arabidopsis. Plant Cell. 2005;17:2204–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Li X, Xie X, Li J, Cui Y, Hou Y, Zhai L, et al. Conservation and diversification of the miR166 family in soybean and potential roles of newly identified miR166s. BMC Plant Biol. 2017;17(1):32.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  55. Palatnik JF, Allen E, Wu X, Schommer C, Schwab R, Carrington JC, et al. Control of leaf morphogenesis by microRNAs. Nature. 2003;425(6955):257–63.

    Article  CAS  PubMed  Google Scholar 

  56. Ori N, Cohen AR, Etzioni A, Brand A, Yanai O, Shleizer S, et al. Regulation of LANCEOLATE by miR319 is required for compound-leaf development in tomato. Nat Genet. 2007;39(6):787–91.

    Article  CAS  PubMed  Google Scholar 

  57. Johansson E, Malik AH, Hussain A, Rasheed F, Newson WR, Plivelic T, et al. Wheat gluten polymer structures: the impact of genotype, environment, and processing on their functionality in various applications. Cereal Chem. 2013;90:367–76.

    Article  CAS  Google Scholar 

  58. Chen X, Yang Y, Ran L, Dong Z, Zhang E, Yu X, et al. Novel insights into miRNA regulation of storage protein biosynthesis during wheat caryopsis development under drought stress. Front Plant Sci. 2017;8:1707.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Every D, Griffin WB, Wilson PE. Ascorbate oxidase, protein disulfide isomerase, ascorbic acid, dehydroascorbic acid and protein levels in developing wheat kernels and their relationship to protein disulfide bond formation. Cereal Chem. 2003;80:35–9.

    Article  CAS  Google Scholar 

  60. Jian HJ, Wang J, Wang TY, Ly W, Li JN, Liu LZ. Identification of rapeseed microRNAs involved in early stage seed germination under salt and drought stresse. Front Plant Sci. 2016;13:00658.

    Google Scholar 

  61. Fujii H, Chiou TJ, Lin SI, Aung K, Zhu JK. A miRNA involved in phosphate-starvation response in Arabidopsis. Curr Biol. 2005;15(22):2038–43.

    Article  CAS  PubMed  Google Scholar 

  62. Aung K, Lin SI, Wu CC, Huang YT, Su CL, Chiou TJ. pho2, a phosphate overaccumulator, is caused by a nonsense mutation in a microRNA399 target gene. Plant Physiol. 2006;141:1000–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Baek D, Kim MC, Chun HJ, Kang S, Park HC, Shin G, et al. Regulation of miR399f transcription by AtMYB2 affects phosphate starvation responses in Arabidopsis. Plant Physiol. 2013;161:362–73.

    Article  CAS  PubMed  Google Scholar 

  64. Zhao M, Tai H, Sun S, Zhang F, Xu Y, Li WX. Cloning and characterization of maize miRNAs involved in responses to nitrogen deficiency. PLoS One. 2012;7:e29669.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Zhao M, Ding H, Zhu JK, Zhang FS, Li WX. Involvement of miR169 in the nitrogen-starvation responses in Arabidopsis. New Phytol. 2011;190(4):906–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Lee HG, Seo PJ. The Arabidopsis MIEL1 E3 ligase negatively regulates ABA signalling by promoting protein turnover of MYB96. Nat Commun. 2016;7:12525.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Lin SI, Santi C, Jobet E, Lacut E, El Kholti N, Karlowski WM, et al. Complex regulation of two target genes encoding SPX-MFS proteins by rice miR827 in response to phosphate starvation. Plant Cell Physiol. 2010;51:2119–31.

    Article  CAS  PubMed  Google Scholar 

  68. Kant S, Peng M, Rothstein SJ. Genetic regulation by NLA and microRNA827 for maintaining nitrate-dependent phosphate homeostasis in Arabidopsis. PLoS Genet. 2011;7:e1002021.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Ma D, Zhang J, Hou J, Li Y, Huang X, Wang C, et al. Evaluation of yield, processing quality, and nutritional quality in different-colored wheat grains under nitrogen and phosphorus fertilizer application. Crop Sci. 2018;58:402–15.

    Article  CAS  Google Scholar 

  70. Liu W, Sun Q, Wang K, Du Q, Li WX. Nitrogen limitation adaptation (NLA) is involved in source-to-sink remobilization of nitrate by mediating the degradation of NRT1.7 in Arabidopsis. New Phytol. 2017;214(2):734–44.

    Article  CAS  PubMed  Google Scholar 

  71. Wang X, Goshe MB, Soderblom EJS, Phinney BS, Kuchar JA, Li J, et al. Identification and functional analysis of in vivo phosphorylation sites of the Arabidopsis BRASSINOSTEROID-INSENSITIVE1 receptor kinase. Plant Cell. 2005;17:1685–703.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Griffiths J, Saini HK, van Dongen S, Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36:D154–8.

    Article  CAS  Google Scholar 

  73. Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, Bowman JL, et al. Criteria for annotation of plant MicroRNAs. Plant Cell. 2008;20(12):3186–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Livak KJ, Schmitthen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2 −△△CT method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  75. Liu ZH, Cheng FM, Cheng WD, Zhang GP. Positional variations in phytic acid and protein content within a panicle of japonica rice. J Cereal Sci. 2005;41:297–303.

    Article  CAS  Google Scholar 

Download references

Acknowledgments

We would like to thank Prof. Yongchun Li for help with the bioinformatic analysis. We also thank the native English-speaking scientists from Elixigen Company (Huntington Beach, California) for proofreading our manuscript. Additionally, we appreciate the valuable comments from the reviewers.

Funding

The work was funded by the National Key Research and Development Program of China with grant 2016YFD0300404 and grant 31571651 from the National Natural Science Foundation of China.

Author information

Author notes

  1. Gege Hou and Chenyang Du contributed equally to this work.

Authors and Affiliations

  1. College of Agronomy/National Engineering Research Center for Wheat, Henan Agricultural University, Zhengzhou, 450002, China

    Gege Hou, Chenyang Du, Honghuan Gao, Sujun Liu, Wan Sun, Hongfang Lu, Juan Kang, Yingxin Xie, Dongyun Ma & Chenyang Wang

  2. The National Key Laboratory of Wheat and Maize Crop Science, Henan Agricultural University, Zhengzhou, 450002, China

    Dongyun Ma & Chenyang Wang

Authors
  1. Gege Hou
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chenyang Du
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Honghuan Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Sujun Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Wan Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Hongfang Lu
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Juan Kang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Yingxin Xie
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Dongyun Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Chenyang Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

CYW and YXX helped design the experiments; DYM designed the experiments and revised the manuscript. GGH carried out the experiments and wrote the draft. CYD and HHG performed the qPCR and prepared some of the figures. HFL and JK carried out the field trials and sampling. WS and SJL performed the grain quality analyses. All authors have contributed to the corrections and approval of the manuscript.

Corresponding authors

Correspondence to Dongyun Ma or Chenyang Wang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Illumina DNA sequencing data for the small RNA libraries.

Additional file 2: Table S2.

miRNAs identified in this study.

Additional file 3: Table S3.

Novel miRNAs identified in developing wheat grains at two different N levels.

Additional file 4: Table S4.

Conserved and novel miRNAs that are differentially expressed during wheat grain development.

Additional file 5: Table S5.

Target function of differentially expressed miRNAs during grain filling under HN and LN treatments.

Additional file 6: Table S6.

Wheat grain miRNAs that are differentially expressed between the HN and LN treatments.

Additional file 7: Figure S1.

Relationships between differentially-expressed miRNAs identified between the HN and LN treatments and their predicted targets as determined by Cytoscape.

Additional file 8: Table S7.

Oligonucleotide primers used for qRT-PCR assays in this study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hou, G., Du, C., Gao, H. et al. Identification of microRNAs in developing wheat grain that are potentially involved in regulating grain characteristics and the response to nitrogen levels. BMC Plant Biol 20, 87 (2020). https://0-doi-org.brum.beds.ac.uk/10.1186/s12870-020-2296-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12870-020-2296-7

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us