Transcriptomic Proling of Fe-Responsive lncRNAs and Their Regulation in Rice

Background Iron (Fe) plays a vital role in various cellular processes in plants, including biosynthesis of chlorophyll, photosynthesis and respiration. Fe deciency directly affects crop growth and development, ultimately resulting in reduced crop yield and quality. Long non-coding RNAs (lncRNAs) have recently been demonstrated to play critical regulatory roles in a multitude of pathways across numerous species. However, systematic screening of lncRNAs responding to Fe deciency in plants has not been reported. Results In this work, lncRNAs responsive to Fe deciency were identied across the rice genome by strand-specic RNA sequencing. In total, 6,477 lncRNAs were identied. In Fe-decient conditions, 47 lncRNAs were up-regulated and 33 lncRNAs were down-regulated in shoots, while 89 lncRNAs were up-regulated and 32 lncRNAs were down-regulated in roots, compared to normal conditions. Two lncRNAs (XLOC_010112 and XLOC_053944) were identied as potential miRNA precursors and another two (XLOC_012715 and XLOC_054182) as miRNA target mimics that may participate in Fe regulation. A number of differentially expressed lncRNAs (DE-lncRNAs) are likely to modulate the expression of Fe-related genes via a cis- or trans-regulation mode, including 3 DE-lncRNAs (XLOC_034336, XLOC_037283 and XLOC_043545) located nearby OsbHLH156 and OsHRZ2 genomic regions. Seventy-six DE-lncRNAs were found to be regulated by bHLH156 at the transcriptional level. Conclusions This study provides a rst prole of lncRNA expression as well as identies the lncRNAs likely to play important roles in the regulation of Fe homeostasis. This identication and characterization form an important basis for understanding Fe regulatory networks in rice.


