Skip to main content
Download PDF

Integrated analyses of miRNAome and transcriptome reveal zinc deficiency responses in rice seedlings

BMC Plant Biology volume 19, Article number: 585 (2019) Cite this article

Abstract

Background

Zinc (Zn) deficiency is one of the most widespread soil constraints affecting rice productivity, but the molecular mechanisms underlying the regulation of Zn deficiency response is still limited. Here, we aim to understand the molecular mechanisms of Zn deficiency response by integrating the analyses of the global miRNA and mRNA expression profiles under Zn deficiency and resupply in rice seedlings by integrating Illumina’s high-throughput small RNA sequencing and transcriptome sequencing.

Results

The transcriptome sequencing identified 360 genes that were differentially expressed in the shoots and roots of Zn-deficient rice seedlings, and 97 of them were recovered after Zn resupply. A total of 68 miRNAs were identified to be differentially expressed under Zn deficiency and/or Zn resupply. The integrated analyses of miRNAome and transcriptome data showed that 12 differentially expressed genes are the potential target genes of 10 Zn-responsive miRNAs such as miR171g-5p, miR397b-5p, miR398a-5p and miR528-5p. Some miRNA genes and differentially expressed genes were selected for validation by quantitative RT-PCR, and their expressions were similar to that of the sequencing results.

Conclusion

These results provide insights into miRNA-mediated regulatory pathways in Zn deficiency response, and provide candidate genes for genetic improvement of Zn deficiency tolerance in rice.

Background

Zinc is one of the essential micronutrients which acts as a catalytic, regulatory or structural co-factor for a lot of enzymes and regulatory proteins in plants and animals [1]. The well-known Zn-containing proteins in plants include the enzymes like carbonic anhydrase, alcohol dehydrogenase and copper (Cu)/Zn superoxide dismutase (CSD), and lots of Zn-finger domain-containing proteins, which function in transcriptional regulation [1, 2]. Zn deficiency (−Zn) is a serious agricultural problem in arable soils worldwide due to the low availability of Zn that is caused by factors such as high pH, prolonged flooding, low redox potential, and high contents of bicarbonate, organic matter and phosphorus [3]. Zn deficiency in soils and plants can also result in human malnutrition through the intake of food that contains low concentrations of Zn and other micronutrients [4]. It is estimated that one-third of the human population, especially children and women suffer from Zn-deficiency-related health problems [5]. Thus, studies to better understand the mechanisms of Zn deficiency response and tolerance in plants can help to develop crops with increased Zn use efficiency, Zn concentration, and Zn deficiency tolerance [4, 6,7,8].

Rice (Oryza sativa) is the second most cultivated crop and the main staple food for three billion people in the world. However, Zn deficiency is becoming one of the most widespread soil constraints affecting rice productivity [9, 10]. Much effects have been made to elucidate the physiological and biochemical mechanisms associated with Zn deficiency response and tolerance in plants [11,12,13]. These mechanisms include increasing Zn availability for root uptake by adjusting the root system architecture, releasing of phytosiderophores and organic acids, and formation of arbuscular mycorrhiza symbiosis [6, 14, 15]; expanding the generation of crown roots [16]; involvement of transporters such as ZRT, IRT-related proteins (ZIPs), metal tolerance proteins and heavy metal tolerance family proteins in Zn uptake and translocation [17,18,19,20,21]; adjusting the expression of Zn-requiring enzymes [2, 13]. Recently, transcriptomic profilings were analyzed to identify Zn-responsive genes in plants [22,23,24]. Four genes up-regulated in Zn-efficient rice varieties under Zn deficiency were reported to be candidate genes conferring Zn efficiency [22]. A large number of Zn deficiency responsive genes were found to be associated with calcium, sugar and hormonal signals in soybean [23]. Although some progress has been made, the molecular regulatory mechanisms underlying Zn deficiency response and tolerance is still limited.

microRNAs (miRNAs) are endogenous small noncoding RNAs which are 20 to 24 nucleotides in size and generated from a single-stranded RNA precursor with a hairpin secondary structure [25]. They are known to direct target mRNA cleavage, translational repression, and DNA methylation on the basis of sequence complementarity with their target genes. In plants, miRNAs mainly direct mRNA cleavage and have important functions in the regulation of growth and development and various environmental stress responses [26,27,28,29,30,31,32]. For example, miR156 is involved in controlling rice grain size, panicle branching, crown root development and flowering time by targeting SPL transcription factor genes [27, 33, 34]; miR164 regulates lateral root development, leaf senescence, grain yield and drought stress tolerance [35,36,37]. miRNAs also regulate plant responses to multiple nutrient deprivation stresses and are involved in nutrient homeostasis regulation [38,39,40]. For instance, miR399, miR827 and miR778 are induced by phosphate deficiency and are involved in regulating phosphate homeostasis by regulating ubiquitin-conjugating E2 enzyme PHO2, RING-type ubiquitin E3 ligase NLA, and histone H3 lysine 9 (H3K9) methyltransferase SUVH6, respectively [40,41,42,43]; miR397, miR398, miR408 and miR857 are involved in Cu homeostasis regulation by controlling the expression of Cu-containing proteins [39, 44]. By using miRNA microarray, eight miRNA families were identified to be responsive to Zn deficiency in Sorghum bicolor, and two Cu/Zn superoxide dismutase genes SbCSD1 and SbCSD2 were found be targeted by miR398 and miR528, respectively [2]. However, to our knowledge, there is still no report on whether and how miRNAs regulate Zn deficiency responses in rice. In this study, we investigated the molecular mechanisms in response to Zn deficiency by integrating miRNAome and transcriptome analyses in rice seedlings. A comprehensive and integrated analysis of these different datasets have identified potential miRNA-mRNA interactions under Zn deficiency in rice.

Results

Transcriptome profilings in response to Zn deficiency and Zn resupply by RNA sequencing

After experiencing Zn deficiency for 14 days, the Zn concentrations in shoots and roots were significantly decreased, and the Zn concentrations in shoots and roots of Zn-deficient seedlings were then significantly increased after Zn resupply for 3 days (Fig. 1a, b-d). No significant difference was observed in shoot and root biomass under Zn deficiency (Fig. 1e). Consistent with the phenotype observed in Sorghum [2], the primary root length were dramatically increased by Zn deficiency (Fig. 1c, f). These results showed that physiological adaptive responses to -Zn and Zn resupply have been exhibited.

Fig. 1
figure1

Physiological responses of rice seedlings to Zn deficiency and Zn resupply. a Schematic representation of the experimental design showing the duration of seedling growth and treatments of Zn deficiency and Zn resupply. Morphological appearance of the whole seedlings (b) and the roots (c) after Zn deficiency treatment for 14 days. d Zn concentration in the shoots and roots of rice seedlings after Zn deficiency for 14 days and Zn resupply for 3 days. Fresh weight of shoots and roots (e) and root length (f) of rice seedlings after Zn deficiency treatment for 14 days. Data are means ± SD from 3 samples for (d), and 10 samples for (e and f). Different letters indicate that the values are significantly different at P < 0.05 (Student’s t-test)

Full size image

