Genome-wide Identification of LncRNAs
To systematically identify and characterize lncRNAs in rice, ssRNA sequencing (ssRNA-seq) was performed on shoot and root samples from rice seedlings grown in Fe-sufficient and -deficient conditions. After 10 days of Fe-deficient growth, rice plants showed significant chlorosis and lower chlorophyll content in the young leaves (Fig. 1a and 1b). The expression of typical Fe-deficiency responsive genes, such as the iron-related bHLH transcription factor 2 (IRO2), nicotianamine synthases 1 and 2 (NAS1/2), Fe(III)-DMA transporters (YSL15/16) and Iron-Regulated Transporter 1 (IRT1), were significantly increased (Fig. 1c), indicating that the rice seedlings were under iron deficiency at the sampling time.
The pipeline for lncRNA identification and characterization is shown in Figure 2a (see methods). Using this pipeline, approximately 700 million 150-bp, pair-end reads were assembled into 31,947 transcripts using Cufflinks. The Coding Potential Calculator (CPC) was used to evaluate the protein-coding potential of the transcripts to distinguish protein coding transcripts and lncRNAs. Transcripts more than 200 bp in length with CPC scores <0 were defined as lncRNAs, the remaining transcripts were classified as protein-coding transcripts (mRNAs). Using this method, 25,470 mRNAs and 6,477 lncRNAs were identified. Based on their relative position to protein-coding genes, lncRNAs can be classified into three types: Intergenic lncRNAs have no overlap with any protein-coding sequences, while sense lncRNAs and anti-sense lncRNAs overlap with one or more exons of another transcript on the same or opposite DNA strand, respectively [20]. Among the 6,477 lncRNAs identified in this work, 3,730 (58%) were intergenic lncRNAs, 1,696 (26%) were cis-lncRNAs, and 1,051 (16%) were antisense lncRNAs (Fig. 2b).
Fe-deficiency responsive LncRNAs and mRNAs in rice shoot and root
To identify the lncRNAs and mRNAs that are differentially expressed in response to Fe deficiency, the normalized expression levels (in fragments per kilobase of exon per million fragments mapped, FPKM) of the lncRNAs or the mRNAs were compared between the Fe-deficient and Fe-sufficient treatments. In shoots, 80 DE-lncRNAs were identified. Among them, 47 lncRNAs were up-regulated and 33 were down-regulated under Fe deficiency (Fig. 2c; Table S1a). In roots, 89 lncRNAs were up-regulated and 32 were down-regulated under Fe deficiency (Fig. 2c; Table S1b). In addition, 394 and 841 mRNAs were differentially expressed in either roots or shoots due to Fe deficiency, respectively. In shoots, 240 mRNAs were up-regulated and 154 were down-regulated (Fig. 2c; Table S1c), while in roots, 536 mRNAs were up-regulated and 305 mRNAs were down-regulated (Fig. 2c; Table S1d).
The DE-lncRNAs and -mRNAs were used to generate a heat map (Fig. 3). Classes I and III contained lncRNA and mRNA transcripts that were expressed significantly higher in Fe-sufficient than in Fe-deficient conditions in either roots (Class I) or shoots (Class III), respectively. In contrast, transcripts in Classes Ⅱ and Ⅳ had higher expression under Fe-deficient conditions in roots or shoots, respectively. Transcripts in Class Ⅴ were more highly expressed in both shoots and roots under Fe-deficient conditions. Among the five groups, Class Ⅱ (the transcripts induced in Fe-deficient roots) contained the largest number of both lncRNAs (Fig. 3a) and mRNAs (Fig. 3b). In total, 171 lncRNAs and 1,001 mRNAs were differentially expressed under the different Fe supply conditions (Fig. 3; Table S2).
Verification of LncRNAs responding to Fe deficiency using RT-qPCR
Quantitative real-time PCR (RT-qPCR) was performed to verify the accuracy of the RNA-seq data for the lncRNAs. Nine intergenic lncRNAs responding to Fe-deficiency were picked for the verification. The RT-qPCR results showed that lncRNAs XLOC_006153 and XLOC_028199 from Class IV were induced in shoots but not detected in roots regardless of the Fe supply status. LncRNAs XLOC_052823 and XLOC_007199 from Class II were up-regulated by Fe deficiency in the roots. The remaining 5 lncRNAs belonged to Class Ⅴ, which were induced upon Fe deficiency in both shoots and roots (Fig. 4). The strong correlation between the RNA-Seq and RT-qPCR result indicated the reliability of our transcriptomic profiling data.
Identification of LncRNAs as potential miRNA precursors and miRNA target mimics
MicroRNAs (miRNAs) regulate key aspects of development, cell signaling, and responses to various biotic and abiotic stresses via binding to specific complementary transcripts, including protein coding or non-coding sequences, resulting in the degradation or translational repression of the target. LncRNAs have been shown to function as precursors of miRNA in many studies [17, 29]. By aligning miRNA precursors to the 171 DE-lncRNAs, 2 of the lncRNAs, XLOC_010112 and XLOC_053944, were identified as potential miRNA precursors, namely of miR398a and miR164f, respectively (Fig. 5a). XLOC_010112 is located in a region that overlaps with a coding gene (LOC_Os10g18150) that so far seems to not be expressed (Fig. 5a). Under Fe deficiency, XLOC_010112 was down-regulated in shoot and up-regulated in root (Fig. 5b). XLOC_053944 is an intergenic lncRNA (Fig. 5a), specifically expressed in root and significantly induced by Fe starvation (Fig. 5b). Interestingly, miR398a and miR164f are both involved in regulation of Fe homeostasis in Arabidopsis [30, 31]. The predicted target genes of miR398a and miR164f include LOC_Os06g23650, LOC_Os06g46270 (ONAC11; OsY37; OMTN4), LOC_Os12g41680 (OMTN3; OsNAC60), and LOC_Os07g11360 (RAL3) (Table 1). Among them, OsNAC11 was significantly induced, while OsNAC60 was suppressed in root under Fe deficiency (Fig. 5c). We speculated that XLOC_053944 is induced in rice root under Fe-deficiency, and consequently generates miR164f. The increased amount of miR164f would reduce the transcript abundance of OsNAC11 and OsNAC60, which regulate Fe-related genes as transcription factors (Fig. 5d).
In addition to generating miRNAs, lncRNAs are also targets of miRNAs. In this case, lncRNAs function as target mimics of the sequestered transcript, known as an endogenous target mimic (eTM), to inhibit miRNA activity [26]. In order to further verify whether miRNA target mimicry is involved in Fe regulation in rice, the potential interactions between the Fe-responsive lncRNAs and known Fe-related microRNAs were investigated. Two endogenous target mimics (eTMs), eTM159 and eTM408, were identified. The lncRNAs XLOC_012715 (up-regulated in shoot and down-regulated in root under iron deficiency) and XLOC_054182 (only expressed in shoot, and slightly induced by Fe starvation), were predicted to bind miR159 and miR408, respectively (Fig. 5e-f). The potential target genes of miR159 / miR408 are listed in Table 1 and include the MYB transcription factors OsGAMYB and OsGAMYBL1 and the calmodulin-like protein OsCML27. The results demonstrated that target mimicry might be a part of the regulation of Fe homeostasis.
Interactions of DE-LncRNAs with mRNAs
Recent studies have shown that lncRNAs regulate the expression of protein-coding genes in two ways, those that are encoded nearby the coding genes (cis-regulation) on the same chromosome and those that are encoded elsewhere (trans-regulation) [32]. The genomic locations of the DE-lncRNAs and DE-mRNAs were mapped to each chromosome of the rice genome. The results indicated that both DE-lncRNAs and DE-mRNAs were evenly distributed to each chromosome, other than two regions on chromosome 9 and 12 that showed higher degrees of clustering of DE-lncRNAs. It is interesting that there were also many Fe-related DE-mRNAs that mapped near the region on chromosome 9 (Fig. 6a). For cis-target analysis, 12 DE-mRNAs spaced less than 10 kb away from 15 DE-lncRNAs in shoot (Table S3a) and less than 10 kb away from 37 DE-lncRNAs in root (Table S3b). The coding genes nearby Fe-regulated/responsive lncRNAs included bHLH transcription factors, E3 ubiquitin ligases and tyrosine protein kinases. Interestingly, 3 lncRNAs were found to be located nearby OsbHLH156 and OsHRZ2 (Oryza sativa haemerythrin motif-containing really interesting new gene (RING)- and zinc-finger protein 2), which are two important regulators involved in Fe homeostasis in rice (Fig. 6b) [14, 33]. LncRNA XLOC_037283 and XLOC_034336 are located in or nearby OsbHLH156. XLOC_037283 is likely a natural antisense transcript (NAT) located near part of promoter and coding region of OsbHLH156 in the opposite orientation, while XLOC_034336 is within 8 kbs upstream of the start codon of OsbHLH156. Both lncRNAs showed similar expression patterns to OsbHLH156 under Fe deficiency (Fig. 6c). XLOC_043545 is within 8 kbs upstream of the start codon of OsHRZ2 and was mainly expressed in rice root under iron deficiency, while OsHRZ2 was induced in both the shoot and root (Fig. 6d). The results demonstrated that lncRNAs might play roles in the Fe signaling pathway as cis-regulators and that they are likely involved in transcriptional or post-transcriptional regulation of OsbHLH156 and OsHRZ2.
For trans-target analysis, 478 interaction nodes between DE-lncRNAs and DE-mRNAs in shoot and 1,516 in root were inferred according to the complementary pairing of bases (Table S4). GO enrichment analysis was performed to identify the potential functions of the trans-target genes. As shown in Figure 7a, we found 14 GO terms that were significantly enriched in root, but only 2 in shoot. Among them, “Response to iron ion”, “Metal ion transport”, “Iron ion transport” and “Iron ion homeostasis” were all associated with response to Fe-deficient stress. Interaction nodes among the DE-lncRNAs and Fe-related genes were built into interaction networks. There were 6 DE-lncRNAs that were predicted to interact with more than 5 Fe-related genes (Fig. 7b). The results implied that a complex regulation network between lncRNAs and mRNAs might contribute to Fe homeostasis regulation in rice.
DE-LncRNAs were transcriptionally regulated by the transcription factors bHLH156 and IRO2
To test whether expression of DE-LncRNAs could be regulated by Fe-related transcription factors, an ssRNA-seq was performed on shoot and root samples of a knock-out mutant of bHLH156 grown in Fe-deficient and -sufficient conditions. bHLH156 acts as a core transcription factor in regulating Fe homeostasis together with IRO2 [14]. The number of DE-lncRNAs in bhlh156 shoot (145 up-regulated and 89 down-regulated) and bhlh156 root (419 up-regulated and 177 down-regulated) was greater than that in WT under normal conditions (Fig. 8a). Under Fe deficiency, 495 lncRNAs were up-regulated and 168 were down-regulated in bhlh156 root when compared with WT (Fig. 8a). The lncRNAs responding in an antagonistic manner in rice roots under Fe-deficiency condition are most likely regulated by bHLH156. In comparison to WT, Fe-deficiency-induced lncRNAs in shoots (14) and in roots (50) were suppressed in the bhlh156 mutant. Conversely, lncRNAs down-regulated in response to Fe deficiency in shoots (3) and roots (12) of wild type showed significantly higher expression in bhlh156 (Figure 8b and 8c). To verify that these lncRNAs were truly regulated by bHLH156, four (XLOC_011962, XLOC_018668, XLOC_043504 and XLOC_056321) were chosen for analysis by RT-qPCR in both the bhlh156 and iro2 mutants. IRO2 is a necessary interacting partner for bHLH156 to activate downstream genes [14]. XLOC_011962, XLOC_018668, XLOC_043504 and XLOC_056321 were all specifically expressed in root and dramatically induced under Fe-deficiency in WT, but were barely detectable in either the bhlh156 or iro2 mutants (Fig. 8d). The results demonstrated that a number of DE-LncRNAs could be regulated by bHLH156 and IRO2 at the transcriptional level.