Barley long non-coding RNAs (lncRNA) responsive to excess boron

Genomics xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Genomics journal homepage: www.elsevier.com/locate/ygeno

Barley long non-coding RNAs (lncRNA) responsive to excess boron ⁎

Turgay Unvera, , Huseyin Tombuloglub a

Ficus Biotechnology, Ostim Teknopark, No: 1/1/76, 06378, Yenimahalle, Ankara, Turkey Department of Genetics Research, Institute for Research and Medical Consultations (IRMC), Imam Abdulrahman Bin Faisal University, P.O. Box 1982, Dammam 31441, Saudi Arabia

b

A R T I C LE I N FO

A B S T R A C T

Keywords: Boron Endogenous target mimicry Hordeum vulgare Long non-coding RNA Novel transcript discovery and annotation

Long non-coding RNA (lncRNA) has a misleading name, since although they do not encode proteins, they may encode small peptides. Such transcripts are emerging as regulatory molecules. With the advent of next-generation sequencing technologies and novel bioinformatics tools, a tremendous amount of lncRNAs have been identified in several plant species. Recent reports demonstrated roles of plant lncRNAs such as development and environmental response. Here, we reported a genome-wide discovery of ~8000 barley lncRNAs and measured their expression pattern upon excessive boron (B) treatment. According to the tissue-based comparison, leaves have a greater number of B-responsive differentially expressed lncRNAs than the root. Functional annotation of the coding transcripts, which were co-expressed with lncRNAs, revealed that molecular function of the ion transport, establishment of localization, and response to stimulus significantly enriched only in the leaf. On the other hand, 32 barley endogenous target mimics (eTM) as lncRNAs, which potentially decoy the transcriptional suppression activity of 18 miRNAs, were obtained. Also, six lncRNAs, differentially expressed upon B-treatment, were selected and quantitatively analyzed in both B-sensitive and B-tolerant cultivars treated by excess B-level. Quantitative real-time polymerase chain reaction (qRT-PCR) analysis confirmed the B-responsive expressional changes obtained by RNA sequencing. Notably, some lncRNAs (i.e., TCONS_00045190 and TCONS_00056415) over-expressed only in B-tolerant cultivar upon excess B treatment. Presented data including identification, expression measurement, and functional characterization of barley lncRNAs suggest that B-stress response might also be regulated by lncRNA expression, via cooperative interaction of miRNA-eTM-coding target transcript modules.

1. Introduction Long non-coding RNAs (lncRNA) are known as non-protein coding transcripts longer than 200 nt. Although the information about their functions are so limited, studies revealed that they have several direct and indirect roles in the transcriptional, post-transcriptional or posttranslational processes such as gene expression, chromatin modification, transcriptional regulation, and conformational changes in proteins (reviewed by [33]). They act as the precursor of micro-RNA (miRNA) and short interfering-RNA (siRNA). For instance, five lncRNA in Arabidopsis (npc34, npc351, npc375, npc520, and npc523) matched with 24-nt siRNAs from both strands, suggesting these lncRNAs are siRNA precursor. Moreover, it was reported that plant miRNAs, miR869a and miR160c, which were derived from lncRNAs of npc83 and npc521, respectively [4]. Additionally, it was also discovered that miR675 is derived from a mouse lncRNA, H19, and extensively expressed in embryonic liver [14,27].

⁎