We then analyzed the global transcriptome profiles of rice shoots and roots in responses to -Zn and Zn resupply by RNA sequencing. The Zn resupply treatments were used to investigate the recovery of Zn-responsive genes or miRNAs under Zn stress. By Illumina’s deep sequencing, a total of 28.6 to 55.2 million reliable clean reads were obtained from each library after excluding the low-quality reads, and most of the clean reads (86.8–97.3%) from each library could be mapped to the rice reference genome (https://rapdb.dna.affrc.go.jp/) (Additional file 1: Table S2). The Pearson’s correlation (R value) of the three biological replicates of each sample was around 90%, indicating the high reliability of the replicates (Additional file 2: Figure S1). The abundance of mapped transcripts was measured in terms of FPKM, and a total of 34,716 gene loci were detected in all these samples (Additional file 1: Table S3).

Differential expression analysis (fold change ≥2 and FDR ≤ 0.05) showed that a total of 151, 227, 123, and 205 genes were differentially expressed in Zn-deficient shoots (ZMS/ZPS), Zn-deficient roots (ZMR/ZPR), Zn-resupply shoots (ZRS/ZMS) and Zn-resupply roots (ZRR/ZMR), respectively (Fig. 2a; Additional file 1: Tables S4-S7). Venn diagram analysis showed that a total of 97 Zn-deficiency-responsive genes were recovered by Zn resupply in shoots and/or roots (Fig. 2b and c), suggesting a causal relationship between the expression levels of these genes and the levels of Zn in the growth media. Among the 89 up-regulated and 62 down-regulated genes in shoots (ZMS/ZPS up and down), 32 (36.0%) and 14 (22.6%) genes were recovered after Zn resupply treatment (ZRS/ZMS down and up), respectively (Table 1). Among the 117 up-regulated and 110 down-regulated genes in roots (ZMR/ZPR up and down), 40 (34.2%) and 18 (16.4%) genes were recovered after Zn resupply treatment (ZRR/ZMR down and up), respectively (Table 2). Seven genes (Os10g0328600, Os04g0280500, Os04g0561500, Os06g0566201, Os06g0566300, Os08g0207500, and Os07g0232800) were commonly recovered by Zn resupply in both roots and shoots; three of them encode Zn transporters (OsZIP4, OsZIP8 and OsZIP10). GO enrichment analysis (FDR < 0.05) revealed that zinc ion transmembrane transport (GO:0071577) and zinc ion transmembrane transporter activity (GO:0005385) were significantly enriched in Zn-deficiency up-regulated DEGs and Zn-resupply down-regulated DEGs in both shoots and roots, while threonine synthase activity (GO:0004795) and L-methionine biosynthetic process from methylthioadenosine (GO:0019509) were significantly enriched in Zn-deficiency down-regulated DEGs in shoots and roots, respectively (Additional file 1: Table S8).

Fig. 2
figure2

Overview of the differentially expressed genes (DEGs) in responses to 14 d of Zn deprivation and/or 3 d of Zn resupply. a The number of DEGs in responses to Zn deprivation and/or Zn resupply in rice roots and leaves (P < 0.05, FDR < 0.05). b Venn diagram representing the overlap of the Zn deficiency up-regulated DEGs and Zn resupply down-regulated DEGs in roots and leaves. c Venn diagram representing the overlap of the Zn deficiency down-regulated DEGs and Zn resupply up-regulated DEGs in roots and leaves

Full size image
Table 1 Zn deficiency-induced or repressed genes returning to basal level after Zn resupply in the shoots. The FPKM value represents mean ± SD (standard deviation) of three biological replicates. Note: MSTRG.18312, MSTRG.5686, and MSTRG.5688 are three new loci. “-” means no description
Full size table
Table 2 Zn deficiency-induced or repressed genes returning to basal level after Zn resupply in the roots. The FPKM value represents mean ± SD of three biological replicates. Note: MSTRG.6046, MSTRG.22477, and MSTRG.30070 are three new loci. “-” means no description
Full size table

miRNA profilings in rice responses to Zn deficiency and Zn resupply

Sequencing results showed that 8.8 to 15.3 million raw reads were obtained for each of the 18 libraries, and after discarding low-quality, junk and adaptor reads, repeats, and Rfam RNA and mRNA sequences, 38.3 to 56.0% of the total reads were valid reads in these libraries (Additional file 1: Table S9). The size distribution analysis showed that small RNAs from the six samples were primarily enriched in sequences with 21- to 24-nt, and 24-nt was the most dominant segment (Additional file 1: Table S10; Additional file 2: Figure S2), which is in concurrence with previous studies [45]. Pearson correlation analysis of three biological replicates of each sample revealed high correlation coefficients ranging from 0.80 to 0.99 (Additional file 2: Figure S3), indicating the high reliability of the replicates. The clean reads were searched against the miRNAs of rice and other plant species from miRBase, and a total of 498 known/conserved miRNAs were identified from all samples. The reads which could not be mapped to miRBase were subjected to novel miRNA prediction. After the removal of small RNAs which do not meet the plant miRNA criteria (eg. length < 20 or > 24, no hairpin structure, less than 10 reads in a sample), a final set of 370 novel miRNA genes was obtained (Additional file 1: Table S11).

A total of 68 miRNAs (including 38 novel miRNAs) were found to be differentially expressed under Zn deficiency (ZMS/ZPS in shoots and ZMR/ZPR in roots) and/or Zn resupply (ZRS/ZMS in shoots and ZRR/ZMR in roots); 38 miRNAs were differentially expressed under Zn deficiency and 44 miRNAs were differentially expressed under Zn resupply; 14 miRNAs were commonly responsive to Zn deficiency and Zn resupply (Additional file 1: Table S12; Fig. 3; Additional file 2: Figures S4 and S5). Of these Zn-responsive miRNAs, 12 and 11 were up-regulated by Zn deficiency in shoots and roots, respectively; 2 and 12 were down-regulated by Zn deficiency in shoots and roots, respectively (Additional file 2: Figure S4). However, none of these miRNAs were commonly responsive to Zn deficiency in both shoots and roots, suggesting the distinct miRNA-mediated regulatory networks in shoots and roots under Zn deficiency. Seventeen and 27 miRNAs were found to be responsive to Zn resupply in shoots and roots, respectively. Interestingly, four and seven Zn-deficiency-responsive miRNAs could be recovered by Zn resupply in shoots and roots, respectively (Additional file 2: Figure S4). For example, osa-miR398a-5p, nov-miR7, nov-miR21, and nov-miR38 were induced by Zn deficiency but were repressed by Zn resupply in shoots; osa-miR818a-3p, osa-miR1428a-3p, and osa-miR5532-5p were induced by Zn deficiency but were repressed by Zn resupply in roots (Fig. 3). Two miRNAs showed similar expression pattern under Zn deficiency and Zn resupply; miR397b-5p was induced by both Zn deficiency and Zn resupply in roots, and nov-miR27 was repressed by Zn deficiency and Zn resupply in roots and shoots, respectively (Fig. 3), suggesting the complex regulations of these miRNAs under Zn deficiency and resupply.

Fig. 3
figure3

Heatmap representation for expression profiles of miRNAs in response to 14 d of Zn deprivation and/or 3 d of Zn resupply in the roots and shoots. The intensities of the color represent the relative magnitude of fold changes in log2 values according to small RNA high-throughput sequencing data. The asterisks indicate significantly differentially expressed miRNAs with an absolute fold change ≥2 and P-value ≤0.05. Red color indicates induction, blue color indicates repression.

Full size image

Identification of potential target genes of Zn-responsive miRNAs

To understand the potential biological function of these Zn-deficiency- and/or Zn-resupply-responsive miRNAs, miRNA-targeted genes were predicted using an online tool psRNATarget [46]. A total of 799 potential target genes were identified for these 68 Zn-responsive miRNAs (Additional file 1: Table S13). At least 19 target genes of the four Zn-responsive miRNAs (miR398, miR408, miR528, miR397) have been previously validated through experiments like RNA ligase-mediated 5′ rapid amplification of cDNA ends (RLM 5′-RACE), degradome sequencing, and genetic analysis (Table 3) [47,48,49]; most of them were Cu-containing proteins and were involved in Cu homeostasis, reactive oxygen species (ROS) homeostasis, photosynthesis and stress tolerance [49,50,51,52,53,54,55,56,57,58].

Table 3 The validated target genes of some Zn-responsive miRNAs
Full size table

To further elucidate the potential biological functions of the predicted target genes in Zn deficiency response, GO enrichment analysis was performed. In total, 26 GO terms of biological process group, 19 GO terms of molecular function group, and 18 GO terms of cellular component group were significantly enriched (FDR < 0.05) ( Fig. 4; Additional file 2: Figure S6 and S7). Among the enriched molecular functions, the most significant GO terms were oxidoreductase activity (GO:0016682), laccase activity (GO:0008471), serine-type peptidase activity (GO:0008236), Cu ion binding (GO:0005507), and serine hydrolase activity (GO:0017171) (Fig. 4). Among the enriched biological processes, the most enriched GO terms were lignin metabolic process (GO:0009808), phenylpropanoid catabolic process (GO:0009808), cellular amino acid derivative catabolic process (GO:0042219), secondary metabolic process (GO:0019748), and proteolysis (GO:0006508) (Additional file 2: Figure S6). In cellular component group, apoplast (GO:0048046), chromatin (GO:0000785), and photosynthetic membrane (GO:0034357) were the most over-represented terms (Additional file 2: Figure S7).

Fig. 4
figure4

Gene ontology (GO) representation of the overrepresented GO terms of molecular functions in the potential target genes of the differentially expressed miRNAs under Zn deprivation and/or Zn resupply. The GO representation was generated using single enrichment analysis (SEA) tool on AgriGO (http://bioinfo.cau.edu.cn/agriGO/) (Fisher’s test, P < 0.05, FDR < 0.05). The number in parenthesis represents the FDR value

Full size image

Integrated analysis of miRNAome and transcriptome

Given that miRNAs negatively regulate the expression of their target mRNAs by target cleavage, the expression patterns of miRNAs generally show an inverse correlation with those of their target genes. By integrating transcriptome and miRNAome data, twelve DEGs under -Zn and/or Zn resupply were predicted to be targeted by 10 miRNAs (including four known and six novel miRNAs) that were responsive to -Zn and/or Zn resupply (Table 4). Half of the potential target DEGs showed a different response to Zn deficiency or Zn resupply as compared with the corresponding miRNAs, suggesting the negative regulation of these miRNAs on their potential target genes. For example, miR397b-5p was up-regulated by Zn resupply, while its potential target genes OsLAC11 (Os03g0273200) and OsLAC29 (Os12g0258700) were up-regulated by -Zn or down-regulated by Zn resupply; osa-miR398a-5p was up-regulated by Zn deficiency, while its potential target gene Os07g0213800 was down-regulated; nov-miR4 was down-regulated by Zn resupply, while its potential target gene Os01g0358700 was up-regulated; nov-miR6 was down-regulated by Zn deficiency, while it potential target gene Os04g0623901 was up-regulated; nov-miR23 was down-regulated by Zn deficiency, while its potential target gene Os01g0923900 was down-regulated by Zn-resupply. However, half of the potential target genes showed similar expression patterns with the miRNAs (Table 4). This may be caused by the following reasons: these genes are inauthentic targets of the miRNAs; there are complicated regulation of the potential target genes; the potential target genes are regulated by translational repression but not by transcript cleavage.

Table 4 DEGs under Zn deficiency and/or Zn resupply that are potential targets of Zn-responsive miRNAs
Full size table

Confirmation of miRNAome and transcriptome results by qRT-PCR

To confirm the small RNA sequencing results, seven Zn-responsive miRNAs were randomly selected for expression validation by qRT-PCR. Primary miRNAs (pri-miRNAs), which contain a characteristic hairpin structure, are transcribed from MIR genes through polymerase II in plants and can be used to indicate the expression level of mature miRNAs [59]. Here, the primary transcripts of these miRNAs were analyzed using gene-specific primers based on the sequences of the corresponding miRNA genes. The expression of pri-miRNAs were almost consistent with the small RNA sequencing results of mature miRNAs. For example, pri-miR398a was up-regulated by Zn deficiency and down-regulated by Zn resupply; pri-miR408 was up-regulated by Zn resupply; pri-miR528 was up-regulated by Zn deficiency (Fig. 5). However, pri-miR397b and pri-nov-miR28 were repressed by Zn deficiency, which is inconsistent with the expression of mature miR397b and nov-miR28, suggesting the different responses of these pri-miRNAs and mature miRNAs under Zn deficiency. To confirm the RNA sequencing results, four Zn-responsive genes (Os05g0472400, Os03g0191900, Os06g0639800, and Os06g0131300) were randomly selected and investigated by qRT-PCR. The expression patterns of these genes under Zn deficiency and recovery were similar to the results of RNA-seq, with the exception of Os06g0131300 (Fig. 5). Significant expression changes of Os06g0131300 in roots under Zn deficiency and recovery could not be observed by qRT-PCR, but its expression was found to be repressed by Zn deficiency in shoots as observed by qRT-PCR. The expression levels of Os03g0191900 in shoots could not be detected by qRT-PCR, which is consistent with the low FPKM value in shoots as revealed by RNA-seq (Fig. 5). These results suggest that the data of small RNA sequencing and mRNA sequencing obtained in this study is reliable.

Fig. 5
figure5

Heatmap of the expression of differentially expressed miRNAs and DEGs revealed by deep sequencing and qRT-PCR. Seven miRNAs and four DEGs were randomly selected and analyzed by qRT-PCR. The primary transcripts of miRNAs (pri-miRNAs) were analyzed to for miRNA expression validation. The intensities of the color represent the fold changes in log2 values according to qRT-PCR results of three replicates. The asterisks indicate an absolute fold change ≥2 and P-value ≤0.05. Red color indicates induction, blue color indicates repression, and gray color indicates no expression detected

Full size image

We also analyzed the expression of several potential target genes of Zn-responsive miRNAs by qRT-PCR. The potential target genes of miR397b-5p showed differential responses to Zn deficiency; OsLAC7 (Os01g0850550) and OsLAC9 (Os01g0850800) were repressed by Zn deficiency in roots, which is contrary to the expression of miR397b-5p; OsLAC7 was induced by Zn deficiency and repressed by Zn resupply in shoots (Fig. 6a). OsCDS2, the potential target gene of miR398a, was found to be repressed by Zn deficiency in roots, which is contrary to the expression of miR398a (Fig. 6b). The potential targets of miR408-3p also showed differential responses to Zn deficiency and Zn resupply, suggesting their different functions in Zn deficiency response. OsUCL16 (Os06g0218600) was up-regulated by Zn deficiency and down-regulated by Zn resupply; OsUCL30 (Os08g0482700) was down-regulated by Zn deficiency in roots; OsORC6 (Os07g0628600) was down-regulated by Zn deficiency in both roots and shoots (Fig. 6c).

Fig. 6
figure6

The relative expression levels of three Zn-responsive miRNAs and their potential target genes. The relative expression of potential target genes of miR397b-5p (a), miR398a (b), and miR408-3p (c) were analyzed by qRT-PCR. The normalized expression levels of miRNAs revealed by small RNA sequencing were shown with line charts. The expression level of +Zn control in roots and shoots was normalized to one. Error bar represents standard deviation of three biological replicates. The asterisk above error bar indicates an absolute fold change ≥2 (−Zn/+Zn, or Zn resupply/−Zn) and P-value ≤0.05. ND means “not detected”

Full size image

Discussion

Identification of DEGs in response to Zn deficiency and recovery

In this study, a total of 569 loci were found to be responsive to Zn deficiency and/or Zn resupply; 360 of them were responsive to Zn deficiency in roots and/or shoots, and 316 of them were responsive to Zn resupply in roots and/or shoots (Fig. 2). By comparing with earlier studies [24], 49 of the DEGs were reported to be responsive to Zn deficiency and/or Zn resupply (Additional file 1: Table S14). Five of the common DEGs are putative Zn transporters, and three of the common DEGs are metallothionein-like proteins, which are implicated in Zn homeostasis [7, 60]. In addition, 97 DEGs were found to be recovered by Zn resupply treatment in shoots and/or roots (Tables 1 and 2), suggesting their specific responses to Zn deficiency. These DEGs were enriched in GO terms of cellular metabolic process, primary metabolic process and biosynthetic process (Additional file 2: Figure S8), suggesting that these processes are affected by Zn deficiency and can be recovered by Zn resupply. Four of the DEGs (Os03g0108300, Os10g0528400, Os12g0570700, Os12g0571000) were also reported to be recovered by Zn resupply [24].

Twelve DEGs were predicted to be targeted by 10 Zn-responsive miRNAs, and half of them showed a differential expression with miRNAs under Zn deficiency and/or Zn resupply (Table 4), suggesting the interactions between miRNAs and their potential target genes. However, this cannot exclude the responsiveness of other potential target genes of Zn-responsive miRNAs. As revealed by qRT-PCR, several potential target genes of miR397b, miR408 and miR398a were also found to be induced or repressed by Zn deficiency and/or Zn resupply (Table 3, Fig. 5b). The number of DEGs potentially regulated by Zn-responsive miRNAs in this study is far below our expectation. One of the possible reasons is that the accepted standard used for the definition of DEGs may miss some interactions between miRNAs and their potential target genes, and more interaction pairs may be identified by lowering the threshold.

miRNAs potentially involved in the regulation of Zn deficiency response in rice

miRNAs have been demonstrated to be involved in signaling and regulation of nutrient stress responses, like nitrogen deficiency, phosphate starvation, sulfate deprivation, and Cu deficiency [38,39,40, 61, 62]. In this study, we analyzed the miRNA expression profilings under Zn deficiency and Zn resupply, and 68 miRNAs were found to be responsive to Zn deficiency and/or Zn resupply in rice (Fig. 3). Some of the miRNA families were also found to be responsive to Zn deficiency in Sorghum bicolor, like miR171, miR398, miR399, miR408, and miR528 [2], suggesting their conserved roles in plant response to Zn deficiency.

Among the Zn-responsive miRNAs, 38 (56%) of them were novel miRNAs. Most of these novel miRNAs (36/38) are encoded by a single locus (Additional file 1: Table S12), suggesting their recent evolutionary origin. Most of the novel miRNAs (30/38) contained an opposite strand (miRNA*) that can be detected by small RNA sequencing, and most of the novel miRNAs (28/38) were 24 nt in length (Additional file 1: Table S12; Additional file 2: Figure S5). In rice, a proportion of pri-miRNAs can be processed into canonical miRNAs (21 nt) and long miRNAs (24 nt) by Dicer-like 1 (DCL1) and DCL3, respectively [63, 64]. The 24 nt long miRNAs are sorted into the effector Argonaute 4 (AGO4) and possibly direct DNA methylation by base pairing with target loci-derived transcripts [26, 63]. DNA methylation modulates gene expression and represses the transcription and movement of transposable elements. The expression of some key nutrient stress-responsive genes has been suggested to be regulated by DNA methylation under specific nutrient stress conditions [65]. In this study, many of the Zn-responsive novel miRNAs (20/38) were down-regulated by Zn deficiency and/or up-regulated by Zn resupply (Fig. 3). The potential target genes of these Zn-responsive novel miRNAs encode various kinds of proteins, such as transcription factors, transporters, protein kinases, Zn-containing protein, and oxidoreductases (Additional file 1: Table S13). Whether these recently-evolved and possibly species-specific Zn-responsive novel miRNAs are associated with acclimation to submergence and low Zn conditions in rice, and are involved in epigenetic regulation of Zn deficiency response deserve further investigations.

Potential role of miRNAs in Cu-Zn interaction in plants

Among the Zn-responsive miRNAs identified in this study, some of the miRNA families were previously documented to be responsive to Cu deficiency, such as miR397, miR398, and miR408, which are all induced by Cu deficiency and their target genes encode various Cu-containing proteins [39]. In rice, nineteen potential target genes of the Cu/Zn-responsive miRNAs have been validated, most of them encode Cu-containing proteins, like CSDs, laccases, plantacyanin like proteins, and uclacyanin-like proteins (Table 3). The potential target genes of Zn-responsive miRNAs were significantly enriched in the GO terms of Cu ion binding and laccase activity (Fig. 4). These results suggest the involvement of miRNAs in the interaction between Cu and Zn nutrition in plants.

It has been demonstrated that Cu and Zn could interact in several ways: Cu competitively inhibits Zn absorption, Zn strongly depresses Cu absorption, and Cu affects the redistribution of Zn within plants [66, 67]. The uptakes of Cu and iron (Fe) were found to be increased by Zn deficiency in rice seedlings [68]. In this study, miR398 and miR528 were found to be induced by Zn deficiency. The induction of miR398 and miR528 could depress the expressions of Zn-containing target gene CSDs [2, 50], which can be functionally replaced by Fe superoxide dismutase and Mn superoxide dismutase, thus further allow the allocation of limited Zn to other essential Zn-containing proteins in order to adapt to the low Zn environment. The miRNA-mediated regulation of Zn homeostasis could be similar to that of the miRNA-mediated regulation of Cu homeostasis as reported previously [69, 70]. Interestingly, another two Cu-responsive miRNAs, miR397 and miR408 were found to be induced by Zn resupply (Fig. 3). The induction of these two miRNAs may be a direct effect of Zn resupply, and may also be an indirect effect of Zn-resupply-induced antagonism on Cu nutrition. The primary transcripts of miR397b and miR408 were found to be repressed by Zn deficiency (Fig. 5). Whether the repressions of miR397b and miR408 would be associated with the Zn deficiency-induced Cu accumulation in plant tissues deserve further studies.

miRNA-mediated oxidative stress in response to Zn deficiency

Zn deficiency can interfere with membrane-bound NADPH oxidase producing ROS and further cause oxidative damage to critical cell compounds by the excessive ROS, which is a major factor affecting plant growth under Zn deficiency [11]. It has been reported that Zn deficiency causes cellular damage in leaves and roots and maintaining a higher ROS defense level is associated with reduced cellular damage in rice leaves [2, 71, 72]. In this study, GO terms of oxidoreductase activity was significantly enriched in the potential target genes of Zn-responsive miRNAs (Fig. 4). miR398 was reported to be involved in the regulation of oxidative stress response by silencing the expression of CSD1 and CSD2 in Arabidopsis [73]. Here, miR398 was induced by Zn deficiency in the shoots. The induction of miR398 could repress the expression of its target genes (CSDs, SODX, CCSD) and then promote the accumulation of ROS. miR398 was also reported to be induced by other stresses, such as salt, heat and infection of blast fungus Magnaporthe oryzae [29, 52]. The induction of miR398 has been suggested to be required for thermotolerance in Arabidopsis and resistance to blast disease in rice [51, 52]. Another Zn-responsive miRNA, miR528 was also reported to be involved in the regulation of ROS accumulation [55, 56]. One of the target gene of miR528, CSD4 (Os08g0561700) is possibly involved in ROS scavenging. Another target gene of miR528, Os06g0567900, encoding a L-ascorbate oxidase (AO), is involved in producing ROS by regulating the redox state of the apoplast by oxidizing apoplastic L-ascorbate acid [56]. Overexpression of osa-miR528 can decrease AO activity and increase ROS scavenging and further improve plant tolerance to salt and nitrogen deficiency stresses in creeping bentgrass [55]. In this study, miR528 was found to be induced by Zn deficiency, suggesting its role in balancing ROS accumulation and scavenging under Zn stress.

Conclusions

The present study is the first attempt to integrate global miRNA and mRNA expression profiles to identify regulatory miRNA/target modules in response to Zn deficiency and recovery in rice. A set of 68 miRNAs including 38 novel miRNAs and five Cu-responsive miRNAs were identified to be differentially expressed. At least 12 potential target genes of Zn-responsive miRNAs were found to be altered by Zn deficiency and recovery. Ninety-seven Zn-deficiency-responsive genes were found to be recovered by Zn resupply. This global analysis of the rice miRNAome and transcriptome in responses to Zn deficiency and recovery identifies candidate miRNAs and genes that can be further characterized to increase our understanding of the molecular mechanisms of Zn deficiency response. The identified differentially expressed miRNA/target modules and genes can also be considered as candidates for improving Zn deficiency tolerance by genetic engineering approaches.

Methods

Plant material and growth conditions

Seeds of the rice (Oryza sativa L. ssp. japonica cv. Nipponbare) from our own lab were surface sterilized in 10% (v/v) H2O2 for 20 min and rinsed with sterile water [74]. Seeds were then transferred to seedling trays floating on sterile water for germination. After growing in the sterile water for 10 days, seedlings were transferred to a plastic container with 2.5 L Yoshida solution [75]. After growing in the nutrient solution (with 1.5 μM ZnSO4) for 1 week, uniform seedlings were used for Zn deficiency (−Zn) treatment. For -Zn treatment, no ZnSO4 was added to the nutrient solution. The fresh nutrient solution was replenished every 2 day, and the pH was adjusted to 5.6. Plants were grown in a growth chamber under controlled conditions. The light intensity was approximately 180 μmol m− 2 s− 1 at shoot height with a day/night cycle of 16 h/8 h at 26 °C/24 °C. The relative humidity was 50%. Roots and shoots of rice seedlings were collected after Zn-deficient treatment for 2 weeks. For Zn-resupply treatment, Zn-deficient seedlings were transferred to normal nutrient solution containing 1.5 μM ZnSO4 for 3 days. After Zn-resupply treatment, roots and shoots were collected for analyses at the same time point as collecting the Zn-deficient samples (4 hours after illumination).

Zn content determination

Zn content was analyzed as described previously [23, 76]. Plant tissues were collected and dried in an oven at 80 °C for 3 days. The roots were rinsed with deionized water for three times before collection. The dried roots or shoots were ground, and about 0.10 g sample was used for analysis. After digestion with HNO3/HClO4 (4:1, v/v) at 180 °C for 4 h, the solution was diluted five times with 1% HNO3, and then the concentration of Zn was determined by Inductively Coupled Plasma Mass Spectroscopy (Perkin-Elmer NexION 300X, Waltham, MA).

cDNA library construction and RNA sequencing

Rice shoots and roots were collected after Zn-deficiency and Zn-resupply treatments. Three rice seedlings were collected and mixed together as one sample. Eighteen samples (Zn plus shoot, ZPS; Zn minus shoot, ZMS; Zn plus root, ZPR; Zn minus root, ZMR; Zn resupply shoot, ZRS; Zn resupply root, ZRR, each with three biological replicates) were harvested for RNA library construction and sequencing. Total RNA was extracted using Trizol reagent (Invitrogen, CA, USA), and the quantity and quality were analyzed by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA). The RIN value (RNA integrity number) of all RNA samples were more than 8.0. Ten μg of total RNA of each sample was subjected to isolate Poly (A) mRNA with poly-T oligo-attached magnetic beads (Invitrogen, Carlsbad, CA). After purification, the mRNA is fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments were reverse-transcribed to create cDNA libraries according to the instructions of the mRNA-Seq sample preparation kit (Illumina, San Diego, USA). The average insert size of the paired-end libraries was 300 ± 50 bp. The paired-end sequencing was performed on an Illumina Hiseq4000 at the LC-BIO (Hangzhou, China). The RNA-seq datasets of this study were deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO accession number GSE130980 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE130980).