Background
Iron (Fe) is an essential micronutrient for plants but is often limited due to low availability in the soil [1].
To overcome Fe de ciency, plants have evolved two strategies to optimize Fe acquisition and uptake, i.e., the reduction strategy (Strategy I) for non-gramineous plants and the chelation strategy (Strategy II) for gramineous plants [2,3]. Rice, which is adapted to grow submerged in a paddy where the reduced form of Fe is available, utilizes both Fe uptake strategies I and II [4][5][6]. A large number of genes are known for their involvement in Fe uptake and homeostasis, including those that encode transcription factors for regulating expression of downstream Fe-responsive genes, enzymes for synthesis of phytosiderophores (MAs), and transporters of MA-Fe(III) or Fe(II) in rice [3,4,[7][8][9][10][11][12][13][14]. However, regulation of Fe-responses by long noncoding RNAs (lncRNAs) has not been reported.
LncRNAs are transcripts of more than 200 nucleotides in length but without coding potential that have recently gained widespread attention [15]. LncRNAs play roles in numerous crucial biological processes across many species by regulating expression of mRNAs at epigenetic, transcriptional, posttranscriptional, translational and post-translational levels [15][16][17]. LncRNAs are classi ed as sense, antisense, intronic, and intergenic, according to their position in relation to neighboring coding genes [18,19]. In plants, lncRNAs were reported to be involved in development and stress responses [20,21]. For instance, the lncRNAs COLDAIR (Cold-Assisted Intronic noncoding RNA) and COOLAIR (Cold-Induced Long Antisense Intragenic RNA) are both located in the Flowering Locus C(FLC) gene, which regulates owering time. COLDAIR and COOLAIR regulate expression of FLC at the epigenetic level by interacting with an evolutionarily conserved repressive complex, PRC2 (Polycomb Repressive Complex 2) [22][23][24].
Another lncRNA, the 1,236-nucleotide long LDMAR (Long-day-speci c male-fertility-associated RNA), regulates photoperiod-sensitive male sterility (PSMS) in rice [25]. A number of lncRNAs have been reported to regulate phosphate homeostasis. IPS1 (Induced by phosphate starvation 1) reduces phosphorus acquisition by inhibiting the activity of miR399, through the target mimicry mechanism [26].
The cis-natural antisense RNA (cis-NATPHO1;2), transcribed from OsPHO1;2, was found to be a translational enhancer of its sense gene (OsPHO1;2) [27]. In yeast, prt (pho1-repressing transcript), generated from the promoter region of the pho1 gene, regulates expression of pho1 responding to different phosphate levels [28]. These two studies show that a complicated network involving lncRNAs regulates phosphate homeostasis. In contrast, there are no reported lncRNAs responding to Fe de ciency.
In this study, the transcriptome of rice was surveyed to systematically identify and characterize any lncRNAs that respond to Fe de ciency.

Genome-wide Identi cation 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-su cient and -de cient conditions. After 10 days of Fe-de cient growth, rice plants showed signi cant chlorosis and lower chlorophyll content in the young leaves ( Fig. 1a and 1b). The expression of typical Fe-de ciency 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 signi cantly increased ( Fig. 1c), indicating that the rice seedlings were under iron de ciency at the sampling time.
The pipeline for lncRNA identi cation 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 Cu inks. 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 de ned as lncRNAs, the remaining transcripts were classi ed as proteincoding transcripts (mRNAs). Using this method, 25,470 mRNAs and 6,477 lncRNAs were identi ed. Based on their relative position to protein-coding genes, lncRNAs can be classi ed 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 identi ed 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-de ciency responsive LncRNAs and mRNAs in rice shoot and root To identify the lncRNAs and mRNAs that are differentially expressed in response to Fe de ciency, 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-de cient and Fe-su cient treatments. In shoots, 80 DE-lncRNAs were identi ed. Among them, 47 lncRNAs were up-regulated and 33 were downregulated under Fe de ciency ( Fig. 2c; Table S1a). In roots, 89 lncRNAs were up-regulated and 32 were down-regulated under Fe de ciency ( Fig. 2c; Table S1b). In addition, 394 and 841 mRNAs were differentially expressed in either roots or shoots due to Fe de ciency, respectively. In shoots, 240 mRNAs were up-regulated and 154 were down-regulated ( Fig. 2c; Table S1c), while in roots, 536 mRNAs were upregulated 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 signi cantly higher in Fe-su cient than in Fe-de cient conditions in either roots (Class I) or shoots (Class III), respectively. In contrast, transcripts in Classes and had higher expression under Fe-de cient conditions in roots or shoots, respectively. Transcripts in Class were more highly expressed in both shoots and roots under Fe-de cient conditions. Among the ve groups, Class (the transcripts induced in Fe-de cient 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).
Veri cation of LncRNAs responding to Fe de ciency 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-de ciency were picked for the veri cation. 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 de ciency in the roots. The remaining 5 lncRNAs belonged to Class , which were induced upon Fe de ciency in both shoots and roots (Fig. 4). The strong correlation between the RNA-Seq and RT-qPCR result indicated the reliability of our transcriptomic pro ling data.
Identi cation 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 speci c complementary transcripts, including protein coding or noncoding 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 identi ed 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 de ciency, XLOC_010112 was down-regulated in shoot and up-regulated in root (Fig. 5b).
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 identi ed. The lncRNAs XLOC_012715 (up-regulated in shoot and down-regulated in root under iron de ciency) 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 Ferelated 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-nger 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 de ciency (Fig. 6c). XLOC_043545 is within 8 kbs upstream of the start codon of OsHRZ2 and was mainly expressed in rice root under iron de ciency, 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 posttranscriptional 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 signi cantly 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-de cient 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 Fede cient and -su cient 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 de ciency, 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-de ciency condition are most likely regulated by bHLH156. In comparison to WT, Fe-de ciency-induced lncRNAs in shoots (14) and in roots (50) were suppressed in the bhlh156 mutant. Conversely, lncRNAs down-regulated in response to Fe de ciency in shoots (3) and roots (12) of wild type showed signi cantly 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 speci cally expressed in root and dramatically induced under Fe-de ciency 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.

Discussion
LncRNAs have roles in a wide range of biological processes, including development, stress responses, and plant nutrition. In this work, lncRNAs that respond to Fe de ciency in rice roots and shoots were identi ed. The results generated in this study promote our understanding of how rice plants respond to Fe de ciency.
LncRNAs arise from intergenic, intronic, or coding regions in the sense and antisense directions but at lower expression levels than mRNAs. Thus, identi cation of lncRNAs requires the use of an ssRNA-seq strategy. In this study, 6477 lncRNAs were identi ed and characterized (Fig. 2a). Differentially expressed lncRNAs and mRNAs were identi ed by comparing their expression levels between +Fe and -Fe conditions. The expression patterns divided the differentially expressed (DE) RNA molecules into ve classes (Fig. 3). Class , the molecules up-regulated in rice root with Fe de ciency, had the greatest number of transcripts. Among the differentially expressed RNA molecules, the number of lncRNAs and mRNAs responding to Fe de ciency showed a similar trend, with more lncRNAs and mRNAs up-regulated under Fe de ciency in roots. Moreover, a greater number of lncRNAs and mRNAs were detected in roots in Fe-de cient condition.
MiRNAs are small RNAs that regulate target genes at both the transcriptional and post-transcriptional levels and that are generated by sequential cleavage of long precursor transcripts. Some lncRNAs could also act as primary transcripts of miRNAs [29]. In this study, two miRNA precursors were identi ed as generated from lncRNAs, of which XLOC_053944 might produce miR164f, which degrades OsNAC mRNA ( Fig. 5a-d). The targeting of NAC by miR164 acts as a negative regulator of drought tolerance in rice [34]. In addition, NAC genes have been found to play important roles in Fe homeostasis [10]. Further investigation should be conducted to verify whether miR164-targeting of NAC also participates in Fe regulation, which might serve as a link between drought tolerance and response to Fe de ciency.
LncRNAs have been shown to regulate phosphate homeostasis in plants by a novel mechanism called target mimicry [26,35]. Two endogenous lncRNA target mimics (eTMs) were identi ed in rice, namely eTM159 and eTM408, which target two Fe-related miRNAs, miR159 and miR408 (Fig. 5e) [31]. Therefore, a target mimicry mechanism similar to the IPS1-mi399 regulation of phosphate homeostasis might also exist in Fe regulation. We speculated that iron de ciency in rice would prompt XLOC_012715 and XLOC_054182 to target miR159 and miR408, thus preventing them from degrading their target coding genes, including a zinc nger protein, calmodulin-like protein 27, OsGAMYB, and OsGAMYB (Table 1). Most target genes identi ed were not Fe-related, except for OsCML27 (LOC_Os03g21380) which was upregulated in root. It is possible that this target mimicry might only happen in speci c cell types, and any changes in expression levels of target genes could not be detected using whole shoot or root tissue samples.
LncRNAs have been shown to either regulate expression of adjacent genes via recruitment of regulatory complexes through RNA-protein interactions or to correlate with expression of neighboring genes through acting as local regulators [36,37]. In order to study whether this mechanism is involved in Fe regulation, we compared the locations of DE-lncRNAs to the locations of 64 known Fe-related genes (up to 10 kb upand down-stream) (Table S5). Three DE-lncRNAs were identi ed nearby OsbHLH156 and OsHRZ2 (Fig.  6b), which are two important regulators involved in Fe signaling. OsbHLH156 regulates Strategy II iron acquisition as a core transcription factor [14], and OsHRZ2 is a putative iron-binding sensor that negatively regulates iron acquisition under Fe su ciency [33]. Both OsbHLH156 and OsHRZ2 were strongly induced by Fe de ciency, however the transcriptional and post-transcriptional regulation of these two genes are largely unknown. Two of the DE-lncRNAs, XLOC_034336 and XLOC_043545, were located upstream of OsbHLH156 and OsHRZ2, respectively. XLOC_037283 is a NAT that overlapped within the genomic sequence of OsbHLH156 and showed a synchronous expression pattern with OsbHLH156 (Fig.  b-d). A series of case studies have shown that NATs can either positively or negatively regulate expression of a cognate loci [38]. For example, an Arabidopsis gene of unknown function named SRO5, which overlaps the P5CDH (Δ 1 -pyrroline-5-carboxylate dehydrogenase) gene in the antisense orientation, generates both 24-nt and 21-nt siRNAs that regulate P5CDH at the posttranscriptional level. The P5CDH-SRO5 gene pair de nes a mode of siRNA function that may be applied to other cis-antisense gene pairs [39]. For example, the Flowering Locus C (FLC) contains multiple coldinduced long antisense intragenic RNAs (COOLAIR) that are transcribed in the antisense orientation in relation to FLC. The COOLAIR antisense lncRNAs have an early role in the epigenetic silencing of the FLC gene, acting to silence FLC transcription transiently [22]. In mammals, NATs can increase mRNA stability by forming a duplex with the sense gene, similar to what occurs at the BACE1 locus [40]. To determine if a similar regulation mechanism exists between XLOC_037283 and OsbHLH156, a more detailed analysis of the spatio-temporal expression of XLOC_037283 and OsbHLH156 upon Fe de ciency and of the methylation level at the OsbHLH156 locus should be conducted. These DE-lncRNAs will enrich the further study of OsbHLH156 and OsHRZ2 under iron de cient condition and our understanding of the total Fe regulation network in rice.
The DE-LncRNAs described above mainly function as upstream regulators of Fe-related genes, however the transcriptional regulation of these Fe-de ciency responsive lncRNAs remains unknown. To study whether DE-lncRNAs could also be regulated by Fe-related genes, an additional ssRNA-seq was performed using the bhlh156 mutant. A total of 76 DE-lncRNAs (in shoot, root, or both) were identi ed whose expression might be activated or inhibited by the transcription factor bHLH156 (Fig. 8c), which, in conjucntion with IRO2, is required for induction of nearly all Strategy II iron acquisition genes in rice [14].
RT-PCR was performed to verify the expression of these DE-lncRNAs in the bhlh156 and iro2 mutants. The results indicate that bHLH156 and IRO2 might also regulate Fe homeostasis via activating downstream lncRNAs.

Conclusions
Our study provides insight into the potential functions and regulatory interactions of mRNA and lncRNA molecules when rice plants are grown in Fe-de cient conditions. We believe our study will serve as an initial reference for understanding the function of lncRNAs in regulating iron homeostasis in Oryza sativa and provides the genic identities needed to design the next wave of experiments aimed at understanding this additional layer of regulation.

Plant growth condition
Oryza sativa L. cv. Nipponbare (Nip) obtained from BIOGLE GeneTech (Hangzhou, China) was used in this study as the wild type (WT). The WT seeds were germinated in the dark for 3 days, and then placed on a net oating on a solution with or without iron (1.43 mM NH 4  Strand-speci c RNA library construction and sequencing Shoots and roots were separately collected from seedlings grown hydroponically for 10 days after germination with or without Fe, and immediately frozen in liquid nitrogen. Three biological replicates were used for each sample. Total RNA was extracted from these tissues using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Strand-speci c RNA (ssRNA) library construction and RNA sequencing were performed by the Beijing Genomics Institute in Shenzhen (BGI, Shenzhen, China). To construct the ssRNA library, rRNA was removed with the Ribo-Zero Gold rRNA Removal Kit (Epicentre, Madison, WI, USA) from the pooled RNA. The RNA was fragmented into 200 -500 nts in length using fragmentation buffer. After synthesis of rst-strand and second-strand cDNA, adapters were added to both sides of the short fragments. The second strand was degraded by Uracil-N-Glycosylase. The resulting single strand cDNA was PCR ampli ed and then sequenced by Illumina HiSeq PE151. The RNAseq data is available in the NCBI (Accession number: PRJNA527175).

Identi cation and characterization pipeline of LncRNAs
The raw data obtained by Illumina sequencing was ltered into clean data by removing adaptor sequences, low-quality reads and rRNA-containing reads with SOAPnuke and SOAP [41]. The dataset was aligned to the rice genome (Rice Genome Annotation Project) using the improved TopHat v 2.0 [42]. Cu inks was used to reconstruct the transcripts. After ltering background noise transcripts, the nal expression data was produced [43]. Transcripts shorter than 200 bp were discarded. For the remaining sequences, the transcript coding potential values were predicted by the Coding Potential Calculator (CPC) [44]. Each transcript with a CPC score <0 was considered a long non-coding RNA. mRNA transcripts (CPC scores >0) were also identi ed from the transcriptome in this work. Differentially expressed mRNAs or lncRNAs were identi ed using the R package NOISeq 2.31 (https://bioconductor.org/packages/release/bioc/html/NOISeq.html).

Validation of several LncRNAs using RT-qPCR
Tissues were collected from shoots and roots grown with or without iron for 10 days. Total RNA was isolated using TRIzol (Invitrogen). cDNA was synthesized from total RNA using a cDNA Synthesis Kit (TIANGEN), and RT-qPCR was performed on a LightCycler480 machine (Roche) with SYBR Green Supermix (CWBIO). ACTIN mRNA was used as the internal control for sample normalization. Means ± SD were calculated by three biological repeats. The RT-qPCR primers (synthesized by TSINGKE) are shown in Table S8.

Prediction of the LncRNA-derived miRNAs and target genes
For miRNA precursor analysis, the miRNA sequences and their locations in the genome were acquired from PmiREN (http://www.pmiren.com/) [45]. An miRNA was de ned as a lncRNA-derived miRNA if the pre-miRNA region in the genome was located within a lncRNA. The online software psRNATarget (http://plantgrn.noble.org/psRNATarget/) was used to predict target genes of miRNAs with a maximum expectation of 2.0 [46]. Less than two mismatches and G/U pairs were allowed within the mRNA and miRNA pairing regions.

Prediction and annotation of DE-LncRNA targets
The potential target genes of the differentially expressed lncRNAs (DE-lncRNA) were predicted based on the possibilities of cis and trans interaction nodes between the lncRNAs and mRNAs. For cis-target analysis, we searched the coding genes located within 10 kb upstream or downstream of the DE-lncRNAs. For trans-target analysis, interactions between DE-lncRNA and DE-mRNA were predicted thorugh complementary pairing of bases. The LncTar [47] tool was used for predicting target genes of the lncRNAs. The free energy and standard free energy of paired sites were calculated, and the target genes with standard free energy threshold <-0.1 were considered trans-target genes, while those <-0.2 were considered cis-target genes of lncRNAs. The online software agriGO (http://systemsbiology.cau.edu.cn/agriGOv2) was used to do the GO enrichment, and only those biological process terms with p <0.001 were considered as signi cantly enriched GO terms. The ironde ciency, lncRNA-mRNA response networks were built using Cytoscape [48], which only contained the DE-lncRNAs and the trans-targets which had been reported as important Fe regulators.

Availability of data and materials
All sequencing data generated inthis study has been submitted to NCBI database (Accession number: PRJNA527175).  Tables   Due to technical limitations, table 1 is only available as a download in the Supplemental Files section. Table S1. LncRNAs and mRNA signi cantly up-or down-regualted by Fe de ciency in shoot and root.