Recently, short peptide-coding sequences were discovered in the non-coding regions of plant primary-miRNAs (pri-miRNA) and called as miRNA-encoded peptide (miPEP), which increases the transcription of pri-miRNA [30,57]. Therefore, it was suggested that some plant lncRNAs might have peptide coding potential [33]. In addition, miRNA activity can be regulated by endogenous target mimicry (eTM) molecules, being a type of lncRNA [26]. Such an example, an endogenous lncRNA called Induced by Phosphate Starvation 1 (IPS1) of Arabidopsis thaliana binds to miR399 to inhibit the cleavage of the miR399 target transcript [18]. Other than these functions, plant lncRNAs were found to be involved in many regulatory mechanisms, such as histone modeling [21], promoter modification [15,62], protein re-localization [46], and alternative splicing [3]. To date, lncRNAs have been identified extensively in mammals, in which the human genome includes > 56,000, and mice have almost 46,000 lncRNAs [58]. Publicly available databases such as LncRNAdb (http://lncrnadb.com), a database for functional lncRNAs, harbor

Corresponding author. E-mail address: turgay.unver@ficusbio.com (T. Unver).

https://doi.org/10.1016/j.ygeno.2019.11.007 Received 28 September 2019; Received in revised form 6 November 2019; Accepted 11 November 2019 0888-7543/ © 2019 Published by Elsevier Inc.

Please cite this article as: Turgay Unver and Huseyin Tombuloglu, Genomics, https://doi.org/10.1016/j.ygeno.2019.11.007

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

transcriptome libraries each including pooled RNAs from tree biological replicate previously generated by our group (SUB337351 and SUB2170217) [50]. Briefly, a total of 208,249,690 clean sequencing reads from four paired-end libraries (50_leaf; 52, 422,032, 50_root; 52,168,358, 1000_leaf; 52,305,062, and 1000_root, 51,354,238) were utilized in this study. We first removed the adapter sequences and lowquality reads from the sequencing reads with Trimmomatic v0.36 [5], and these clean reads were aligned to the barley reference genome (ASM32608v1 assembly) by TopHat2 v2.1.1 with default parameters [29]. Afterward, genome-aligned reads were assembled ab initio using popular transcriptome analysis suit Cufflinks v2.2.1 [55] to build potential transcript structures. All gene transfer format (GTF) files produced in the assembly step were merged with Cuffmerge tool and transcript features were queried against the Ensembl Plant database (release 33) [28] by Cuffcompare to discover previously unannotated transcript sequences. Among all un-annotated transcripts, including coding and non-coding sequences, we obtained lncRNAs as follows: (i) we first tested the coding potential of each transcript individually using TransDecoder (https://transdecoder.github.io/) and filtered out those with an open reading frame having > 100 amino acids, (ii) we then removed transcripts which were shorter than 200 nucleotides in length, (iii) to remove housekeeping lncRNA species, we queried all potential lncRNAs against to non-coding RNA family database, Rfam (v12.1) [37], with Infernal tool (v1.1.1; cutoff E-value ≤0.001) [36]. Then, potential miRNA precursors, tRNAs, etc. were removed and were not included for further analysis, and (iv) we aligned all potential lncRNA transcripts to the Swiss-Prot database (release 2017_01) [49] using Blastx (v2.5.0; cutoff E-value ≤0.001) [8] to eliminate transcripts with potential protein-coding ability.

functionally annotated lncRNAs, of the majority, belong to the human. A few databases were so far released for plant lncRNAs; such as the Green Non-Coding Database (GreeNC, http://greenc.sciencedesigners. com) [39] and CANTATAdb (http://cantata.amu.edu.pl) [47], which provides information for around 45,117 lncRNAs from several plant species. RNA-sequencing or deep transcriptome analysis is an important technology, which provides information not only for protein-coding transcripts but also for non-coding RNAs (such as miRNA, siRNA, piwiinteracting RNA (piRNA), and small nucleolar RNA (snoRNA) as well as lncRNAs). In addition, it allows distinguishing lncRNAs expressed in different tissues or cells. Boron (B) is an essential micronutrient for plants and its unfavorable concentration negatively affects plant growth and productivity where the soils having insufficient or excess B [56]. The range of sufficient B concentration in soil is so limited, thus generally plants suffer from either B-deficiency or -toxicity problem that is common in agricultural soils around the world [9]. However, plants growing in B-contaminated soils must tolerate the excess level of B to survive. In the last decade, many studies conducted to understand the cellular mechanisms underlying to balance cellular B content in plants [35,50]. Facilitated transport of boric acid (a regular form of B in soil) by transporter channels were suggested to be the molecular regulatory mechanism. Several B-importer and exporter proteins have recently been identified as B-transporters to regulate its cellular homeostasis [35]. On the other hand, novel sequencing-based approaches to discover the transcriptional response at genome-wide level are being extensively utilized in plants faced with unfavorable environmental conditions. In this context, to quantify gene expression and to annotate coding-transcripts, we performed a high-throughput genome-wide transcriptome analysis of barley tissues treated with excess B (1 mM), previously [50]. In addition, we also screened B-responsive miRNA expression pattern [38]; identified MYB type transcription factors (TF) [51] and water channel Aquaporins (AQP) [52] to understand B homeostasis of barley. These studies helped to observe the main or possible players involved in B regulation. Besides the emerging evidence suggested that the molecular regulatory mechanism is so complicated and not only limited to the activity of those coding transcripts. LncRNAs have been discovered to be novel players, but also as regulatory molecules in the regulation of gene expression. To date, a large set of RNA-seq libraries were used to identify lncRNAs in genome-wide scale or tissue/condition/inoculation-specific manner [11,31,32]. Stress-responsive lncRNAs were examined from the RNA-seq data of the plants under the distinct type of stress conditions as well (reviewed by [10,44,61]). Qi et al. [40] identified 584 lncRNAs, which were responsive to drought stress in foxtail millet. 125 putative lncRNAs were identified in wheat, responsive to powdery mildew infection and heat stress [59]. Detailed examination of a large set of poplar (Populus trichocarpa) RNA-seq data revealed 504 lncRNAs in response to drought [45]. Additionally, Huang et al. [23] reported over 12,000 barley lncRNAs, of them 604 were Fusarium head blight inoculation responsive [23]. However, no such a study to profile the expression level of lncRNAs under the B-excess as one of the abiotic stress conditions was conducted until now. In this study, we identified and quantitatively compared the expression of B-responsive barley lncRNAs from four transcriptome datasets. Tissue-specific (root and leaf) and excess B-responsive lncRNAs, which were co-expressed with coding-transcripts, discovered and comprehensively analyzed. Besides, the expressional status of some selected lncRNAs were analyzed by qRTPCR both in B-sensitive and B-tolerant cultivars of barley.

2.2. Expression pattern analysis of coding and non-coding RNAs Transcript abundances in each library were measured with Kallisto v0.43.0 [7]. Then the transcripts expressed < 1 Transcripts Per Million (TPM) in all libraries were considered as transcriptional noise and were removed from further downstream analysis steps. Differentially expressed transcripts within each group (leaf and root samples) were determined by calculating fold changes of TPM values in RNA-seq datasets. Transcripts with differential expression values ≥2 fold-changes in compared datasets were classified as boron responsive. The Gene Ontology (GO) analysis of differentially expressed genes was performed using an online GO analysis toolkit, agriGO applying default parameters [16].

2.3. Co-expression analysis of lncRNAs with coding mRNAs and prediction of endogenous target mimicry (eTM) sequences We predicted the putative functions of differentially expressed lncRNAs with “guilt-by-association” approach, which employed in earlier studies for lncRNA annotation [43] [20] [12]. To show potential lncRNA-mRNA associations, we identified co-expressed mRNA-lncRNA pairs with Spearman's correlation test in R v3.1.0 statistical computation environment [48]. Then, co-localized mRNA-lncRNA pairs on the reference genome were identified with Bedtools v2.25.0 [41]. We considered the mRNA-lncRNA pair as co-expressed if the Spearmen's rho is equal or > 0.90 (p-val < 0.01) between the expression values of coding and lncRNA transcripts, and as co-localized when the distance between two transcripts were < 10 kb. To dissect putative eTM sequences among the transcripts annotated as lncRNA, we employed our analysis pipeline previously introduced by our group [26]. In the eTM sequences analysis pipeline, we utilized mature miRNA sequences of barley collected from miRBase (release 21) [19]. To identify potential target transcripts of barley miRNAs, we used psRNATarget, an online miRNA target analysis tool [13] as previously described [1,2,17,24,60].

2. Materials and methods 2.1. Identification of barley lncRNAs To study global expression profiling of barley (Hordeum vulgare cv Sahara) lncRNAs under boron stress condition, we utilized four 2

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

B

15 log2 TPM with 1 offset

900

600

300

10

5

f

0_ r 00 X1

X5

Chromosomes

lea

7

0_

6

00

5

X1

4

ro ot

f

3

lea

2

0_

1

oo t

0

0

X5 0_

Total number of lncRNAs

A

50_root

1000_root

1000_leaf

Color Key

50_leaf

C

0 5 10 15 log2 TPM values

Fig. 1. Chromosomal distribution of barley lncRNAs and their expression in normal and stress conditions. A. Total number of lncRNAs localized to each barley chromosome (from chromosome 1 to 7). B. Box plot representation of expression levels of lncRNAs before and after boron exposure in root and leaves. C. Expressional changes of lncRNAs across samples are being illustrated as heat map graph.

synthesize the first strand of cDNA. The procedure was adapted according to the suggested protocol of the supplier (Maxima First Strand cDNA Synthesis Kit for RT-qPCR, Thermo Scientific). The reaction mixture including 5xReaction Mixture, Maxima Enzyme Mix, deionized water, and RNA sample was kept in a thermal block for 10 min at 25 °C followed by 15 min at 50 °C. Afterward, the reaction was stopped by incubating the tubes for 5 min at 85 °C. The samples were stored in a refrigerator (−20 °C) [53]. Transcriptome datasets revealed hundreds of B-responsive differentially expressed lncRNAs. Among these genes, six lncRNAs were randomly selected for the validation of their expression by qRT-PCR. A real-time qPCR instrument (Applied Biosystems, 7500 Fast Real-Time PCR System) was used to determine the transcript abundances. Genespecific primers were designed by using Primer Express 3.0.1 (Life Technologies) software (Table S6). Gene amplification was carried out using Power SYBR Green PCR Master Mix kit (Thermo Scientific). The qPCR reaction was initiated by incubating the reaction mixture at 95 °C for 10 min. Then, 40 cycles of 95 °C (15 s) following with 65 °C (45 s) was applied. Finally, a melting curve analysis was conducted by heating the tubes from 60 °C to 95 °C. 2-ΔΔCt method was used to compute the relative expression level of the transcripts [42]. The expression of genes was normalized by using barley 18S rRNA [6,54].

2.4. Plant growth and stress treatment Barley seeds (Hordeum vulgare cv. Sahara (B-tolerant) and cv. Clipper (B-sensitive)) were kindly provided from Prof Julie Hayes and Prof Tim Sutton from Adelaide University, Australia. The seeds were surface-sterilized by immersing in a 4% H2O2 solution (3 min) followed by washing them treble with deionized (DI) water and placed in between two layers of sterile tissue papers in Petri plates. After dark incubation for four days, the seedlings (n = 5) were transferred into hydroponic growing system filled with modified Hoagland solution [22] including either 0.05 mM (control) or 1 mM (B-excess) of boric acid. The solution was aerated for a regime of 3 h per day by an aquarium pump and renewed for three days. The seedlings were grown for three weeks in greenhouse conditions: relative humidity of 60–70%, temperature interval of 23–26 °C, and 16/8 h light/dark cycle. Afterward, the leaf and root tissues were harvested, immediately frozen by liquid nitrogen, and stored at −80 °C until RNA extraction. 2.5. qRT-PCR validation assay Total RNA was extracted using Trizol® reagent (Invitrogen) according to the recommendations of the supplier. Quality assessment of the RNA (integrity and quantity) was carried out by using NanoDrop 2000 (Thermo Scientific). Five μg of RNA samples were treated by five units of RNase-free dsDNase enzyme (Thermo Sci.) in order to remove contaminating genomic DNA. One μg of RNA sample was used to 3

Tr ph an os sfe ph ra or se Pr us ac ot −c tiv ei on ity n se ta , tr in an rin in s e/ g fe th gr rr re ou ing on ps in e ki na se ac tiv Pr Ph ity ot ei os n ph k in ot as ra e ns ac fe tiv ra se ity ac Pr tiv ot i t ei Li y, n ga al bi se co nd ho ac in lg g tiv ro ity up ,f or as m in ac g ce ca pt rb or on −n itr og en bo nd s Li ga se ac tiv ity K C i n ys as te e in ac e− tiv ty ity pe pe pt id as e ac tiv ity

Counts

Root

517

FDR

Leaf

5092

Counts

Root

A

FDR

Leaf

R e R es spo esp po ns o ns e nse Tra Pu n e to rin to ot to spo e Re o h st r Pu ribo spo R rga er o imu t rin nu ns es nic rg lus e cle e po s an nu o to n ub is cl tid ex se st m eo e te to an tid me rn fu ce e ta al ng Pr m bo st u et li im s ot Pr ab c u e ot i ol pro lus ei P n m ic c n ro o am te di P pr ess Po in fic r oc st i n −t m a o e ra Pr o a et tion teo ss ns im cid ab p lys la ar p ol ro is tio y ho ic ce na m sp pr ss Ph l p eta ho oc os ro bo ry es te lic la s p Ph ho in t os ru P m pro ion ph s m ho od ce at e sp ific ss e ta h a m bo or tio y e M M O tab lic p lati n xi o r on ac a da lic oc ro cr m om M tio pr es ol o et n oc s ec le ab re e ul cu ol du ss e le ic c m m p tio et o ro n ab d ce C Es C el el ol ific ss lu tab lu ic at la lis la C L pr ion rp h rm e ro me Io oca oce ac llu t l ei nt n liza ss ro ar n C mo ma m of l tran tion el le c C et oc sp lu cu ro el a a o C lar lar mo lula C boli liza rt el m c le r e c lu ac o c m llu p tion la r m ul e la ro r c om pl e ta r ce ar o ex me bo pro ss bo lec su ta lic c hy ul bu bo pr ess dr ar ni lic oc at co t o p es e r bi mp rga oce s os le ni s yn x a za s th ss tio C etic em n at p b io ro ly n tra ces ns s po rt R

T. Unver and H. Tombuloglu

Genomics xxx (xxxx) xxx–xxx

Leaf Root

B

2340

3. Results and discussion

3.1. Barley lncRNA identification

After the adapter sequence trimming and removal of low-quality

4

Leaf Root

1799 531 1242

C

25 50 75 0.04

0.03

0.02

0.01

D 10 15 20 0.04

0.03

0.02

0.01

Gene Ontology Terms

Fig. 2. Differentially expressed transcripts and GO term enrichment analysis. A-B. Venn diagram of differentially expressed coding (A) and lncRNA (B) transcripts in leaf and root. C-D. Significant biological process (C) and molecular function (D) GO terms associated with differentially expressed coding transcripts.

reads, the mean library size of four sequencing libraries included in the study was over 20 million (min: 18,223,418, max: 22,171,430). Additionally, we observed an average of 82.35% (min: 78.7%, max: 90.9%) overall read mapping for the alignment step and considered all sequencing libraries had sufficient quality to perform an ab initio

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

A

B TCONS_00056131 TCONS_00002116 MLOC_65770.2

TCONS_00056131

TCONS_00062668

MLOC_62874.1

TCONS_00002116

TCONS_00061958 TCONS_00080191

MLOC_62874.1 MLOC_32738.1

TCONS_00061521 MLOC_6187.1

MLOC_6187.1

TCONS_00062668

MLOC_62200.3

TCONS_00080191 MLOC_65770.2

TCONS_00061958 MLOC_25879.3

MLOC_25879.3 MLOC_62200.3

TCONS_00061521

1000_leaf

50_leaf

1000_root

50_root

MLOC_32738.1

Color key

0

C

D

2 4 6 8 log2 TPM values

TCONS_00049524

TCONS_00056185

MLOC_58327.1 TCONS_00021410

TCONS_00026166

TCONS_00000654 MLOC_63022.1

MLOC_81784.1

TCONS_00049499

MLOC_10584.1

MLOC_10584.1

MLOC_36532.2

MLOC_62887.1 MLOC_62874.1

TCONS_00023381

TCONS_00015537 TCONS_00026166 TCONS_00043598

TCONS_00000654 TCONS_00056185

MLOC_79529.1 MLOC_69539.4

MLOC_38182.2

TCONS_00049524

MLOC_62887.1

MLOC_81784.1

MLOC_23387.1 MLOC_38182.2

TCONS_00061958

MLOC_23387.1

TCONS_00049523 TCONS_00054322 MLOC_58327.1

TCONS_00055156

MLOC_36532.2

TCONS_00049523 TCONS_00005023 TCONS_00055156

MLOC_79529.1

TCONS_00023381

MLOC_62874.1 TCONS_00021410 TCONS_00049499 TCONS_00054322 TCONS_00061958 TCONS_00043598 MLOC_69539.4 TCONS_00015537 MLOC_63022.1

1000_leaf

50_leaf

1000_root

50_root

TCONS_00005023

Fig. 3. Ion transport and response to stimulus associated lncRNAs. A-C. LncRNAs and coding transcripts related to ion transport (A) and response to stimulus (C) GO terms. LncRNAs and coding transcripts are represented as circles and round rectangles, respectively. Differentially expressed lncRNAs before and after boron exposure are shown as colored. Green colored lncRNAs are differentially expressed only in leaves. However, lcnRNAs are colored as red if they are differentially regulated in both root and leaves. LncRNAs and associated coding transcripts (co-expressed and co-localized) are linked to each other with network edges. B-D. Heat map representation of expression levels both lncRNAs and coding transcripts associated with ion transport (B) and response to stimulus (D) GO terms. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

levels of most lncRNAs compared to protein-coding genes make it more difficult to detect lncRNAs [34]. Generally, they are excluded from the total lncRNA pool resulting in fluctuations of total lncRNA number. But it is important to note that lncRNAs with low expression may have a big impact, thus extensive and a deep pipeline is required to extract lncRNA, which may have important biological functions. In general, the distribution of lncRNAs to the barley chromosomes is proportional with its chromosome sizes, except chromosome 2, which includes the highest number (Fig. 1A).

transcriptome reconstruction analysis. The run of Cufflinks pipeline was led to the identification of as many as 34,000 previously unannotated intergenic transcripts, of which 10,439 were the lack of coding potential and > 200 bp in length (Table S1). When we filtered out the transcripts expressed at low levels (< 1 TPM in all samples), we obtained 8009 intergenic putative lncRNAs in the final list, which were distributed almost equally to all barley chromosomes (Fig. 1A). However, chromosome 2 was observed to be the richest one in terms of the total number of lncRNAs it harbors. In this study, a total number of ~8000 barley lncRNAs were identified which is smaller than that of human (~56,000) [58] and mouse (~46,000); higher than fruit fly (~3300) [11], and poplar (2542) [45] (Table S1). Actual numbers of lncRNAs can be altered depending on the sample examined. In this analysis, four RNA-Seq libraries were used to detect total lncRNAs. More lncRNAs can be found from the barley genome by increasing the number of RNA-Seq sets from distinct tissues and/or conditions. In general, low expression

3.2. Expression pattern of barley lncRNAs and coding transcripts upon excess B-treatment As the expression profiles of lncRNAs in root and leaf samples were examined, it was determined that expression levels of lncRNAs in the samples collected from the same tissue were similar to one another 5

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

A TCONS_00038547

MLOC_15016.7 MLOC_65770.2

MLOC_23387.1 MLOC_81784.1

TCONS_00042959 MLOC_62200.3

MLOC_25879.3

TCONS_00002116 TCONS_00056185

TCONS_00061958 TCONS_00062668

TCONS_00003436 MLOC_62874.1 MLOC_6187.1

TCONS_00049499

TCONS_00061521 TCONS_00080191

MLOC_17577.1 TCONS_00056131

MLOC_63022.1 TCONS_00026166 MLOC_32738.1

B TCONS_00056131 TCONS_00002116 TCONS_00062668 TCONS_00061958 TCONS_00038547 TCONS_00042959 TCONS_00061521 MLOC_81784.1 MLOC_17577.1

MLOC_15016.7 TCONS_00026166 TCONS_00003436

0

MLOC_63022.1

Color key

TCONS_00080191

2 4 6 8 log2 TPM values

MLOC_6187.1

MLOC_23387.1 MLOC_62200.3 TCONS_00056185 MLOC_32738.1 MLOC_65770.2 MLOC_25879.3 TCONS_00049499

1000_leaf

50_leaf

1000_root

50_root

MLOC_62874.1

Fig. 4. Establishment of localization associated lncRNAs. A. Network representation of establishment of localization related coding transcripts and lncRNAs (the same color scheme as in Fig. 3 is used). B. Expression levels of establishment of localization GO term associated coding and lncRNAs.

3.3. Functional annotation of lncRNAs co-expressed with coding transcripts

(Fig. 1B). The hierarchical clustering analysis, however, revealed that particular lncRNA clusters were expressed at relatively higher levels specific to tissue types (Fig. 1C). Differential expression analysis of both coding and lncRNA transcripts showed that there were 2-fold or more change (up- or down-regulation) in the expression of the vast amount of transcripts in response to boron stress in leaves and root tissues (Table S2). We observed that the number of common coding transcripts that were differentially regulated in both tissues was 517; in addition, the total number of differentially expressed coding transcripts in the leaf tissue was roughly doubled as compared to root tissue (Fig. 2A). Similar to the differential expression analysis of coding transcripts, we detected a greater number of up- or downregulated lncRNAs in the leaves than the root samples (Fig. 2B).

The GO enrichment analysis of differentially expressed coding transcripts (Table S3) revealed that ion transport (GO:0006811), establishment of localization (GO:0051234), and response to stimulus (GO:0050896) terms significantly enriched (FDR < 0.05) only in the leaf samples (Fig. 2C). In addition to this, molecular function terms of ligase activity (GO:0016874) and cysteine-type peptidase activity (GO:0008234) were significant and specific to the leaves (Fig. 2D). Based on the co-expression and co-localization analysis of barley lncRNAs and coding transcripts, we observed potential lncRNAs in association with the ion transports, establishment of localization, and response to stimulus-related transcripts (Figs. 3 and 4). We determined six lncRNAs, which strongly linked to ion transport-related genes (Fig. 3A). However, only one lncRNA (TCONS_00061958) was showing ≥2 fold expressional change in both the leaf and root samples (Fig. 3B). 6

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

both leaf and root tissues. Studies also revealed the differential expression pattern of plant lncRNAs such as Arabidopsis, wheat, and poplar in response to biotic and abiotic stresses [32,45,59]. Here, we identified that barley boron-responsive lncRNAs are expressed in a tissue-specific manner (Figs. 2–4 and Table S2). According to GO term enrichment analysis, differentially expressed coding transcripts were categorized into three biological processes: response to stimulus (50 transcripts), ion transport (17 transcripts), and establishment of localization (52 transcripts), which were leaf specific in this analysis (Fig. 2C and Table S3). These findings were consistent with our earlier report [50], where the biological process, ion transport, and establishment of localization categories were found to be enriched in both leaf and root tissues.

Table 1 Computationally identified putative eTM sequences having potential to act as miRNA sponge. miRNA name

Corresponding eTM sequence

hvu-miR1130 hvu-miR159a hvu-miR159b hvu-miR166a hvu-miR166b hvu-miR166c hvu-miR399 hvu-miR5049a hvu-miR5049f hvu-miR6180 hvu-miR6186 hvu-miR6189 hvu-miR6191 hvu-miR6196

TCONS_00051546 TCONS_00004839 TCONS_00004839 TCONS_00012368 TCONS_00012368 TCONS_00012368 TCONS_00043651 TCONS_00032652 TCONS_00022657, TCONS_00073195, TCONS_00034493 TCONS_00020181 TCONS_00022441, TCONS_00002822, TCONS_00076928 TCONS_00054070 TCONS_00000261, TCONS_00007423 TCONS_00001077, TCONS_00014009, TCONS_00063637, TCONS_00072190

hvu-miR6197 hvu-miR6202 hvu-miR6212 hvu-miR6214

TCONS_00022661, TCONS_00022662 TCONS_00074388

TCONS_00045158 TCONS_00011514, TCONS_00058908,

3.4. eTM sequence discovery Here, we predicted 32 barley eTMs, which might decoy the transcriptional suppression activity of 18 miRNAs, including conserved barley miRNAs such as miR159a, miR166a, and miR399 (Table 1 and Table S4). In the GO enrichment analysis of miRNA target genes using the agriGO tool, we found 102 significant (FDR < 0.05) different GO terms in three domains, consisting of molecular function, biological process, and cellular component (Table S5). Among the terms, the most significant one was “protein amino acid phosphorylation” term (GO:0006468, FDR = 4.70E-58). Our findings in eTM analysis suggest that barley lncRNAs might regulate several distinct cellular and molecular processes via mimicking specific miRNA target transcripts. As it was previously reported by recent studies, lncRNAs might behave as a mimicry-transcript that targeted by miRNAs and fate it to degradation [25]. It was first reported in Arabidopsis, over-expression of non-coding gene IPS1 suppressed the miR399 expression that resulted in higher expression of the miR399 target [18]. On the other hand, we determined and measured the boron-responsive barley miRNAs [38].

TCONS_00078999 TCONS_00009121, TCONS_00012725, TCONS_00044421, TCONS_00058276, TCONS_00059667, TCONS_00061965,

A great majority of response to stimulus-related lncRNAs increased their expression (≥ 2 fold) in leaf samples upon the B exposure (Fig. 3C–D). Additionally, we detected 12 lncRNA sequences, which had similar expression patterns and co-localized with the “establishment of localization” related genes (Fig. 4). We also detected that some of the lncRNAs were only differentially expressed in tissue-specific manner such as TCONS_00002116 expression in leaf under B-excess (Fig. 4B). However, TCONS_00061958 was found to be differentially regulated in

Fig. 5. Relative expression level of selected lncRNAs in B-tolerant (Sahara) and B-sensitive (Clipper) cultivars of barley. Barley seedlings were subjected by 0.05 mM (control) and 1 mM of boric acid. The gene abundances were analyzed in both root and leaf tissues by qRT-PCR. The data was normalized by using 18S rRNA genes. 7

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

to the tolerant variety, thus effectively contributing to boron-excess tolerance in them. This has obvious implications for barley breeding, to improve tolerance to boron-excess. Supplementary data to this article can be found online at https:// doi.org/10.1016/j.ygeno.2019.11.007.

Accordingly, miR5049 was down-regulated in B-stressed leaf (three times than control leaf). Also, miR399 was over-expressed in leaf and suppressed in root tissue upon B-exposure (three times of each tissue than that of control). In this study, miR5049 and miR399 were also found to be as regulated miRNAs upon B-treatment by eTM analysis (Table S4) where TCONS_00032652 and TCONS_00043651 have predicted as the target mimic sequences, which are able to decoy the miRNA activities, respectively. Thus, the expression of transcripts targeted by these miRNAs might be altered due to differential expressions of lncRNAs. Phosphate transporter 2 (PHO2) and putative ubiquitin-conjugating enzyme (UBC) were found to be the target genes of miR399. Also, tubby protein-like transcript was determined as the miR5049 target. In the transcriptome analysis, ubiquitin carboxyl-terminal hydrolase gene was highly up-regulated in leaf tissue upon excess B treatment [50]. Interestingly, expression profiles of miRNA and its corresponding lncRNA target mimic transcript provide insights into the regulation of B-stress in plants. For instance, lncRNA TCONS_00043651, a potential target mimic sequence of miR399, upregulated in roots (three times than that of control) upon B-exposure. Oppositely, miR399 expression was reflected with the same pattern in a negative direction (three times down-regulated). Similarly, five times increase of lncRNA TCONS_00043651 in leaf tissues may prevent the expression of miR399, which was up-regulated only three times (Table S2). These preliminary results suggest that boron regulation can be cooperatively controlled by the interaction of miRNA-eTM (lncRNA)coding target transcript modules.

Author contributions TU organized the study. TU and HT performed analyses, and interpreted the data. TU and HT wrote the manuscript. Declaration of Competing Interest The authors declare no conflict of interest. Acknowledgements The authors greatly appreciate Prof Julie Hayes and Prof Tim Sutton for providing B-tolerant and B-sensitive cultivars. References [1] G. Akdogan, E.D. Tufekci, S. Uranbey, T. Unver, miRNA-based drought regulation in wheat, Funct Integr Genomics 16 (2016) 221–233, https://doi.org/10.1007/ s10142-015-0452-1. [2] Y. Bakir, V. Eldem, G. Zararsiz, T. Unver, Global transcriptome analysis reveals differences in gene expression patterns between nonhyperhydric and hyperhydric peach leaves, Plant Genome 9 (2016), https://doi.org/10.3835/plantgenome2015. 09.0080. [3] F. Bardou, F. Ariel, C.G. Simpson, N. Romero-Barrios, P. Laporte, S. Balzergue, J.W. Brown, M. Crespi, Long noncoding RNA modulates alternative splicing regulators in Arabidopsis, Dev. Cell 30 (2014) 166–176, https://doi.org/10.1016/j. devcel.2014.06.017. [4] B. Ben Amor, S. Wirth, F. Merchan, P. Laporte, Y. d’Aubenton-Carafa, J. Hirsch, A. Maizel, A. Mallory, A. Lucas, J.M. Deragon, H. Vaucheret, C. Thermes, M. Crespi, Novel long non-protein coding RNAs involved in Arabidopsis differentiation and stress responses, Genome Res 19 (2009) 57–69, https://doi.org/10.1101/gr. 080275.108. [5] A.M. Bolger, M. Lohse, B. Usadel, Trimmomatic: a flexible trimmer for Illumina sequence data, Bioinformatics 30 (2014) 2114–2120, https://doi.org/10.1093/ bioinformatics/btu170. [6] S.M. Bostancioglu, G. Tombuloglu, H. Tombuloglu, Genome-wide identification of barley MCs (metacaspases) and their possible roles in boron-induced programmed cell death, Mol. Biol. Rep. 45 (3) (2018) 211–225. [7] N.L. Bray, H. Pimentel, P. Melsted, L. Pachter, Near-optimal probabilistic RNA-seq quantification, Nat. Biotechnol. 34 (2016) 525–527, https://doi.org/10.1038/nbt. 3519. [8] C. Camacho, G. Coulouris, V. Avagyan, N. Ma, J. Papadopoulos, K. Bealer, T.L. Madden, BLAST+: architecture and applications, BMC Bioinform 10 (2009) 421, https://doi.org/10.1186/1471-2105-10-421. [9] J.J. Camacho-Cristobal, E.M. Martin-Rejano, M.B. Herrera-Rodriguez, M.T. Navarro-Gochicoa, J. Rexach, A. Gonzalez-Fontes, Boron deficiency inhibits root cell elongation via an ethylene/auxin/ROS-dependent pathway in Arabidopsis seedlings, J. Exp. Bot. 66 (2015) 3831–3840, https://doi.org/10.1093/jxb/erv186. [10] J.A. Chekanova, Long non-coding RNAs and their functions in plants, Curr. Opin. Plant Biol. 27 (2015) 207–216, https://doi.org/10.1016/j.pbi.2015.08.003. [11] M.J. Chen, L.K. Chen, Y.S. Lai, Y.Y. Lin, D.C. Wu, Y.A. Tung, K.Y. Liu, H.T. Shih, Y.J. Chen, Y.L. Lin, L.T. Ma, J.L. Huang, P.C. Wu, M.Y. Hong, F.H. Chu, J.T. Wu, W.H. Li, C.Y. Chen, Integrating RNA-seq and ChIP-seq data to characterize long non-coding RNAs in Drosophila melanogaster, BMC Genomics 17 (2016) 220, https://doi.org/10.1186/s12864-016-2457-0. [12] E. D'Haene, E.Z. Jacobs, P.J. Volders, T. De Meyer, B. Menten, S. Vergult, Identification of long non-coding RNAs involved in neuronal development and intellectual disability, Sci. Rep. 6 (2016) 28396, , https://doi.org/10.1038/ srep28396. [13] X. Dai, P.X. Zhao, psRNATarget: a plant small RNA target analysis server, Nucleic Acids Res. 39 (2011) W155–W159, https://doi.org/10.1093/nar/gkr319. [14] B.K. Dey, K. Pfeifer, A. Dutta, The H19 long noncoding RNA gives rise to microRNAs miR-675-3p and miR-675-5p to promote skeletal muscle differentiation and regeneration, Genes Dev. 28 (2014) 491–501, https://doi.org/10.1101/gad.234419. 113. [15] J. Ding, J. Shen, H. Mao, W. Xie, X. Li, Q. Zhang, RNA-directed DNA methylation is involved in regulating photoperiod-sensitive male sterility in rice, Mol. Plant 5 (2012) 1210–1216, https://doi.org/10.1093/mp/sss095. [16] Z. Du, X. Zhou, Y. Ling, Z.H. Zhang, Z. Su, agriGO: a GO analysis toolkit for the agricultural community, Nucleic Acids Res. 38 (2010) W64–W70, https://doi.org/ 10.1093/nar/gkq310. [17] V. Eldem, U. Celikkol Akcay, E. Ozhuner, Y. Bakir, S. Uranbey, T. Unver, Genomewide identification of miRNAs responsive to drought in peach (Prunus persica) by

3.5. Expression of selected lncRNAs in B-tolerant and B-sensitive barley To validate lncRNA expressions obtained by transcriptome analysis, six differentially expressed genes were selected and analyzed in the root and leaf tissues of barley by qRT-PCR (Fig. 5). In addition to B-resistant cultivar (Sahara), B-susceptible (Clipper) barley was treated by excess-B and the genes were quantitatively analyzed. Consistent with the data from RNA-seq, all selected genes in cv. Sahara was found to be B-responsive at least in one tissue (i.e., root or leaf). Compared to the untreated control, the abundance of TCONS_00045190 and TCONS_00002619 transcripts dramatically increased in the root of Sahara exposed by excess-B. However, the relative expression of these genes in leaves was limited. Similarly, the expression of these genes in cv. Clipper was found to be relatively less than those obtained from cv. Sahara. Notably, the expression level of TCONS_00056415 reflected a reverse pattern in between cv. Sahara and cv. Clipper, which was respectively down-regulated and up-regulated against excess-B treatment. Obtained results suggest that lncRNAs possess expressional differences in the varieties of barley, which may contribute to acquiring B-toxicity tolerance. 4. Conclusion With the development of next-generation sequencing technologies and advancement in bioinformatics, more transcriptional datasets were generated, including the units with little or no protein-coding potential. In recent years, the lncRNA has been considered as regulatory molecules in several bioprocesses. Though a large number of lncRNA transcripts have been identified in plants using genome-wide approaches, they have not been previously applied to an important crop like barley. Also, the possible involvement of lncRNA in boron-response mechanism has not been previously reported. Here, we report the genome-wide discovery of ~8000 barley lncRNA after excessive boron treatment, further measuring their expression pattern. Furthermore, we functionally annotated the coding transcripts, which are co-expressed with lncRNAs. Curiously, we have shown that cooperative interaction of miRNA-eTM (lncRNA)-coding target transcript modules might regulate the B-response in barley. Interestingly, expression analyses of lncRNA in sensitive and tolerant varieties revealed that some lncRNA are exclusive 8

Genomics xxx (xxxx) xxx–xxx

T. Unver and H. Tombuloglu

[18]

[19] [20]

[21]

[22] [23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35] [36]

[37]

[38]

[39]

nar/gkv1215. [40] X. Qi, S. Xie, Y. Liu, F. Yi, J. Yu, Genome-wide annotation of genes and noncoding RNAs of foxtail millet in response to simulated drought stress by deep sequencing, Plant Mol. Biol. 83 (2013) 459–473, https://doi.org/10.1007/s11103-013-0104-6. [41] A.R. Quinlan, BEDTools: The Swiss-Army tool for genome feature analysis, Curr Protoc Bioinformatics 47 (2014), https://doi.org/10.1002/0471250953.bi1112s47 11 12 1–34. [42] M.W. Pfaffl, A new mathematical model for relative quantification in real-time RT–PCR, Nucleic Acids Res. 29 (9) (2001) e45. [43] J.L. Rinn, H.Y. Chang, Genome regulation by long noncoding RNAs, Annu. Rev. Biochem. 81 (2012) 145–166, https://doi.org/10.1146/annurev-biochem-051410092902. [44] S. Shafiq, J. Li, Q. Sun, Functions of plants long non-coding RNAs, Biochim. Biophys. Acta 1859 (2016) 155–162, https://doi.org/10.1016/j.bbagrm.2015.06. 009. [45] P. Shuai, D. Liang, S. Tang, Z. Zhang, C.Y. Ye, Y. Su, X. Xia, W. Yin, Genome-wide identification and functional prediction of novel and drought-responsive lincRNAs in Populus trichocarpa, J. Exp. Bot. 65 (2014) 4975–4983, https://doi.org/10.1093/ jxb/eru256. [46] C. Sousa, C. Johansson, C. Charon, H. Manyani, C. Sautter, A. Kondorosi, M. Crespi, Translational and structural requirements of the early nodulin gene enod40, a shortopen reading frame-containing RNA, for elicitation of a cell-specific growth response in the alfalfa root cortex, Mol. Cell. Biol. 21 (2001) 354–366, https://doi. org/10.1128/MCB.21.1.354-366.2001. [47] M.W. Szczesniak, W. Rosikiewicz, I. Makalowska, CANTATAdb: a collection of plant long non-coding RNAs, Plant Cell Physiol 57 (2016) e8, https://doi.org/10.1093/ pcp/pcv201. [48] R.C. Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2016 2015. URL http.www.R-project.org. [49] C. The UniProt, UniProt: the universal protein knowledgebase, Nucleic Acids Res. 45 (2017) D158–D169, https://doi.org/10.1093/nar/gkw1099. [50] G. Tombuloglu, H. Tombuloglu, M.S. Sakcali, T. Unver, High-throughput transcriptome analysis of barley (Hordeum vulgare) exposed to excessive boron, Gene 557 (2015) 71–81, https://doi.org/10.1016/j.gene.2014.12.012. [51] H. Tombuloglu, G. Kekec, M.S. Sakcali, T. Unver, Transcriptome-wide identification of R2R3-MYB transcription factors in barley with their boron responsive expression analysis, Mol. Gen. Genomics. 288 (2013) 141–155, https://doi.org/10.1007/ s00438-013-0740-1. [52] H. Tombuloglu, I. Ozcan, G. Tombuloglu, S. Sakcali, T. Unver, Aquaporins in borontolerant barley: identification, characterization, and expression analysis, Plant Mol. Biol. Report. 34 (2016) 374–386. [53] H. Tombuloglu, Genome-wide identification and expression analysis of R2R3, 3Rand 4R-MYB transcription factors during lignin biosynthesis in flax (Linum usitatissimum), Genomics. (2019), https://doi.org/10.1016/j.ygeno.2019.05.017 (In press). [54] H. Tombuloglu, Y. Slimani, G. Tombuloglu, M. Almessiere, A. Baykal, Uptake and translocation of magnetite (Fe3O4) nanoparticles and its impact on photosynthetic genes in barley (Hordeum vulgare L.), Chemosphere 226 (2019) 110–122. [55] C. Trapnell, B.A. Williams, G. Pertea, A. Mortazavi, G. Kwan, M.J. van Baren, S.L. Salzberg, B.J. Wold, L. Pachter, Transcript assembly and quantification by RNASeq reveals unannotated transcripts and isoform switching during cell differentiation, Nat. Biotechnol. 28 (2010) 511–515, https://doi.org/10.1038/nbt.1621. [56] T. Unver, O. Bozkurt, M.S. Akkaya, Identification of differentially expressed transcripts from leaves of the boron tolerant plant Gypsophila perfoliata L, Plant Cell Rep. 27 (2008) 1411–1422, https://doi.org/10.1007/s00299-008-0560-7. [57] P.M. Waterhouse, R.P. Hellens, Plant biology: coding in non-coding RNAs, Nature 520 (2015) 41–42, https://doi.org/10.1038/nature14378. [58] C. Xie, J. Yuan, H. Li, M. Li, G. Zhao, D. Bu, W. Zhu, W. Wu, R. Chen, Y. Zhao, NONCODEv4: exploring the world of long non-coding RNA genes, Nucleic Acids Res. 42 (2014) D98–103, https://doi.org/10.1093/nar/gkt1222. [59] M. Xin, Y. Wang, Y. Yao, N. Song, Z. Hu, D. Qin, C. Xie, H. Peng, Z. Ni, Q. Sun, Identification and characterization of wheat long non-protein coding RNAs responsive to powdery mildew infection and heat stress by using microarray analysis and SBS sequencing, BMC Plant Biol. 11 (2011) 61, https://doi.org/10.1186/14712229-11-61. [60] H. Yanik, M. Turktas, E. Dundar, P. Hernandez, G. Dorado, T. Unver, Genome-wide identification of alternate bearing-associated microRNAs (miRNAs) in olive (Olea europaea L.), BMC Plant Biol. 13 (2013) 10, https://doi.org/10.1186/1471-222913-10. [61] H. Zhang, Z. Chen, X. Wang, Z. Huang, Z. He, Y. Chen, Long non-coding RNA: a new player in cancer, J. Hematol. Oncol. 6 (2013) 37, https://doi.org/10.1186/17568722-6-37. [62] H. Zhou, Q. Liu, J. Li, D. Jiang, L. Zhou, P. Wu, S. Lu, F. Li, L. Zhu, Z. Liu, L. Chen, Y.G. Liu, C. Zhuang, Photoperiod- and thermo-sensitive genic male sterility in rice are caused by a point mutation in a novel noncoding RNA that produces a small RNA, Cell Res. 22 (2012) 649–660, https://doi.org/10.1038/cr.2012.28.

high-throughput deep sequencing, PLoS One 7 (2012) e50298, , https://doi.org/10. 1371/journal.pone.0050298. J.M. Franco-Zorrilla, A. Valli, M. Todesco, I. Mateos, M.I. Puga, I. Rubio-Somoza, A. Leyva, D. Weigel, J.A. Garcia, J. Paz-Ares, Target mimicry provides a new mechanism for regulation of microRNA activity, Nat. Genet. 39 (2007) 1033–1037, https://doi.org/10.1038/ng2079. S. Griffiths-Jones, miRBase: the microRNA sequence database, Methods Mol. Biol. 342 (2006) 129–138, https://doi.org/10.1385/1-59745-123-1:129. X. Guo, L. Gao, Q. Liao, H. Xiao, X. Ma, X. Yang, H. Luo, G. Zhao, D. Bu, F. Jiao, Q. Shao, R. Chen, Y. Zhao, Long non-coding RNAs function annotation: a global prediction method based on bi-colored networks, Nucleic Acids Res. 41 (2013) e35, https://doi.org/10.1093/nar/gks967. J.B. Heo, S. Sung, Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA, Science 331 (2011) 76–79, https://doi.org/10.1126/science. 1197349. D.R. Hoagland, D.I. Arnon, The water-culture method for growing plants without soil, Circular California Agricultural Experiment Station, 2nd edit, 347 1950. Y. Huang, L. Li, K.P. Smith, G.J. Muehlbauer, Differential transcriptomic responses to Fusarium graminearum infection in two barley quantitative trait loci associated with Fusarium head blight resistance, BMC Genomics 17 (2016) 387, https://doi. org/10.1186/s12864-016-2716-0. B. Inal, M. Turktas, H. Eren, E. Ilhan, S. Okay, M. Atak, M. Erayman, T. Unver, Genome-wide fungal stress responsive miRNA expression in wheat, Planta 240 (2014) 1287–1298, https://doi.org/10.1007/s00425-014-2153-8. L. Juan, G. Wang, M. Radovich, B.P. Schneider, S.E. Clare, Y. Wang, Y. Liu, Potential roles of microRNAs in regulating long intergenic noncoding RNAs, BMC Med. Genet. 6 (Suppl. 1) (2013) S7, https://doi.org/10.1186/1755-8794-6-S1-S7. G. Karakulah, K. Yucebilgili Kurtoglu, T. Unver, PeTMbase: a database of plant endogenous target mimics (eTMs), PLoS One 11 (2016) e0167698, , https://doi. org/10.1371/journal.pone.0167698. A. Keniry, D. Oxley, P. Monnier, M. Kyba, L. Dandolo, G. Smits, W. Reik, The H19 lincRNA is a developmental reservoir of miR-675 that suppresses growth and Igf1r, Nat. Cell Biol. 14 (2012) 659–665, https://doi.org/10.1038/ncb2521. P.J. Kersey, J.E. Allen, I. Armean, S. Boddu, B.J. Bolt, D. Carvalho-Silva, M. Christensen, P. Davis, L.J. Falin, C. Grabmueller, J. Humphrey, A. Kerhornou, J. Khobova, N.K. Aranganathan, N. Langridge, E. Lowy, M.D. McDowall, U. Maheswari, M. Nuhn, C.K. Ong, B. Overduin, M. Paulini, H. Pedro, E. Perry, G. Spudich, E. Tapanari, B. Walts, G. Williams, M. Tello-Ruiz, J. Stein, S. Wei, D. Ware, D.M. Bolser, K.L. Howe, E. Kulesha, D. Lawson, G. Maslen, D.M. Staines, Ensembl Genomes 2016: more genomes, more complexity, Nucleic Acids Res. 44 (2016) D574–D580, https://doi.org/10.1093/nar/gkv1209. D. Kim, G. Pertea, C. Trapnell, H. Pimentel, R. Kelley, S.L. Salzberg, TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions, Genome Biol. 14 (2013) R36, https://doi.org/10.1186/gb-2013-14-4-r36. D. Lauressergues, J.M. Couzigou, H.S. Clemente, Y. Martinez, C. Dunand, G. Becard, J.P. Combier, Primary transcripts of microRNAs encode regulatory peptides, Nature 520 (2015) 90–93, https://doi.org/10.1038/nature14346. L. Li, S.R. Eichten, R. Shimizu, K. Petsch, C.T. Yeh, W. Wu, A.M. Chettoor, S.A. Givan, R.A. Cole, J.E. Fowler, M.M. Evans, M.J. Scanlon, J. Yu, P.S. Schnable, M.C. Timmermans, N.M. Springer, G.J. Muehlbauer, Genome-wide discovery and characterization of maize long non-coding RNAs, Genome Biol. 15 (2014) R40, https://doi.org/10.1186/gb-2014-15-2-r40. J. Liu, C. Jung, J. Xu, H. Wang, S. Deng, L. Bernad, C. Arenas-Huertero, N.H. Chua, Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis, Plant Cell 24 (2012) 4333–4345, https://doi.org/10.1105/tpc.112. 102855. X. Liu, L. Hao, D. Li, L. Zhu, S. Hu, Long non-coding RNAs and their biological roles in plants, Genom Proteom Bioinform 13 (2015) 137–147, https://doi.org/10.1016/ j.gpb.2015.02.003. T.R. Mercer, D.J. Gerhardt, M.E. Dinger, J. Crawford, C. Trapnell, J.A. Jeddeloh, J.S. Mattick, J.L. Rinn, Targeted RNA sequencing reveals the deep complexity of the human transcriptome, Nat. Biotechnol. 30 (2011) 99–104, https://doi.org/10. 1038/nbt.2024. K. Miwa, T. Fujiwara, Boron transport in plants: co-ordinated regulation of transporters, Ann. Bot. 105 (2010) 1103–1108, https://doi.org/10.1093/aob/mcq044. E.P. Nawrocki, D.L. Kolbe, S.R. Eddy, Infernal 1.0: inference of RNA alignments, Bioinformatics 25 (2009) 1335–1337, https://doi.org/10.1093/bioinformatics/ btp157. E.P. Nawrocki, S.W. Burge, A. Bateman, J. Daub, R.Y. Eberhardt, S.R. Eddy, E.W. Floden, P.P. Gardner, T.A. Jones, J. Tate, R.D. Finn, Rfam 12.0: updates to the RNA families database, Nucleic Acids Res. 43 (2015) D130–D137, https://doi.org/ 10.1093/nar/gku1063. E. Ozhuner, V. Eldem, A. Ipek, S. Okay, S. Sakcali, B. Zhang, H. Boke, T. Unver, Boron stress responsive microRNAs and their targets in barley, PLoS One 8 (2013) e59543, , https://doi.org/10.1371/journal.pone.0059543. A. Paytuvi Gallart, A. Hermoso Pulido, I. Anzar Martinez de Lagran, W. Sanseverino, R. Aiese Cigliano, GREENC: a Wiki-based database of plant lncRNAs, Nucleic Acids Res 44 (2016) D1161–D1166, https://doi.org/10.1093/

9