Mapping of RNA-seq reads and identification of differentially expressed genes

The initial base calling and quality filtering of the reads generated with the Illumina analysis pipeline (Fastq format) were conducted using a custom Perl script and the default parameters of the Illumina pipeline (http://www.illumina.com). Additional filtering for poor-quality bases was performed using the FASTX-toolkit available in the FastQC software package (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The rice reference genome (version IRGSP-1.0) (https://rapdb.dna.affrc.go.jp/) was indexed by Bowtie2 to facilitate the read mapping. The read mapping was performed using the Tophat software (version 2.1.1) package [77]. Tophat allows multiple alignments per read (up to 40) and a maximum of two mismatches when mapping the reads to the reference genome. Reads were first mapped directly to the genome using indexing and the unmapped reads were used to identify novel splicing events. The aligned read files were processed by Cufflinks to measure the relative abundances of the transcripts by using the normalized RNA-seq fragment counts. The estimated gene abundance was measured in terms of the fragments per kilobase of transcript per million mapped reads (FPKM). The DEGs between the two treatments were identified using Cuffdiff (version 2.1.1). Only the genes with a log2 fold change ≥1 or ≤ − 1, and an adjusted P-value (FDR ≤ 0.05) were considered as significantly DEGs. The Pearson correlation was calculated and plotted using the R Stats package (version 3.5.0) based on the normalized counts of mRNAs of each sample.

Functional annotation and gene ontology (GO) enrichment

Zn-deficiency responsive genes and potential target genes of miRNAs were annotated for gene ontology (GO) terms [78] and categorized into molecular function, cellular component, and biological process categories. GO term enrichment was conducted using single enrichment analysis tool on AgriGo (http://bioinfo.cau.edu.cn/agriGO/) [79]. GO category (http://geneontology.org/) with an adjusted P-value ≤0.05 was regarded as significantly enriched.

Construction of small RNA libraries and sequencing

Eighteen samples (ZPS, ZMS, ZPR, ZMR, ZRS, and ZRR, each with three biological replicates) were harvested for small RNA library construction and sequencing. Total RNAs were extracted from samples using Trizol reagent (Invitrogen, Carlsbad, CA). Approximately 1.0 μg of total RNA was used to prepare a small RNA library according to the instruction of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA). Single-end sequencing with a length of about 50 bp was performed on an Illumina Hiseq2500 at the LC-BIO (Hangzhou, China). The small RNA sequencing datasets of this study were deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO accession number GSE131003 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE131003).

