Genome-wide identification of Arabidopsis long noncoding RNAs in response to the blue light

To annotate the relationships between lncRNAs and miRNAs in Arabidopsis under blue light, we performed lncRNA sequencing of three biological replicates, each with three samples using 4 day dark-grown WT seedlings exposed to continuous blue light for 0 h, 2 h and 8 h (Fig. 1A). 1,095,370–1,100,180 M paired-end reads with 100 bp in length for the mRNAs and lncRNAs, and 11,600,936–12,413,799 single-end reads with 50 bp in length for the miRNAs were acquired from nine samples (Table S1 in Supporting Information). All the clean reads were aligned to Arabidopsis reference genome assembly (UMD3.1.80) for further analysis. We calculated the correlation coefficient (r²) of the sequencing data among the three individual samples of biological replicates in each period based on the FPKM and TPM values, and found that the r² values were 0.8499–0.9958, 0.9774–0.9995 and 0.8955–0.9907 for the mRNAs, lncRNAs and miRNAs, respectively, indicating that the similarities of the mRNA, lncRNA or miRNA biological replicates were sufficiently high (Fig. S1).

To annotate lncRNAs, we extracted the total RNA from samples. Non-coding RNAs and mRNAs were enriched by removing the ribosomal RNAs (rRNAs) from the total RNA. Due to the instability of rRNA removing efficiency, we also removed rRNA-contained reads by alignment. We choose Coding Poteintial Calculator (CPC)²⁸ and iSeeRNA²⁹ which perform well compared to other softwares in both accuracy and efficiency for lncRNA prediction. Cuffmege was used to merge several assemblies together and filter the transfrags that are probably artifacts. A single annotation file was then generated for downstream differential analysis. We obtained the whole parsimonious set of transcripts to detect the novel transcripts from the initial assemblies, and compared the assembled transcripts to the reference annotation. HISAT reference annotation based assembly method was used to reconstruct the transcripts, and background noise was reduced by using FPKM and coverage threshold^30,31. The novel transcripts were identified by calculating the coding potential of these transcripts (Fig. 1B).

We also calculated the distribution of lncRNA loci induced under blue light treatments for 2 h and 8 h in NONCODE v3.0³² with Cuffcompare categories. We found that the majority of lncRNA loci (46.15%) are antisense, while 28.85% are intergenic, 10.58% divergent, 9.62% intronic and 4.81% sense (Fig. 1C).

To identify the regulatory networks of mRNAs, lncRNAs and miRNAs under blue light, we analyzed the differentially expressed mRNAs, lncRNAs and miRNAs (DEmRNAs, DElncRNAs, DEmiRNAs). Compared to the dark condition, 481 or 545 DElncRNAs were found in Arabidopsis seedlings under blue light treatment for 2 h or 8 h respectively, (Fig. 2A,B), of which 316 or 308 were up-regulated and 165 or 237 down-regulated. In addition, 4197 or 5207 DEmRNAs were identified in Arabidopsis seedlings under blue light treatment for 2 h or 8 h respectively, with 3222 or 3970 up-regulated and 975 or 1237 down-regulated more than two fold. In total, 375 or 286 DEmiRNAs were found in Arabidopsis seedlings under blue light treatments for 2 h or 8 h respectively, with 233 or 175 up-regulated and 142 or 111 down-regulated for more than two fold (Fig. 2C,D).

