Genome-wide identification of long non-coding RNAs in tomato plants irradiated by neutrons followed by infection with Tomato yellow leaf curl virus

Neutron irradiation, plant growth, and TYLCV infection

We used seeds of the tomato cultivar ‘Seokwang,’ which is susceptible to TYLCV. Two different seed conditions, dry and presoaked, were subjected to neutron irradiation using the MC_50 cyclotron at the Korea Institute of Radiological and Medical Science (KIRAMS) in Seoul, Korea. The presoaked seeds contain a higher portion of oxygen and hydrogen as compared to dry seeds. In particular, it is known that the oxygen can interact with the neutron during irradiation process producing the rapid reactive oxygen species (ROS). ROS plays a pivotal role as a signaling molecule in plants involved in pathogen defense (Apel & Hirt, 2004). Thus, it will be interesting to examine possible effects of ROS in the plant transcriptome. The tomato seeds were subjected to neutron irradiation by proton bombardment of beryllium at 40 MeV energy and a 20 µA current. Two different irradiation times of 30 min and 90 min were applied. Neutron-irradiated tomato seeds were grown in a growth chamber under a 16 h light/8 h dark illumination time (Cheng & Edwards, 1991). The five-week-old tomato seedlings were infected by TYLCV using a whitefly (Bemisia tabaci) vector. TYLCV infected tomato plants were grown for seven weeks in a growth chamber. The general experimental scheme is drawn in Fig. 1A.

Sample collection and total RNA extraction

Leaf samples from individual tomato mutants showing TYLCV disease symptoms were harvested and immediately frozen in liquid nitrogen. Each sample was examined for TYLCV infection by PCR using TYLCV specific primers. Total RNA was extracted using a Qiagen RNeasy Plant Mini Kit (Qiagen GmbH, Hilden, Germany) following the manufacturer’s instructions. DNase I was used to digest genomic DNA in the extracted RNA. The quantity and quality of RNA were measured by an Eppendorf BioPhotometer (MedWOW Ltd., Istanbul, Turkey) and agarose gel electrophoresis.

Library preparation and RNA sequencing

The poly(A) RNA libraries were prepared using an NEB Next Ultra RNA Library Prep Kit for Illumina (New England BioLabs Inc., Ipswich, MA, USA) according to manufacturer’s instructions. The prepared libraries were analyzed by a 2100 Bioanalyzer instrument (Agilent Genomics, Waldbronn, Germany) to measure quality. A total of eight libraries were paired-end (101 bp × 2) sequenced by an Illumina HiSeq 2000 system at Macrogen, Seoul, Korea. The raw sequenced data from this study were deposited in the SRA database at the National Center for Biotechnology Information with the following accession numbers, SRR6019475 –SRR6019483. The detailed information of samples is provided in Table 1.

Assembly of tomato transcripts using RNA-Seq data

All raw data from eight libraries were aligned using Tophat 2.1.1 to the tomato reference genome (ITAG2.4_genomic.fasta) from the International Tomato Genome Sequencing Project (Trapnell et al., 2012). In addition, the Gene Transfer Format (ITAG2.4_gene_models.gff3) file was used for alignment using Tophat resulting in a single BAM file. The BAM file was subjected to Cufflinks 1.3.0 (Trapnell et al., 2010) to assemble transcripts. Finally, Cuffcompare was used to annotate lncRNAs (Trapnell et al., 2010). We obtained assembled transcripts using Cufflinks. The assembled transcripts were subjected to Cuffcompare to annotate the assembled transcripts by comparing the assembled tomato transcripts to the tomato reference annotation. After that, all transcripts including lncRNAs were annotated. Based on annotation, identified lncRNAs were categorized as intergenic regions, antisense, overlapping, and intronic regions.

Identification of tomato lncRNAs

As shown in Fig. 1B, several steps were used to predict novel lncRNAs from the assembled transcripts. First, we extracted assembled transcripts in a fasta format using the gffread program (Weirick et al., 2015). The assembled transcripts were subjected to the Coding Potential Calculator 2 (CPC2) to predict the coding capacity of each transcript with default parameters (Kang et al., 2017). Only transcripts labelled with “noncoding” predicted by CPC2 were included for further analyses. Second, the assembled transcripts were subjected to BLASTX (version 2.2.31) against a tomato protein database with an evalue of 1e−3 as a cutoff to exclude transcripts with a shared sequence similarity to known tomato proteins. In addition, the transcripts with lengths of less than 200 nt were discarded (Nagata et al., 2004). Third, the transcripts were subjected to PfamScan (version 31.0) with default parameters against the hmmer protein database (version 3.0b2) to identify those containing conserved protein domains (Li et al., 2016). Fourth, results from CPC2, PfamScan, and BLASTX were compared to exclude transcripts with coding capacities. Fifth, we excluded transcripts that showed sequence similarity to known miRNA precursors. For this, BLASTN with an evalue of 1e−as a cutoff was conducted against the plant miRNA precursors, house-keeping ncRNAs, small nucleolar RNAs (snoRNAs), and non-redundant (NR) protein databases (Liu et al., 2005; Pruitt, Tatusova & Maglott, 2005; Xiao et al., 2009). Finally, we calculated the fragments per kilobase of transcript per million mapped reads (FPKM) to identify lncRNAs using the Binary version of Sequence Alignment (BWA) followed by the BBMap program with default parameters (Li & Durbin, 2010; Rodríguez-García, Sola-Landa & Barreiro, 2017). The transcripts in which FPKM values were greater than one were considered putative lncRNAs.