Identification of miRNAs and prediction of their target genes

The raw reads were performed on an in-house program ACGT101-miR (LC Sciences, Houston, Texas, USA) to remove adaptor sequences, low-complexity sequences, junk reads, repeat sequences, and the reads that matched the common non-coding RNA families (tRNA, rRNA, snoRNA, and snRNA) deposited in Rfam database (http://www.sanger.ac.uk/software/Rfam). Subsequently, unique sequences with length of 20–24 nucleotide (nt) were mapped to rice miRNA precursors in miRBase 22.0 by BLAST (version 2.6.0) to identify known miRNAs and novel 3p- and 5p- derived miRNAs. Length variation at both 3′ and 5′ ends, and one nt mismatch of the miRNA sequence for the alignment were allowed. The unique sequences mapping to rice mature miRNAs in the precursors were identified as known miRNAs. The unique sequences mapping to the opposite arm of known rice miRNA precursor hairpin were considered to be novel 5p- or 3p- derived miRNA candidates. The remaining sequences were mapped to miRNA precursors of other plant species in miRBase 22.0 (http://www.mirbase.org/) by BLAST (version 2.6.0), and the mapped pre-miRNAs were further aligned with the rice genome to determine their genomic locations. The remaining unmapped sequences were further aligned with the rice genome, and the flanking sequences with a length of 120 nt were collected for secondary structure prediction using RNAfold software (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi). The criteria for novel miRNA prediction and characterization were based on previous study [80].

The abundance of each miRNA were normalized to the expression of tags per million (TPM) following the formula: normalized expression = actual miRNA count/total count of clean reads*106. To identify Zn deficiency- or Zn resupply- responsive miRNAs, only the miRNAs with a log2 fold change ≥1 or ≤ − 1 and a p-value ≤0.05 were considered as significantly differentially expressed miRNAs. To predict potential target genes of Zn-responsive miRNAs, psRNATarget (a plant small RNA target analysis server) (http://plantgrn.noble.org/psRNATarget/) (V2, 2017 release) was employed based on the imperfect complementation between the sequences of miRNAs and their target genes [46].

Quantitative RT-PCR analysis

Quantitative RT-PCR (qRT-PCR) was performed on a MyiQ Single Color Real-time PCR system (Bio-Rad) as described previously [81]. qRT-PCR analysis of primary miRNA (pri-miRNA) were performed as described previously [59, 82]. The stem-loop sequences were used for primer design, and if no satisfactory primers could be found, stem-loop sequences were elongated by 100 bp of flanking genomic sequences on each side for primer design. Relative expression levels were normalized to that of two internal control genes OsACTIN1 (Os03g0718100) and OsUBQ2 (Os02g0161900) using the Pfaffl method (Ratio = (Etarget)∆CTtarget(control-sample)/(Eref) ∆CTref(control-sample)) [83]. The calculated efficiency (E) of all primers was between 1.8 and 2.2. The sequences of primers are listed in Additional file 1: Table S1.

Availability of data and materials

The datasets used during the current study are available from the corresponding author on reasonable request.

Abbreviations

AGO4:

Argonaute 4

AO:

L-ascorbate oxidase

bZIP:

Basic leucine zipper

CCSD:

Copper chaperone for superoxide dismutase

CSD:

Copper/zinc superoxide dismutase

Cu:

Copper

DCL1:

Dicer-like 1

DEG:

Differentially expressed gene

Fe:

Iron

FPKM:

Fragments per kilobase of transcript per million mapped reads

GA:

Gibberellin

GO:

Gene ontology

H3K9:

Histone H3 lysine 9

miRNA:

microRNA

NAS:

Nicotianamine synthase

qRT-PCR:

Quantitative RT-PCR

RLM 5′-RACE:

RNA ligase-mediated 5′ rapid amplification of cDNA ends

ROS:

Reactive oxygen species

SODX:

Superoxide dismutaseX

TF:

Transcription factor

ZDRE:

Zinc deficiency response element

ZIFL:

Zinc-induced facilitator

ZIP:

ZRT, IRT-related protein

Zn:

Zinc

References

  1. 1.

    Broadley MR, White PJ, Hammond JP, Zelko I, Lux A. Zinc in plants. New Phytol. 2007;173(4):677–702.

    CAS  PubMed  Article  Google Scholar 

  2. 2.

    Li Y, Zhang Y, Shi D, Liu X, Qin J, Ge Q, Xu L, Pan X, Li W, Zhu Y, et al. Spatial-temporal analysis of zinc homeostasis reveals the response mechanisms to acute zinc deficiency in Sorghum bicolor. New Phytol. 2013;200(4):1102–15.

    CAS  PubMed  Article  Google Scholar 

  3. 3.

    Alloway B. Soil factors associated with zinc deficiency in crops and humans. Environ Geochem Health. 2009;31(5):537–48.

    CAS  PubMed  Article  Google Scholar 

  4. 4.

    Cakmak I. Enrichment of cereal grains with zinc: agronomic or genetic biofortification? Plant Soil. 2008;302(1–2):1–17.

    CAS  Google Scholar 

  5. 5.

    Brown KH, Rivera JA, Bhutta Z, Gibson RS, King JC, Lonnerdal B, Ruel MT, Sandtrom B, Wasantwisut E, Hotz C. International zinc nutrition consultative group (IZiNCG) technical document #1. Assessment of the risk of zinc deficiency in populations and options for its control. Food Nutr Bull. 2004;25(1 Suppl 2):S99–203.

    PubMed  Google Scholar 

  6. 6.

    Rose TJ, Impa SM, Rose MT, Pariasca-Tanaka J, Mori A, Heuer S, Johnson-Beebout SE, Wissuwa M. Enhancing phosphorus and zinc acquisition efficiency in rice: a critical review of root traits and their potential utility in rice breeding. Ann Bot. 2012;112:331–45.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  7. 7.

    Sinclair SA, Krämer U. The zinc homeostasis network of land plants. Biochimica et Biophysica Acta (BBA) – Mol Cell Res. 2012;1823(9):1553–67.

    CAS  Article  Google Scholar 

  8. 8.

    Ishimaru Y, Bashir K, Nishizawa NK. Zn uptake and translocation in rice plants. Rice. 2011;4(1):21–7.

    Article  Google Scholar 

  9. 9.

    White JG, Zasoski RJ. Mapping soil micronutrients. Field Crop Res. 1999;60(1):11–26.

    Article  Google Scholar 

  10. 10.

    Lee J-S, Sajise AGC, Gregorio GB, Kretzschmar T, Ismail AM, Wissuwa M. Genetic dissection for zinc deficiency tolerance in rice using bi-parental mapping and association analysis. Theor Appl Genet. 2017;130(9):1903–14.

    CAS  PubMed  Article  Google Scholar 

  11. 11.

    Cakmak I. Possible roles of zinc in protecting plant cells from damage by reactive oxygen species. New Phytol. 2000;146(2):185–205.

    CAS  Article  Google Scholar 

  12. 12.

    Impa SM, Johnson-Beebout SE. Mitigating zinc deficiency and achieving high grain Zn in rice through integration of soil chemistry and plant physiology research. Plant Soil. 2012;361(1):3–41.

    CAS  Article  Google Scholar 

  13. 13.

    Hacisalihoglu G, Hart JJ, Wang Y-H, Cakmak I, Kochian LV. Zinc efficiency is correlated with enhanced expression and activity of zinc-requiring enzymes in wheat. Plant Physiol. 2003;131(2):595.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Cavagnaro T. The role of arbuscular mycorrhizas in improving plant zinc nutrition under low soil zinc concentrations: a review. Plant Soil. 2008;304(1–2):315–25.

    CAS  Article  Google Scholar 

  15. 15.

    Widodo B, Broadley MR, Rose T, Frei M, Pariasca-Tanaka J, Yoshihashi T, Thomson M, Hammond JP, Aprile A, Close TJ, et al. Response to zinc deficiency of two rice lines with contrasting tolerance is determined by root growth maintenance and organic acid exudation rates, and not by zinc-transporter activity. New Phytol. 2010;186(2):400–14.

    PubMed  Article  CAS  Google Scholar 

  16. 16.

    Mori A, Kirk GJ, Lee JS, Morete MJ, Nanda AK, Johnson-Beebout SE, Wissuwa M. Rice genotype differences in tolerance of zinc-deficient soils: evidence for the importance of root-induced changes in the Rhizosphere. Front Plant Sci. 2015;6:1160.

    PubMed  Google Scholar 

  17. 17.

    Ishimaru Y, Masuda H, Suzuki M, Bashir K, Takahashi M, Nakanishi H, Mori S, Nishizawa NK. Overexpression of the OsZIP4 zinc transporter confers disarrangement of zinc distribution in rice plants. J Exp Bot. 2007;58(11):2909–15.

    CAS  PubMed  Article  Google Scholar 

  18. 18.

    Olsen LI, Hansen TH, Larue C, Osterberg JT, Hoffmann RD, Liesche J, Kramer U, Surble S, Cadarsi S, Samson VA, et al. Mother-plant-mediated pumping of zinc into the developing seed. Nat Plants. 2016;2(5):16036.

    CAS  PubMed  Article  Google Scholar 

  19. 19.

    Cai H, Huang S, Che J, Yamaji N, Ma JF. The tonoplast-localized transporter OsHMA3 plays an important role in maintaining Zn homeostasis in rice. J Exp Bot. 2019;70(10):2717–25.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. 20.

    Sasaki A, Yamaji N, Mitani-Ueno N, Kashino M, Ma JF. A node-localized transporter OsZIP3 is responsible for the preferential distribution of Zn to developing tissues in rice. Plant J. 2015;84(2):374–84.

    CAS  PubMed  Article  Google Scholar 

  21. 21.

    Sinclair SA, Senger T, Talke IN, Cobbett CS, Haydon MJ, Kramer U. Systemic Upregulation of MTP2- and HMA2-mediated Zn partitioning to the shoot supplements local Zn deficiency responses. Plant Cell. 2018;30(10):2463–79.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. 22.

    Nanda AK, Pujol V, Wissuwa M. Patterns of stress response and tolerance based on transcriptome profiling of rice crown tissue under zinc deficiency. J Exp Bot. 2017;68(7):1715–29.

    CAS  PubMed  Article  Google Scholar 

  23. 23.

    Zeng H, Zhang X, Ding M, Zhang X, Zhu Y. Transcriptome profiles of soybean leaves and roots in response to zinc deficiency. Physiol Plant. 2019;167:330–51.

    CAS  PubMed  Article  Google Scholar 

  24. 24.

    Bandyopadhyay T, Mehra P, Hairat S, Giri J. Morpho-physiological and transcriptome profiling reveal novel zinc deficiency-responsive genes in rice. Funct Integr Genomics. 2017;17(5):565–81.

    CAS  PubMed  Article  Google Scholar 

  25. 25.

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

    CAS  PubMed  Article  Google Scholar 

  26. 26.

    Song X, Li Y, Cao X, Qi Y. MicroRNAs and their regulatory roles in plant-environment interactions. Annu Rev Plant Biol. 2019;70:489–525.

    CAS  PubMed  Article  Google Scholar 

  27. 27.

    Tang J, Chu C. MicroRNAs in crop improvement: fine-tuners for complex traits. Nat Plants. 2017;3(7):17077.

    PubMed  Article  Google Scholar 

  28. 28.

    Li S, Castillo-Gonzalez C, Yu B, Zhang X. The functions of plant small RNAs in development and in stress responses. Plant J. 2017;90(4):654–70.

    CAS  PubMed  Article  Google Scholar 

  29. 29.

    Khraiwesh B, Zhu JK, Zhu J. Role of miRNAs and siRNAs in biotic and abiotic stress responses of plants. Biochim Biophys Acta. 2012;1819(2):137–48.

    CAS  PubMed  Article  Google Scholar 

  30. 30.

    Gao S, Yang L, Zeng HQ, Zhou ZS, Yang ZM, Li H, Sun D, Xie F, Zhang B. A cotton miRNA is involved in regulation of plant response to salt stress. Sci Rep. 2016;6:19736.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. 31.

    Meng Y, Shao C, Wang H, Chen M. The regulatory activities of plant microRNAs: a more dynamic perspective. Plant Physiol. 2011;157(4):1583–95.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. 32.

    Chen Z, Hu L, Han N, Hu J, Yang Y, Xiang T, Zhang X, Wang L. Overexpression of a miR393-resistant form of transport inhibitor response protein 1 (mTIR1) enhances salt tolerance by increased osmoregulation and Na+ exclusion in Arabidopsis thaliana. Plant Cell Physiol. 2015;56(1):73–83.

    CAS  PubMed  Article  Google Scholar 

  33. 33.

    Shao Y, Zhou HZ, Wu Y, Zhang H, Lin J, Jiang X, He Q, Zhu J, Li Y, Yu H, et al. OsSPL3, an SBP-domain protein, regulates crown root development in Rice. Plant Cell. 2019;31(6):1257–75.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. 34.

    Hong Y, Jackson S. Floral induction and flower formation--the role and potential applications of miRNAs. Plant Biotechnol J. 2015;13(3):282–92.

    CAS  PubMed  Article  Google Scholar 

  35. 35.

    Fang Y, Xie K, Xiong L. Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot. 2014;65(8):2119–35.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. 36.

    Kim JH, Woo HR, Kim J, Lim PO, Lee IC, Choi SH, Hwang D, Nam HG. Trifurcate feed-forward regulation of age-dependent cell death involving miR164 in Arabidopsis. Science. 2009;323(5917):1053–7.

    CAS  PubMed  Article  Google Scholar 

  37. 37.

    Jiang D, Chen W, Dong J, Li J, Yang F, Wu Z, Zhou H, Wang W, Zhuang C. Overexpression of miR164b-resistant OsNAC2 improves plant architecture and grain yield in rice. J Exp Bot. 2018;69(7):1533–43.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  38. 38.

    Zeng H, Wang G, Hu X, Wang H, Du L, Zhu Y. Role of microRNAs in plant responses to nutrient stress. Plant Soil. 2014;374:1005–21.

    CAS  Article  Google Scholar 

  39. 39.

    Pilon M. The copper microRNAs. New Phytol. 2017;213(3):1030–5.

    CAS  PubMed  Article  PubMed Central  Google Scholar 

  40. 40.

    Liu TY, Lin WY, Huang TK, Chiou TJ. MicroRNA-mediated surveillance of phosphate transporters on the move. Trends Plant Sci. 2014;19(10):647–55.

    CAS  PubMed  Article  Google Scholar 

  41. 41.

    Wang L, Zeng HQ, Song J, Feng SJ. Yang ZM: miRNA778 and SUVH6 are involved in phosphate homeostasis in Arabidopsis. Plant Sci. 2015;238:273–85.

    CAS  PubMed  Article  Google Scholar 

  42. 42.

    Hu B, Zhu C, Li F, Tang J, Wang Y, Lin A, Liu L, Che R, Chu C. LEAF TIP NECROSIS1 plays a pivotal role in the regulation of multiple phosphate starvation responses in rice. Plant Physiol. 2011;156(3):1101–15.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  43. 43.

    Lin WY, Huang TK, Chiou TJ. Nitrogen limitation adaptation, a target of microRNA827, mediates degradation of plasma membrane-localized phosphate transporters to maintain phosphate homeostasis in Arabidopsis. Plant Cell. 2013;25(10):4061–74.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. 44.

    Zhang H, Zhao X, Li J, Cai H, Deng XW, Li L. MicroRNA408 is critical for the HY5-SPL7 gene network that mediates the coordinated response to light and copper. Plant Cell. 2014;26(12):4933–53.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  45. 45.

    Jeong DH, Park S, Zhai J, Gurazada SG, De Paoli E, Meyers BC, Green PJ. Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. Plant Cell. 2011;23(12):4185–207.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. 46.

    Dai X, Zhao PX. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 2011;39(Web Server issue):W155–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  47. 47.

    Li YF, Zheng Y, Addo-Quaye C, Zhang L, Saini A, Jagadeeswaran G, Axtell MJ, Zhang W, Sunkar R. Transcriptome-wide identification of microRNA targets in rice. Plant J. 2010;62(5):742–59.

    CAS  PubMed  Article  Google Scholar 

  48. 48.

    Zhou M, Gu L, Li P, Song X, Wei L, Chen Z, Cao X. Degradome sequencing reveals endogenous small RNA targets in rice (Oryza sativa L. ssp. indica). Front Biol. 2010;5(1):67–90.

    CAS  Article  Google Scholar 

  49. 49.

    Pan J, Huang D, Guo Z, Kuang Z, Zhang H, Xie X, Ma Z, Gao S, Lerdau MT, Chu C, et al. Overexpression of microRNA408 enhances photosynthesis, growth, and seed yield in diverse plants. J Integr Plant Biol. 2018;60(4):323–40.

    CAS  PubMed  Article  Google Scholar 

  50. 50.

    Lu Y, Feng Z, Bian L, Xie H. Liang J: miR398 regulation in rice of the responses to abiotic and biotic stresses depends on CSD1 and CSD2 expression. Funct Plant Biol. 2011;38(1):44.

    CAS  Article  Google Scholar 

  51. 51.

    Guan Q, Lu X, Zeng H, Zhang Y, Zhu J. Heat stress induction of miR398 triggers a regulatory loop that is critical for thermotolerance in Arabidopsis. Plant J. 2013;74(5):840–51.

    CAS  PubMed  Article  Google Scholar 

  52. 52.

    Li Y, Cao XL, Zhu Y, Yang XM, Zhang KN, Xiao ZY, Wang H, Zhao JH, Zhang LL, Li GB, et al. Osa-miR398b boosts H2O2 production and rice blast disease-resistance via multiple superoxide dismutases. New Phytol. 2019;222(3):1507–22.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  53. 53.

    Ma C, Burd S. Lers a: miR408 is involved in abiotic stress responses in Arabidopsis. Plant J. 2015;84(1):169–87.

    CAS  PubMed  Article  Google Scholar 

  54. 54.

    Zhang JP, Yu Y, Feng YZ, Zhou YF, Zhang F, Yang YW, Lei MQ, Zhang YC, Chen YQ. MiR408 regulates grain yield and photosynthesis via a Phytocyanin protein. Plant Physiol. 2017;175(3):1175–85.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  55. 55.

    Yuan S, Li Z, Li D, Yuan N, Hu Q, Luo H. Constitutive expression of Rice MicroRNA528 alters plant development and enhances tolerance to salinity stress and nitrogen starvation in creeping Bentgrass. Plant Physiol. 2015;169(1):576–93.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  56. 56.

    Wu J, Yang R, Yang Z, Yao S, Zhao S, Wang Y, Li P, Song X, Jin L, Zhou T, et al. ROS accumulation and antiviral defence control by microRNA528 in rice. Nat Plants. 2017;3:16203.

    CAS  PubMed  Article  Google Scholar 

  57. 57.

    Zhang YC, Yu Y, Wang CY, Li ZY, Liu Q, Xu J, Liao JY, Wang XJ, Qu LH, Chen F, et al. Overexpression of microRNA OsmiR397 improves rice yield by increasing grain size and promoting panicle branching. Nat Biotechnol. 2013;31(9):848–52.

    CAS  PubMed  Article  Google Scholar 

  58. 58.

    Xue C, Yao JL, Qin MF, Zhang MY, Allan AC, Wang DF, Wu J. PbrmiR397a regulates lignification during stone cell development in pear fruit. Plant Biotechnol J. 2019;17(1):103–17.

    CAS  PubMed  Article  Google Scholar 

  59. 59.

    Pant BD, Musialak-Lange M, Nuc P, May P, Buhtz A, Kehr J, Walther D, Scheible WR. Identification of nutrient-responsive Arabidopsis and rapeseed microRNAs by comprehensive real-time polymerase chain reaction profiling and small RNA sequencing. Plant Physiol. 2009;150(3):1541–55.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  60. 60.

    Lee S, Jeong HJ, Kim SA, Lee J, Guerinot ML, An G. OsZIP5 is a plasma membrane zinc transporter in rice. Plant Mol Biol. 2010;73(4):507–17.

    CAS  PubMed  Article  Google Scholar 

  61. 61.

    Liang G, Yang F, Yu D. MicroRNA395 mediates regulation of sulfate accumulation and allocation in Arabidopsis thaliana. Plant J. 2010;62(6):1046–57.

    CAS  PubMed  Google Scholar 

  62. 62.

    Lu YB, Qi YP, Yang LT, Guo P, Li Y, Chen LS. Boron-deficiency-responsive microRNAs and their targets in Citrus sinensis leaves. BMC Plant Biol. 2015;15:271.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  63. 63.

    Wu L, Zhou H, Zhang Q, Zhang J, Ni F, Liu C, Qi Y. DNA methylation mediated by a microRNA pathway. Mol Cell. 2010;38(3):465–75.

    CAS  PubMed  Article  Google Scholar 

  64. 64.

    Cuperus JT, Fahlgren N, Carrington JC. Evolution and functional diversification of MIRNA genes. Plant Cell. 2011;23(2):431–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  65. 65.

    Secco D, Whelan J, Rouached H, Lister R. Nutrient stress-induced chromatin changes in plants. Curr Opin Plant Biol. 2017;39:1–7.

    CAS  PubMed  Article  Google Scholar 

  66. 66.

    Hafeez B, Khanif Y, Saleem M. Role of zinc in plant nutrition-a review. Am J Exp Agri. 2013;3(2):374–91.

    CAS  Google Scholar 

  67. 67.

    Kausar MA, Chaudhry FM, Rashid A, Latif A, Alam SM. Micronutrient availability to cereals from calcareous soils. Plant Soil. 1976;45(2):397–410.

    CAS  Article  Google Scholar 

  68. 68.

    Dong C, He F, Berkowitz O, Liu J, Cao P, Tang M, Shi H, Wang W, Li Q, Shen Z, et al. Alternative splicing plays a critical role in maintaining mineral nutrient homeostasis in Rice (Oryza sativa). Plant Cell. 2018;30(10):2267–85.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  69. 69.

    Yamasaki H, Abdel-Ghany SE, Cohu CM, Kobayashi Y, Shikanai T, Pilon M. Regulation of copper homeostasis by micro-RNA in Arabidopsis. J Biol Chem. 2007;282(22):16369–78.

    CAS  PubMed  Article  Google Scholar 

  70. 70.

    Abdel-Ghany SE, Pilon M. MicroRNA-mediated systemic down-regulation of copper protein expression in response to low copper availability in Arabidopsis. J Biol Chem. 2008;283(23):15932–45.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. 71.

    Frei M, Wang Y, Ismail AM, Wissuwa M. Biochemical factors conferring shoot tolerance to oxidative stress in rice grown in low zinc soil. Funct Plant Biol. 2010;37(1):74–84.

    CAS  Article  Google Scholar 

  72. 72.

    Lee JS, Wissuwa M, Zamora OB, Ismail AM. Biochemical indicators of root damage in rice (Oryza sativa) genotypes under zinc deficiency stress. J Plant Res. 2017;130(6):1071–7.

    CAS  PubMed  Article  Google Scholar 

  73. 73.

    Sunkar R, Kapoor A, Zhu JK. Posttranscriptional induction of two cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell. 2006;18(8):2051–65.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  74. 74.

    Zeng HQ, Liu G, Kinoshita T, Zhang RP, Zhu YY, Shen QR, Xu GH. Stimulation of phosphorus uptake by ammonium nutrition involves plasma membrane H+ ATPase in rice roots. Plant Soil. 2012;357(1–2):205–14.

    CAS  Article  Google Scholar 

  75. 75.

    Yoshida S, Forno DA, Cock JH. Laboratory manual for physiological studies of rice. Manila: International Rice Research Institute; 1971.

  76. 76.

    Liu XS, Feng SJ, Zhang BQ, Wang MQ, Cao HW, Rono JK, Chen X, Yang ZM. OsZIP1 functions as a metal efflux transporter limiting excess zinc, copper and cadmium accumulation in rice. BMC Plant Biol. 2019;19(1):283.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  77. 77.

    Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  78. 78.

    Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. Gene Ontol Cons Nat Genet. 2000;25(1):25–9.

    CAS  Google Scholar 

  79. 79.

    Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38(Web Server issue):W64–70.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  80. 80.

    Axtell MJ, Meyers BC. Revisiting criteria for plant MicroRNA annotation in the era of big data. Plant Cell. 2018;30(2):272–84.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  81. 81.

    Zeng H, Zhang Y, Zhang X, Pi E, Zhu Y. Analysis of EF-hand proteins in soybean genome suggests their potential roles in environmental and nutritional stress signaling. Front Plant Sci. 2017;8(877):877.

    PubMed  PubMed Central  Article  Google Scholar 

  82. 82.

    Zeng HQ, Zhu YY, Huang SQ, Yang ZM. Analysis of phosphorus-deficient responsive miRNAs and cis-elements from soybean (Glycine max L.). J Plant Physiol. 2010;167(15):1289–97.

    CAS  PubMed  Article  Google Scholar 

  83. 83.

    Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29(9):e45.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

Download references

Acknowledgements

We thank Dr. Jin Xu (Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences) for critical reading of this manuscript.

Funding

This work was supported by the Zhejiang Provincial Natural Science Foundation of China (No. LY20C150002) and the National Natural Science Foundation of China (No. 31471937 and No. 31800201). The funders have no role in the study design, data analysis and interpretation, and manuscript writing, but just provide the financial supports.

Author information

Author notes

  1. Houqing Zeng and Xin Zhang contributed equally to this work.

Affiliations

  1. College of Life and Environmental Sciences, Hangzhou Normal University, Hangzhou, 311121, China

    Houqing Zeng

  2. College of Agriculture and Biotechnology, Hunan University of Humanities, Science and Technology, Loudi, 417000, China

    Xin Zhang

  3. College of Resources and Environmental Sciences, Nanjing Agricultural University, Nanjing, 210095, China

    Ming Ding & Yiyong Zhu

Authors
  1. Houqing Zeng
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xin Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ming Ding
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yiyong Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

H.Z. and Y.Z. conceived the study and designed the experiments. H.Z., X.Z., and M.D. performed the experiments. H.Z. and Y.Z. analyzed the data and wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Houqing Zeng.

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.

Primers used for quantitative RT-PCR analysis. Table S2. Summary of the reads from RNA sequencing in the 18 libraries. Table S3. RNA-seq analysis of rice shoot and root samples under Zn deprivation and Zn resupply. Table S4. List of differentially expressed loci in response to 14 d of Zn deprivation in rice shoots. Table S5. List of differentially expressed loci in response to 14 d of Zn deprivation in rice roots. Table S6. List of differentially expressed loci in response to 3 d of Zn resupply in rice shoots. Table S7. List of differentially expressed loci in response to 3 d of Zn resupply in rice roots. Table S8. GO representation of the over-represented GO terms in the DEGs regulated by Zn deprivation and/or Zn resupply. Table S9. Overview of small RNA sequencing data of the 18 libraries. Table S10. Counts and length distribution of total sRNAs and unique sRNAs in this study. Table S11. Summary of the detected known and predicted miRNAs in this study. Table S12. Differentially expressed miRNAs in response to 14 d of Zn deprivation and/or 3 d of Zn resupply in the roots and shoots. Table S13. Potential target genes of the Zn-responsive miRNAs. Table S14. Common DEGs identified in this study and earlier studies.

Additional file 2: Figure S1.

Pearson’s correlation (R-value) of three biological replicates of shoot and root samples for RNA sequencing under control, Zn deficiency, and Zn resupply conditions. Zn plus shoot, ZPS; Zn minus shoot, ZMS; Zn plus root, ZPR; Zn minus root, ZMR; Zn resupply shoot, ZRS; Zn resupply root, ZRR. Figure S2. Length distribution of the small RNAs from the small RNA sequencing data in the six samples. Each bar indicates the mean ± SD of three biological replicates. Figure S3. Pearson’s correlation (R-value) of three biological replicates of shoot and root samples for small RNA high-throughput sequencing under control, Zn deficiency, and Zn resupply conditions. Figure S4. Venn diagram representing the overlap of the Zn-deficiency- and Zn-resupply- responsive miRNAs (a, known and novel miRNAs; b, novel miRNAs) in roots and shoots. Figure S5. Predicted secondary structures of 38 Zn-responsive novel miRNAs identified in this study. The red lines indicate the mature miRNA sequences that were found to be differentially expressed under Zn deficiency and/or Zn resupply. The blue lines indicate the miRNA* sequences that were detected by the small RNA sequencing. If no blue line was indicated, it means that no miRNA* sequences were detected by the small RNA sequencing. Figure S6. Gene ontology (GO) representation of the overrepresented GO terms of biological processes in the potential target genes of the differentially expressed miRNAs under Zn deficiency and/or Zn resupply. The GO representation was generated using single enrichment analysis (SEA) tool on AgriGO (http://bioinfo.cau.edu.cn/agriGO/) (Fisher’s test, P < 0.05, FDR < 0.05). The number in parenthesis represents the FDR value. Figure S7. Gene ontology (GO) representation of the overrepresented GO terms of cellular component in the potential target genes of the differentially expressed miRNAs under Zn deficiency and/or Zn resupply. The GO representation was generated using single enrichment analysis (SEA) tool on AgriGO (http://bioinfo.cau.edu.cn/agriGO/) (Fisher’s test, P < 0.05, FDR < 0.05). The number in parenthesis represents the FDR value. Figure S8. Gene ontology (GO) representation of the overrepresented GO terms of biological process in the DEGs recovered by Zn resupply in roots and shoots. The GO representation was generated using single enrichment analysis (SEA) tool on AgriGO (http://bioinfo.cau.edu.cn/agriGO/) (Fisher’s test, P < 0.05, FDR < 0.05). The number in parenthesis represents the FDR value.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Zeng, H., Zhang, X., Ding, M. et al. Integrated analyses of miRNAome and transcriptome reveal zinc deficiency responses in rice seedlings. BMC Plant Biol 19, 585 (2019). https://doi.org/10.1186/s12870-019-2203-2

Download citation

Keywords

  • Rice (Oryza sativa)
  • Zinc deficiency
  • microRNA
  • Transcriptome
  • Copper
  • Oxidative stress

BMC Plant Biology

ISSN: 1471-2229

Contact us