To validate the expression patterns of lncRNAs under the blue light treatment, we performed qRT-PCR analysis to examine the expressions and sequences of 6 lncRNAs (which include 4 transcripts present in the Arabidopsis database (http://www.noncode.org/↗) and 2 novel transcripts based on our assembly of RNA-seq data) (Fig. 3A–F). The expression patterns of all lncRNAs were found to be consistent with the RNA-seq data. In addition, we also performed Northern blots of 4 miRNAs to validate the results of miRNA sequencing. We found that the expression patterns of the 4 miRNAs were similar to those shown in the miRNA sequencing data (Fig. 3G).

It was known that lncRNA can play a role in sequestering the specific miRNA by acting as a ceRNA in a type of target mimicry to protect the target mRNA in both plants and animals^33–36. We measured the likelihood of the ceRNAs using StarBase v2.0^37,38. Using the hypergeometric test³⁹, we calculated the p-value for each potential ceRNA pair (lncRNA-mRNA)³⁸. Our bioinformatics analysis indicated that 180 of the identified lncRNAs have the potential to compete with miRNAs, including several conserved miRNAs such as miR160, miR164, miR167, miR169 and miR172 (Fig. 4A,B, Figs. S2 and S3). We then constructed lncRNA-miRNA-mRNA interaction networks (Figs. S2 and S3). We found that 744 DEGs (differentially expressed genes) and 183 DElncRNAs were targeted by 155 miRNAs (Figs. S2 and S3).

We found that miR167 plays a role in modulating hypocotyl elongation²⁷, therefore, we focused on the regulatory network of lncRNA-miRNA167-mRNA. The 14 complementary nucleotides of miR167 to BLIL1 are shown in Fig. 4, including 8 bases at the 5′ end of miR167 (Fig. 4B,C). To confirm the expression profiles of BLIL1 and miR167, we performed qRT-PCRs to examine the transcript levels of BLIL1 and miR167-targeted mRNAs under blue light treatments for 2 h, 8 h and 12 h. As shown in Fig. 5A–C, the transcript levels of BLIL1, ARF6 and ARF8 which are targets of miR167⁴⁰, were induced under blue light treatments. Moreover, Northern blots confirmed that the transcript level of miR167 is opposite to those of BLIL1, ARF6 and ARF8 (Fig. 5A–C).

We then investigated the role of the BLIL1-ARF6/8-miR167 regulatory network in photomorphogenesis by examining the blue light-dependent inhibition of hypocotyl elongation. We analyzed the hypocotyl growth of 4-day-old Col and blil1-1, blil1-2, arf6/8 and miR167 mutant seedlings grown under continuous blue light at 10 µmol s⁻¹m⁻² intensity, or continuous darkness. The blil1-1, blil1-2, arf6/8 mutants display long hypocotyl phenotypes compared with wild type after treatment with continuous blue light, and miR167 mutant shows shorter hypocotyl than Col under blue light (Fig. 5D–F). The hypocotyl lengths of these mutants were similar to that of Col grown under dark condition (Fig. S4). These results revealed the competitive relationships between miR167 and BLIL1 or ARF6/8 in regulation of the blue light-dependent photomorphogenesis.

As miR167 was reported to be involved in stress responses⁴¹, we validated the competitive relationship among BLIL1, ARF6/8 and miR167 under osmotic stresses. We performed qRT-PCRs to examine the induction profile of BLIL1, ARF6/8 and miR167 under 200 mM mannitol stress. When the seedlings were treated by 200 mM mannitol, the levels of BLIL1 and ARF6/8 transcripts increased slightly at 3 h, obviously at 6 h, and reached the maximum at 12 h (Fig. 6A,C), whereas the level of miR167 decreased obviously at 6 h with no significant difference at 3 h (Fig. 6B).

To further confirm the competitive relationship, we investigated the mannitol tolerances of blil1-1, blil1-2, arf6/8 and miR167 mutant seedlings. Five days after germination on normal MS medium, the seedlings were transferred to MS medium plates containing 200 mM mannitol, and the plates were stood-up. The amounts of lateral roots were measured at the 7th day after transferring. As shown in Fig. 6D–F, the results demonstrated that blil1-1, blil1-2, arf6/8 mutants displayed more lateral roots than Col, while miR167 mutant showed less lateral roots than Col under mannitol stress. As a control, no significant differences were observed among the mutants and Col seedlings under normal MS medium. These results indicated that the mutants of blil1-1, blil1-2, and arf6/8 showed reduced sensitivity to mannitol treatment, whereas miR167 displayed elevated sensitivity to mannitol treatment, suggesting that BLIL1 and ARF6/8 acts as positive factors, and miR167 as a negative factor in response to mannitol stress.

Arabidopsis thaliana (ecotype Col-0), miR167b (CS872594↗),arf6 (CS24606), arf8 (CS24608), arf6/8 (CS24632), blil1-1 (SALK_112856) and blil1-2 (CS65633) mutants were obtained from the Arabidopsis Biological Resource Center (Ohio State University). All plants were grown in soil or Murashige and Skoog (MS) medium at 16 hr light/8 hr dark photoperiod unless specifically indicated otherwise.

The in-house Perl script was used to process the original reads in fastq format. Clean reads were obtained by removing the reads with adapters and poly-N, and low-quality reads. Q20 (base ratio > 20), q30 (base ratio > 30) and GC content of the obtained clean data were calculated. The high-quality clean data were used for all the subsequent analysis⁶⁰. The Arabidopsis reference genome and gene model annotation files were downloaded from ftp://ftp.ensembl.org/pub/release-80. We built an index of the reference genome using Bowtie v2.0.680 and aligned paired-end clean reads to the reference genome using TopHat v2.0.981. The mapped reads of each sample were assembled by Stringtie⁶¹.

Fastx_uutoolkit-0.0.13.2 was used to pre-process the original data by removing low-quality reads, 5′ joint contaminated reads, 3′ spliced reads, and reads with length less than 18 nt or with poly-A. The clean reads were then compared with known ncRNAs by using blast (v2.2.26) in RFAM database (v10.1). Clean reads matching any sequence in the RFAM database was deleted. Then, the proportion of four bases (A/T/G/C) in the clean reads at each location was calculated. The Arabidopsis miRNA sequences were downloaded from miRBase (http://www.mirbase.org/↗), and used to search for the clean reads consistent with the Arabidopsis reference genome⁶⁰.

We used CuffDiff (V2.1.1) to calculate the values of lncRNAs and coding mRNAs for FPKM (the number of fragments per kilo-base of exon per million mapped fragments), and TPM (transcripts per million) fractions of miRNAs. The differential expressions of miRNAs were analyzed by EdgeR. We analyzed differential expressions of lncRNAs, mRNAs and miRNAs by comparing blue light treatment with dark condition. In any pairwise comparisons, the expressions of lncRNAs and mRNAs with corrected p-values < 0.05 and absolute fold-change values>2.0, or miRNAs with p-values < 0.05, absolute fold-change values>2.0, and miRDeep2 score ≥1, were considered to be differentially expressed⁶⁰.

The first strand of cDNA is synthesized from total RNA (1 μg) and M-MLV reverse transcriptase (Promega) using random primers and serves as a template for subsequent PCR amplification. Pri-miRNA levels were detected by qRT-PCRs using Bio-Rad CFX Real-Time System. Actin gene was used as the internal control for normalization of DNA templates. Each PCR was repeated at least three times. Data were analyzed by Bio-Rad iCycler iQ Real-Time Detection system. The transcript level was calculated by the relative 2^–ΔΔCt method⁶². The primers used for qRT-PCRs are listed in Table S2 supplemental information.

Small RNAs were isolated from 4-day-old plant seedlings by using the mirVana miRNA Isolation Kit (Ambion, AM1561) and sequenced by high-throughput sequencing of Illumina HiSeq. 2000. Sequencing data of all known miRNAs were hierarchical clustered in an unsupervised manner to analyze the degree of differential expression of small RNA⁶³.

LncRNAs and mRNAs were isolated from 4-day-old plant seedlings. Ribosomal RNAs (rRNAs) were removed using TruSeq Stranded Total RNA with Ribo-Zero Plant kit (Illumina, 20020610) and the rRNA-free residue was obtained by ethanol precipitation. The first-strand cDNA was synthesized using a random Hexamer-primer and reverse transcriptase (Invitrogen). The second strand cDNA was synthesized using RNase H (Invitrogen) and DNA polymerase I (New England BioLabs). Subsequently, the sequencing libraries were generated using the rRNA-free RNA with a NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, MA, USA), and sequenced by high-throughput sequencing of Illumina Hiseq Xten pair-end 150 (PE150).

The seedlings were treated with light as described⁶⁴ with some modifications. The seedlings of each line were cultured in the same 100 mm Petri dish for each repeat. The seeds were soaked in 70% ethanol for 5 minutes, and then in 95% ethanol for 5 minutes. Seeds were vernalizated for 3 days in darkness at 4 °C. The culture dishes were placed in white light for 3 hours to induce germination and in the growth chamber at 21 °C for 21 hours. The dishes were then transferred to experimental light conditions. The seedlings were transplanted into blue light (0.1, 0.5, 1 and 10 µmol s⁻¹m⁻²) in a growth chamber (Percival Scientific Inc., Perry, Iowa, USA) for 4 days at 21 °C, respectively. Dark control seedlings were kept in darkness. The hypocotyl lengths of at least 30 seedlings were measured after blue light or dark treatments.

Five days after germination on normal MS medium, the wild type and blil1-1, blil1-2, arf6/8 and miR167 mutant seedlings were transferred to MS medium plates with or without 200 mM mannitol, and the plates were stood-up. The phenotypes of seedlings were visualized 7 days after transferring.

Small RNA gel blot was performed as previously described²⁴. Small RNAs were isolated from 4-day-old plant seedlings by using mirVana miRNA Isolation Kit (Ambion, AM1561). Small RNAs of about three micrograms were fractionated on a 15% polyacrylamide gel containing 8 M urea and then transferred to a nylon transfer membrane (GE Healthcare). Antisense oligonucleotides (Table S2 Supplemental information) were synthesized and biotin probes were labeled at 3′ end. Hybridization was carried out overnight in Ambion (AM8663) at 42 °C. A probe complementary to U6 (5′CATCCTTGCGCAGGGGCCA 3′) was used as a loading control.

All necessary data generated or analyzed during the present study are included in this published article, its Supplementary Information files and NCBI (the National Center for Biotechnology Information) with the number PRJNA596364.