Functional classification of lncRNAs and gene ontology (GO) enrichment analysis

To obtain a broad functional overview of identified lncRNAs, we predicted targets of identified lncRNAs by BLASTN search against tomato mRNA sequences using an e-value of 1e−as a cutoff. Finally, we selected 794 lncRNAs associated with tomato mRNAs. The selected tomato gene sequences were blasted against Arabidopsis genes (https://www.arabidopsis.org/) to obtain known functions. A majority of tomato genes were converted to the corresponding Arabidopsis locus. Using Arabidopsis loci, we conducted GO enrichment analysis using the GOEAST program with default parameters (Zheng & Wang, 2008). The GO terms obtained for tomato genes possessing lncRNAs were classified according to biological process, molecular function, and cellular component. The REVIGO program was used for visualization of enriched GO terms (Supek et al., 2011).

Expression profiles for 1,563 lncRNAs in response to TYLCV infection

To establish the expression profiles for the 1,563 lncRNAs in this study, we used previously published RNA-Seq data (PRJNA291401) (Wang et al., 2015). The published RNA-Seq data consisted of three mock and three TYLCV infected samples. The identified 1,563 lncRNAs were used as reference sequences. Raw sequence reads from the published RNA-Seq data were mapped on the 1,563 lncRNAs using Tophat (Trapnell et al., 2012). Transcript assembly was conducted by Cufflinks. TYLCV infected samples were compared to mock samples to identify differentially expressed genes (DEGs) using Cuffdiff based on a p-value less than 0.05 and a log₂(fold change) greater than one.

Assembly of tomato transcripts using eight RNA-Seq data

To maximize the identification of lncRNAs, we combined all raw data from eight libraries. Raw sequence reads from eight different libraries were subjected to mapping on the reference tomato genome. Of the 128 million reads from eight libraries, almost 117 million (92.2%) were mapped to the tomato genome, whereas about 11 million (8.8%) were not mapped (Table 2). The number of mapped reads ranged from 14,824,307 (30S1) to 16,273,849 (90D2), whereas the portion of unmapped reads ranged from 7.4% (90D2) to 12.0% (30S1). These mapped reads were used for transcript assembly by Cufflinks, resulting in a total of 39,067 transcripts (Fig. 1B).

Identification of tomato lncRNAs

To identify putative tomato lncRNAs, nine different processes were applied (Fig. 1B). First, 19,780 (50.63%) transcripts with coding capacity were excluded by the CPC2 program, and then a BLASTX search against tomato proteins excluded 17,029 (43.58%) transcripts. In total, 165 transcripts with lengths less than 200 nt were excluded. PfamScan was conducted to exclude three transcripts encoding conserved protein domains. This was followed by a BLASTN search against a miRNA precursor database, which excluded 13 transcripts. In addition, 89 and 10 transcripts were excluded by a BLASTN search for showing sequence similarity to known ncRNAs and snoRNAs, respectively. Finally, the remaining transcripts were subjected to a BLASTX search against the NR protein database in NCBI, which excluded 553 transcripts. As a result, we obtained 1,563 lncRNAs from eight tomato samples (Table S1).

Classification of identified tomato lncRNAs

The lengths of the identified lncRNAs ranged from 201 nt to 4,647 nt. About 26% (416 lncRNAs) of identified lncRNAs were more than 1,000 nt, while 74% of identified lncRNAs were less than 1,000 nt. However, most protein coding transcripts (94.3%) were less than 1,000 nt in size. Among the lncRNAs less than 1,000 nt, lengths ranged from 301 to 800 nt (Fig. 2A). We compared the length distribution of identified lncRNAs between our study and a previous study (Wang et al., 2015). Both studies showed similar length distribution of identified lncRNAs. In particular, the length of lncRNAs ranged from 400 nt to 500 nt was the highest proportion in both studies. In addition, we examined the number of exons in lncRNAs as well as protein coding transcripts (Fig. 2B). Most lncRNAs (90.2%) were derived from a single exon, and there were two lncRNAs, TCONS_00021506 and TCONS_00021507, derived from six and seven exons, respectively. However, the functions of corresponding genes were unknown. In comparison, half of the protein coding transcripts contain a single exon, whereas 11.7% of the protein coding transcripts have more than 10 exons.

Next, we categorized the identified lncRNAs. Most lncRNAs (49%) were derived from intergenic regions, followed by antisense lncRNAs (35%), overlapping (4%), and intronic regions (3%) (Fig. 2C). Apart from the 769 lncRNAs derived from intergenic regions, 794 lncRNAs were associated with a gene. We further classified the 794 lncRNAs according to chromosome. Five lncRNAs were not assigned to any chromosome. It seems that they were located on random scaffold. With the exception of 110 lncRNAs on chromosome 1, the number of lncRNAs on each chromosome ranged from 50 (Chromosome 12) to 77 (Chromosome 3), as shown in Fig. 2D. The positions of identified lncRNAs on each tomato chromosome are indicated in a circos diagram (Fig. 2E). In order to visualize the distribution of identified lncRNAs on different chromosomes, a combination of different graphs including bar plots and volcano plots were used. Interestingly, most lncRNAs were highly enriched at the beginning and the ending of each chromosome. There were relatively few lncRNAs located in the middle of individual chromosome.

Expression profiles of identified lncRNAs in eight samples

We examined the expression of 1,563 lncRNAs in eight different samples by calculating the FPKM values. Some lncRNAs were not expressed in eight conditions; however, most lncRNAs were highly expressed. In order to examine the distribution of the expression values for all identified lncRNAs, we generated a box plot using log₁₀ converted FPKM values (Fig. 3A). Although standard deviation showed a high degree of difference in each sample, minimum, median, and maximum values among samples did not show significant difference. Next, we conducted PCA analysis to cluster eight conditions based on gene expression of lncRNAs. Our result showed that seven conditions except the condition for soaked sample for 30 min (30S1) were clustered together (Fig. 3B).

We next examined the expression profile of 1,563 lncRNAs in response to TYLCV infection using previously published data (Wang et al., 2015). As a result, 915 lncRNAs out of 1,563 were expressed in mock and TYLCV infected sample (Table S2). A volcano plot showed that the number of up-regulated lncRNAs (19 lncRNAs) was higher than that of down-regulated lncRNAs (11 lncRNAs) (Fig. 3C). Expression of 30 differentially expressed lncRNAs was visualized by a heat map (Fig. 3D). Unfortunately, the functions of those 30 differentially expressed lncRNAs are not currently known due to lack of functional studies associated with tomato lncRNAs.

Functional classification of tomato genes associated with identified lncRNAs

Of the 1,563 lncRNAs identified, 794 were associated with tomato genes (Table S3). In order to reveal the functional roles of tomato genes associated with lncRNAs, we examined GO using corresponding Arabidopsis genes. GO enrichment analyses revealed that 47 GO terms (biological process), eight GO terms (cellular component), and 14 GO terms (molecular function) were highly enriched in tomato gene associated lncRNAs (Table S4). In relation to biological processes, GO terms associated with positive regulation of the abscisic acid-activated signaling pathway, organic substance metabolism, gluconeogenesis, organelle organization, and cellular catabolism were highly enriched (Fig. 3E). Based on cellular component analysis, GO terms related to cytosol, membrane, and cell junction were highly enriched (Fig. 3F). Interestingly, according to molecular function, many GO terms were associated with binding. For example, heterocyclic compound binding, phospholipase D activity, and ubiquitin protein ligase binding were frequently identified (Fig. 3G).

Identification of target tomato mRNAs of the 1,563 lncRNAs

GO enrichment analysis revealed that many genes associated with lncRNAs have binding functions important for interaction of lncRNAs with other molecules such as RNAs and proteins. We identified target tomato mRNAs that showed strong sequence similarity to the lncRNAs identified by BLASTN search (Table S5).

Of the lncRNAs identified, 824 showed sequence similarity to at least one tomato mRNA. For instance, 566 lncRNAs has a single target, whereas two lncRNAs showed sequence similarity to 10 different nucleotide sequences (Fig. 4A). The lncRNA TCONS_00000794 (Fig. 4B), which was 2,297 nt in length, displayed 10 different regions of six genes. The lncRNA TCONS_00003273, which was 2,081 nt in length, displayed 10 different regions of six genes (Fig. 4C).

Of 30 differentially expressed lncRNAs, six lncRNAs showed sequence similarity to tomato mRNAs; however, functions of only three corresponding mRNAs are known. For instance, the lncRNA TCONS_00005642, which was down-regulated by TYLCV infection, is associated with BHLH transcription factor (Solyc02g063440.2). Two lncRNAs, TCONS_00020980 and TCONS_00035472, showed sequence similarity to a gene coding pyruvate decarboxylase 2 (Solyc06g082140.2) and a gene coding ariadne-like ubiquitin ligase (Solyc11g008590.1), respectively. Both lncRNAs were up-regulated by TYLCV infection.