Nội dung

Zeng et al. BMC Plant Biology (2019) 19:585 https://doi.org/10.1186/s12870-019-2203-2 RESEARCH ARTICLE Open Access Integrated analyses of miRNAome and transcriptome reveal zinc deficiency responses in rice seedlings Houqing Zeng1*† , Xin Zhang2†, Ming Ding3 and Yiyong Zhu3 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. Keywords: Rice (Oryza sativa), Zinc deficiency, microRNA, Transcriptome, Copper, Oxidative stress 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 Znfinger 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 * Correspondence: zenghq@hznu.edu.cn † Houqing Zeng and Xin Zhang contributed equally to this work. 1 College of Life and Environmental Sciences, Hangzhou Normal University, Hangzhou 311121, China Full list of author information is available at the end of the article 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–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–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 © The Author(s). 2019 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. Zeng et al. BMC Plant Biology (2019) 19:585 [6, 14, 15]; expanding the generation of crown roots [16]; involvement of transporters such as ZRT, IRTrelated proteins (ZIPs), metal tolerance proteins and heavy metal tolerance family proteins in Zn uptake and translocation [17–21]; adjusting the expression of Znrequiring enzymes [2, 13]. Recently, transcriptomic profilings were analyzed to identify Zn-responsive genes in plants [22–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–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–37]. miRNAs also regulate plant responses to multiple nutrient deprivation stresses and are involved in nutrient homeostasis regulation [38–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–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. Page 2 of 18 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. 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, Os06g0566 201, 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) Zeng et al. BMC Plant Biology (2019) 19:585 Page 3 of 18 Fig. 1 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) 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). 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). Zeng et al. BMC Plant Biology (2019) 19:585 Page 4 of 18 Fig. 2 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 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 miRNAmediated 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; osamiR818a-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. 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 Zeng et al. BMC Plant Biology (2019) 19:585 Page 5 of 18 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 Gene Description ZPS (FPKM) ZMS (FPKM) ZRS (FPKM) Os02g0192700 Similar to Thioredoxin peroxidase 0 68.55 ± 60.14 0 Os02g0663100 GRAS transcription factor domain containing protein 0 22.24 ± 4.50 0 Os05g0102000 SAM dependent carboxyl methyltransferase family protein 0 19.82 ± 17.40 0 Os10g0493300 Hypothetical conserved gene 0 13.02 ± 13.99 0 Os04g0382900 Conserved hypothetical protein 0 6.38 ± 9.55 0 Os01g0504100 Protein of unknown function DUF250 domain containing protein 0 6.10 ± 5.29 0 Os10g0528400 Glutathione S-transferase, C-terminal-like domain containing protein 0 5.39 ± 4.72 0 Os12g0103500 Ribosome-inactivating protein domain containing protein 0 4.59 ± 4.11 0 Os09g0116600 Hypothetical gene 0 4.06 ± 5.68 0 Os11g0490300 Hypothetical conserved gene 0 3.26 ± 2.85 0 Os03g0746700 Conserved hypothetical protein 0 3.08 ± 2.72 0 Os01g0532300 Conserved hypothetical protein 0 2.28 ± 2.08 0 Os06g0208951 Conserved hypothetical protein 0 2.00 ± 1.81 0 Os05g0562800 Similar to predicted protein 0 1.90 ± 2.55 0 Os03g0272900 DVL family protein 0 1.69 ± 1.59 0 Os09g0444800 Protein of unknown function DUF716 family protein 0 1.57 ± 1.36 0 Os10g0328600 – 0 1.39 ± 0.88 0 Os02g0550400 Tetratricopeptide-like helical domain containing protein 0 1.10 ± 0.95 0 Os07g0586900 GRAS transcription factor domain containing protein 0 1.02 ± 1.00 0 Os09g0412900 Pentatricopeptide repeat domain containing protein 0 0.64 ± 0.64 0 Os11g0606200 Similar to Leucine Rich Repeat family protein 0 0.49 ± 0.49 0 Os10g0328600 Hypothetical protein 0.02 ± 0.01 8.27 ± 1.81 0.18 ± 0.09 Os04g0280500 Similar to Non-S-locus F-box-like protein 2011 0.03 ± 0.03 10.03 ± 1.61 0.49 ± 0.22 Os04g0561500 Similar to Prolyl endopeptidase (Post-proline cleaving enzyme) 1.78 ± 0.63 204.72 ± 11.10 18.07 ± 8.12 Os06g0566201 – 1.38 ± 0.28 20.93 ± 10.37 1.48 ± 0.19 Os02g0306401 Similar to Nicotianamine aminotransferase A 5.81 ± 1.28 76.61 ± 8.44 15.04 ± 5.84 Os06g0566300 Zinc transporter OsZIP10 2.97 ± 0.26 38.94 ± 10.62 1.84 ± 0.34 MSTRG.18312 – 0.18 ± 0.02 1.92 ± 0.42 0.17 ± 0.08 Os08g0207500 Zinc transporter OsZIP4 1.7 ± 0.42 16.32 ± 8.27 1.9 ± 0.39 Os07g0232800 Zinc transporter OsZIP8 6.84 ± 1.31 50.73 ± 5.05 11.81 ± 6.25 Os03g0293100 Conserved hypothetical protein 7.25 ± 0.67 44.01 ± 14.47 7.78 ± 2.48 Os09g0485900 Similar to 60S ribosomal protein L9 (Gibberellin-regulated protein) 35.05 ± 6.96 164.62 ± 105.32 29.16 ± 13.21 Os11g0577866 Conserved hypothetical protein 9.27 ± 2.32 0 7.93 ± 4.44 Os06g0715200 Conserved hypothetical protein 8.48 ± 7.9 0 15.11 ± 6.49 Os08g0526100 NAD(P)-binding domain containing protein 6.15 ± 5.96 0 4.27 ± 3.76 Os01g0243450 Non-protein coding transcript 3.96 ± 0.68 0 1.66 ± 1.54 Os05g0528701 Non-protein coding transcript 1.81 ± 1.59 0 1.36 ± 1.27 Os01g0856800 Pleckstrin homology-type domain containing protein 0.86 ± 0.8 0 1.14 ± 0.98 Os09g0460300 Alpha/beta hydrolase fold-3 domain containing protein 0.8 ± 0.74 0 9.12 ± 8.68 Os08g0451201 – 0.69 ± 0.67 0 0.85 ± 0.76 Os03g0411900 – 0.64 ± 0.73 0 0.35 ± 0.32 Os01g0880250 Hypothetical protein 6.77 ± 6.57 0.03 ± 0.02 3.64 ± 5.01 Zeng et al. BMC Plant Biology (2019) 19:585 Page 6 of 18 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 (Continued) Gene Description ZPS (FPKM) ZMS (FPKM) ZRS (FPKM) MSTRG.5686 – 78.84 ± 109.73 0.42 ± 0.31 227.64 ± 346.33 Os06g0639800 Cytochrome P450 family protein 5.69 ± 0.85 0.06 ± 0.10 5.2 ± 2.26 MSTRG.5688 – 7.28 ± 9.58 0.20 ± 0.06 23.02 ± 34.31 Os03g0592500 – 66.95 ± 83.98 5.22 ± 1.37 30.1 ± 25.42 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–49]; most of them were Cu-containing proteins and were involved in Cu homeostasis, reactive oxygen species (ROS) homeostasis, photosynthesis and stress tolerance [49–58]. 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). 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. Confirmation of miRNAome and transcriptome results by qRT-PCR To confirm the small RNA sequencing results, seven Znresponsive 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, primiR398a was up-regulated by Zn deficiency and downregulated 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, Zeng et al. BMC Plant Biology (2019) 19:585 Page 7 of 18 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 Gene Description ZPR (FPKM) ZMR (FPKM) ZRR (FPKM) Os08g0235800 Similar to WRKY transcription factor 25 0 23.49 ± 22.4 0 Os01g0896400 Similar to JHL23J11.5 protein 0 11.66 ± 18.26 0 Os01g0693800 Similar to Threonine synthase, chloroplast precursor 0 11.17 ± 9.72 0 Os05g0121500 Similar to structural constituent of ribosome 0 6.57 ± 6.92 0 Os03g0191900 Pathogenesis-related transcriptional factor and ERF domain containing protein 0 5.76 ± 5.66 0 Os11g0107600 Similar to prenylated Rab receptor 2 0 5.35 ± 4.77 0 Os02g0455400 Hypothetical conserved gene 0 4.66 ± 4.04 0 Os03g0852400 Conserved hypothetical protein 0 4.28 ± 3.71 0 Os03g0280750 Plant disease resistance response protein family protein 0 3.48 ± 3.11 0 Os05g0571600 Conserved hypothetical protein 0 3.2 ± 2.77 0 Os02g0310400 Conserved hypothetical protein 0 3.08 ± 2.68 0 Os05g0525900 Similar to Zing finger transcription factor PEI1 0 2.67 ± 2.32 0 Os01g0742300 3-hydroxyacid dehydrogenase/reductase domain containing protein 0 1.96 ± 1.79 0 Os03g0253200 Similar to WAG1; kinase 0 1.91 ± 1.66 0 Os06g0192800 Similar to RING-H2 finger protein ATL1R (RING-H2 finger protein ATL8) 0 1.59 ± 1.78 0 Os01g0234700 Harpin-induced 1 domain containing protein 0 1.25 ± 1.4 0 Os03g0103400 GRAS transcription factor domain containing protein 0 0.99 ± 0.89 0 Os10g0490100 Virulence factor, pectin lyase fold family protein 0 0.98 ± 0.88 0 Os03g0245700 Similar to senescence-associated protein 15 0 0.96 ± 0.84 0 Os05g0578100 Protein of unknown function DUF1645 family protein 0 0.82 ± 0.15 0 Os01g0584300 – 0 2.42 ± 2.23 0 Os04g0685700 Similar to H0723C07.4 protein 0 0.75 ± 0.65 0 Os07g0585900 Hypothetical conserved gene 0 0.75 ± 0.66 0 MSTRG.6046 – 0 0.66 ± 0.7 0 Os02g0242900 Similar to hydroquinone glucosyltransferase 0 0.63 ± 0.6 0 Os11g0256100 Pentatricopeptide repeat domain containing protein 0 0.41 ± 0.35 0 Os12g0137700 Sulfotransferase family protein 0.01 ± 0.01 3.9 ± 3.66 0.03 ± 0.04 Os04g0280500 Similar to Non-S-locus F-box-like protein 2011 0.07 ± 0.02 18.01 ± 5.28 0.31 ± 0.22 Os10g0553600 Exostosin-like family protein 0.01 ± 0.01 1.69 ± 0.89 0.01 ± 0.01 Os08g0207401 Hypothetical gene 0.01 ± 0.03 0.97 ± 0.34 0.04 ± 0.03 Os10g0328600 Hypothetical protein 0.11 ± 0.09 6.52 ± 1.13 0.13 ± 0.03 Os06g0566201 – 0.4 ± 0.03 13.82 ± 8.44 0.88 ± 0.15 Os06g0566300 Zinc transporter OsZIP10 0.88 ± 0.12 23.97 ± 7.22 1.58 ± 0.07 Os04g0561500 Similar to Prolyl endopeptidase (Post-proline cleaving enzyme) 7.38 ± 0.4 155.99 ± 10.59 10.99 ± 4.64 Os05g0472400 Similar to Zinc transporter 9 1.48 ± 0.25 29.79 ± 5.17 5.12 ± 2.05 Os07g0232800 Zinc transporter OsZIP8 0.75 ± 0.23 14.84 ± 1.15 0.99 ± 0.31 Os04g0304400 Similar to MADS-box protein AGL16-II 0.51 ± 0.04 8.26 ± 3.48 1.83 ± 0.09 Os08g0207500 Zinc transporter OsZIP4 2.9 ± 0.24 28.12 ± 9.16 4.29 ± 1.91 MSTRG.22477 – 0.34 ± 0.11 3.12 ± 0.97 0.09 ± 0.06 MSTRG.30070 – 21.6 ± 19.8 91.8 ± 29.43 0.46 ± 0.24 Os06g0131300 Arginine decarboxylase, Chilling stress respons 98.89 ± 73.37 0 123.68 ± 112.96 Os03g0244950 Heat shock protein DnaJ, N-terminal domain containing protein 6.6 ± 7.68 0 3.45 ± 3.04 Os06g0505501 Non-protein coding transcript 3.46 ± 1.52 0 1.61 ± 1.6 Zeng et al. BMC Plant Biology (2019) 19:585 Page 8 of 18 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 (Continued) Gene Description ZPR (FPKM) ZMR (FPKM) ZRR (FPKM) Os02g0756850 Hypothetical protein 2.36 ± 0.31 0 1.33 ± 1.33 Os11g0706200 GRAS transcription factor domain containing protein 1.8 ± 1.63 0 1.98 ± 1.91 Os12g0124700 CDC45-like protein family protein 1.16 ± 1.08 0 1.83 ± 2.03 Os09g0412900 Pentatricopeptide repeat domain containing protein 1.12 ± 0.98 0 0.96 ± 0.84 Os03g0140300 Conserved hypothetical protein 0.99 ± 1.07 0 0.62 ± 0.54 Os04g0277400 Membrane bound O-acyl transferase, MBOAT family protein 0.84 ± 0.71 0 1.14 ± 1.09 Os10g0328700 Tetratricopeptide-like helical domain containing protein 0.58 ± 0.47 0 0.39 ± 0.34 Os07g0172600 Pentatricopeptide repeat domain containing protein 0.55 ± 0.58 0 0.67 ± 0.63 Os01g0163000 Leucine-rich repeat, N-terminal domain containing protein 0.46 ± 0.1 0 0.4 ± 0.38 Os03g0821700 Pentatricopeptide repeat domain containing protein 0.45 ± 0.15 0 0.53 ± 0.19 Os02g0621600 Similar to OSIGBa0157K09-H0214G12.24 protein 0.44 ± 0.59 0 0.62 ± 0.4 Os11g0133500 Serine/threonine protein kinase-related domain containing protein 0.35 ± 0.31 0 0.71 ± 0.8 Os12g0585300 Cyclin-like F-box domain containing protein 1.17 ± 1 0 1.53 ± 1.33 Os11g0111000 Conserved hypothetical protein 1.1 ± 0.92 0.01 ± 0.01 1.07 ± 1.11 Os03g0758500 Similar to blue copper protein 6.31 ± 5.39 0.15 ± 0.1 11.47 ± 3.69 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. 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 miR397b5p; 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 upregulated 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). 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 Znresponsive miRNAs, and half of them showed a differential expression with miRNAs under Zn deficiency and/or Zn resupply (Table 4), suggesting the interactions Zeng et al. BMC Plant Biology (2019) 19:585 Page 9 of 18 Table 3 The validated target genes of some Zn-responsive miRNAs miRNA Target gene Target gene annotation osamiR398a Os07g0665200 Cu/Znsuperoxidase dismutase2/CSD2 Potential function Stress tolerance, antioxidant defense, copper homeostasis Os03g0351500 Cu/Znsuperoxidase dismutase1/CSD1 Os11g0203300 Superoxide dismutaseX/SODX Os04g0573200 Copper chaperone for superoxide dismutase/CCSD osamiR4083p Os08g0482700 Uclacyanin-like protein UCL30 Os03g0709100 Uclacyanin-like protein UCL8 Stress tolerance, photosynthesis, copper homeostasis, Os04g0545600 Uclacyanin-like protein UCL15 Os06g0218600 Uclacyanin-like protein UCL16 Os02g0653200 Plantacyanin like 2/Uclacyanin-like protein UCL4 Os02g0731400 Plantacyanin-like 3/Uclacyanin-like protein UCL5 Os02g0758800 Plantacyanin-like 4/Uclacyanin-like protein UCL6 Os01g0827300 Laccase 3/LAC3 Os03g0297900 Laccase 12/LAC12 osamiR5285p Os06g0567900 Ascorbate oxidase 2 Os08g0137400 Uclacyanin-like protein UCL23 Os07g0570550 Cupredoxin domain containing protein Reactive oxygen species homeostasis, stress tolerance, copper homeostasis, strigolactone signal perception Os06g0154200 F-box/LRR-repeat MAX2 homolog/ DWARF3 Fig. 3 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. 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 Os08g0561700 Superoxide dismutase 4/CSD4 Os05g0458600 Laccase 15/LAC15 osamiR397b5p Copper homeostasis, lignin biosynthesis, brassinosteroid response 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 Zeng et al. BMC Plant Biology (2019) 19:585 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–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 Page 10 of 18 Cu deficiency and their target genes encode various Cucontaining 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-resupplyinduced 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