Strand-specific RNA-seq reveals widespread occurrence of novel cis- natural antisense transcripts in rice

ssRNA-seq and assembly of rice transcripts

In order to comprehensively identify rice cis- NATs, we applied ssRNA-seq to deeply sequence rice cDNAs for assessing transcripts with clear transcriptional orientations. Rice seedlings grown under normal and three abiotic stress conditions (salt, cold and drought treatments) were collected for preparations of mRNA and small RNAs (see Methods), and four strand-specific cDNA libraries were then constructed according to the empirical protocol [37]. In general, a modified RNA-seq method is used for incorporation of deoxy-UTP during second-strand cDNA synthesis and subsequent destruction of the uridine-containing strand in the sequencing library. This enables the identification of transcript orientation. The high-throughput sequencing of strand-specific cDNAs was performed on the Illumina GAIIX. In total, 14.7, 11.5, 14.2 and 13.2 million paired-end reads of 2 × 120 bp with high-quality, which were generated from the untreated and salt, cold and drought stress treatments libraries, respectively, matched unambiguously and uniquely to the rice reference genome [41] (Additional file 1). We estimated the accuracy of transcriptional orientation by comparing the sequencing reads with the annotated gene datasets [42]. About 89.4-95.5% of the mapped reads appeared to be aligned with correct transcriptional orientation, providing strong evidence for ssRNA-seq data in a strand-specific manner. In contrast to ssRNA-seq data, only half of the mapped reads from the previous RNA-seq data were consistent with the gene models at the same transcriptional orientations. In addition, we also deeply sequenced epidermal cells of normal rice seedlings, and about 10.4 million paired-end reads of 2 × 100 bp were generated as high-quality data.

Using the software TopHat and Cufflinks [43–46], we assembled ssRNA-seq data combined with the rice genome annotation [42] into 76,013 transcripts corresponding to 45,844 unique gene loci (including 4873 novel gene loci). Of them, 25,924 were identified as novel transcripts, which were composed of 5063 ncRNA, 16,494 CDS with protein hits and 4367 CDS without any protein hits (Additional file 2).

Identification of putative cis-NATs in rice

Small RNAs and nat-siRNAs in rice

To investigate the complexity of small RNAs in rice, we generated four small RNA libraries from rice seedlings under normal conditions and three abiotic stress treatments (salt, cold and drought), and deeply sequenced the libraries on the Illumina GAIIX. The small RNA libraries were made with an RNA ligation method which produces strand-specific libraries. After removing low-quality reads and those mapped to rRNA, tRNA, sn/snoRNA, mitochondria and chloroplasts, 4,843,040, 3,973,627, 2,894,255 and 5,492,145 distinct small RNAs (representing 17,632,759, 12,923,509, 8,720,251 and 20,069,157 of 48,683,191, 49,254,272, 25,705,840 and 50,192,805 raw reads, respectively) from corresponding untreated, salt, cold and drought conditions, were identified (Additional file 1). Small RNAs of 24-nt were the predominant size class (Additional file 4A). We found that the majority were located in transposon-related regions, followed by upstream, intergenic, downstream, and intron regions (Additional file 4C).

To determine the amount of nat-siRNAs, we searched in the four small RNA libraries that matched uniquely and perfectly to the overlapping regions of the 3819 cis- NATs. A total of 180,239 reads corresponding to 90,977 unique small RNAs (normal: 25,420; salt: 18,598; cold: 18,152 and drought: 28,807) were derived from overlapped regions (Additional files 1 and 5). Here, nat-siRNAs showed wide ranges in size of 18-34 nt, with 21-25 nt the most common size (Additional file 4B). The 5^′-first nucleotide of nat-siRNAs was predominantly adenosine (Additional file 4F), which differed from that of the total small RNAs (Additional file 4E).

cis-NATs and nat-siRNAs with gene expression evidence under normal or abiotic stress conditions

Of 2292 cis- NATs, 1072 (46.7%) cis- NATs were expressed in both normal and abiotic stresses (Figure 1A). The cis-NAT pairs could be divided into five subgroups based on a scatter plot comparing transcripts expressional ratio trends of cis- NAT pairs between normal and cold stressed conditions (Figure 1B). We also investigated the functional bias among each group of cis- NAT pairs through functional domain and expressional profiling analyses. We found that Protein kinase domains were commonly identified in cis- NAT pairs among the five subgroups. However, other functional domains and expression levels of sense/antisense transcripts were detected to be associated with different subgroups of cis-NAT pairs.

Subgroup-1 cis- NATs were predominantly associated with the pentatricopeptide repeat (PPR) and Protein tyrosine kinase domains. We also found that the expression levels of the cis- NAT pairs of Subgroup-1 under stress conditions were similar to that under normal conditions. Subgroup-2 cis- NATs were identified to contain domains of leucine rich repeats and glycosyltransferase. In this subgroup, the expression levels of sense-transcripts of the cis- NAT pairs were significantly increased under stress conditions, while expression levels of antisense transcripts were greatly reduced. Thus, expression levels were higher for sense than for antisense transcripts under stress conditions. Subgroup-3 can be classified as the family of alpha/beta hydrolase fold. The expression levels of sense-transcripts of the Subgroup-3 cis- NAT pairs were significantly reduced under stress conditions, while expression levels of antisense transcripts significantly increased. Thus, expression levels were lower for sense than for antisense transcripts under stress conditions. Subgroup-4 cis- NATs can be classified as the families of the eukaryotic aspartyl protease and sugar (and other) transporter. Subgroup-5 cis- NATs were identified to be the families of glycosyl hydrolases and ubiquitin-conjugating enzyme. The remaining four cis-NAT pairs were not classified into any subgroups, as the expression levels of sense transcripts were not correlated with that of anti-sense partner under normal and abiotic conditions. We identified these four cis- NAT pairs as abnormal values (Figure 1B; indicated as the brown spots). Subgroup-1 was shown to be represented the major part. Ratio values of 913 cis- NAT pairs (Subgroup-1) between normal and cold stressed conditions were highly correlated, with R² > 0.85. Similar results were also observed between normal and salt/drought stressed conditions with R² > 0.93 of 917 pairs and R² > 0.93 of 898 pairs, respectively (Additional file 10). In total, of 1072 cis- NAT pairs, 95.5% (1024) showed positive correlations between sense and antisense transcripts. One cis- NAT pair (Os05t0500000-00 vs. Os05t0500101-01) was taken as an example to compare expression changes under normal and cold stressed conditions (Figure 2). Northern blot analyses also demonstrated their differential expression levels under abiotic stresses (Figure 2C). Sense transcript Os05t0500000-00 is annotated as an UDP-glucoronosyl and an UDP-glucosyl transferase, which is an important enzyme for catalyzing transportation of sugars. Antisense transcript Os05t0500101-01 had no hits to any protein families.

The numbers of cis-NATs expressed under different conditions were quite different (Figure 1A). The shared cis-NATs expressed under normal (control), cold, salt and drought conditions were also identified. Among them, 503 (21.9%) cis-NATs were detected preferentially under abiotic stresses compared to normal conditions. There were 84, 74 and 128 cis-NAT pairs clearly expressed under drought, salt and cold stresses, respectively. As an example, under cold treatment, sense and antisense transcripts of cis- NAT pair Os09t0482800-02 vs. CUFF.14823.1 were both much more expressed than under normal or other stressed conditions (Figure 3). Semi-quantitative RT-PCR and reverse transcription PCR (RT-PCR) were used to validate the expression of the cis-NAT pair under normal and cold conditions (Figure 3B). Functional annotation of this cis- NAT pair was composed of EF-hand protein (Os09t0482800-02) and ncRNA (CUFF.14823.1). Moreover, 1238 of 2292 (54.0%) cis- NATs were expressed in leaf epidermal cells. Of 1238 cis- NATs, 725 belonged to co-expressed pairs (Table 2). We further performed analysis of differentially expressed genes (DEGs). In total, we identified 112 cis- NATs with DEGs of at least either of each pair, which belonged to different subgroups (Additional file 11). Of them, expression levels of 69 cis- NATs (either sense or antisense transcripts) were up-regulated under abiotic stresses, while 46 cis- NATs were down-regulated under abiotic stresses. Moreover, the results showed no clear expression pattern between transcripts and nat-siRNAs (Figure 4). The number of nat-siRNAs notably increased along with higher transcriptional expression level under cold-stress conditions (Figure 4A and B); in contrast, some cis- NATs generated less nat-siRNAs, along with higher expression levels (Figure 4C and D). In addition, some cis- NATs generated nat-siRNAs with more complex expression patterns (Figure 4E-H).

We used an enrichment analysis of PFAM protein families [47] to functionally characterize three sets that could be annotated from 2292 cis- NATs. The three sets were 767 cis- NAT pairs of protein-coding genes and non-protein-coding RNAs (Set I), 818 cis- NAT pairs of protein-coding genes and predicted CDS without any PFAM domain-containing (Set II), and 377 cis- NATs composed of both protein-coding genes (Set III), respectively. Roughly, 43 protein families were the major enriched terms in either of three sets (Figure 1C). Of them, four proteins exhibited distinctively enriched terms with P-values < 0.001. PPR families were enriched in Set III, leucine rich repeat and NB-ARC domain were significantly reduced in Set III, and UDP-glucoronosyl and UDP-glucosyl transferase were significantly enriched in Set II (Figure 1C).

Networks formed by cis-NAT groups

We further investigated networks formed by cis- NATs in rice, named ‘many-to-many’ cis- NATs, i.e. one transcript in a cis- NAT had more than one antisense partner. Networks formed by cis-NATs reflect the complexity of their post-transcriptional regulation. It has been reported that in Arabidopsis, several genes are involved in two cis- NATs, one is convergent and the other divergent [20, 27]. In rice, we found 461 cis- NATs (composed of 685 transcripts) involved in 223 cis- NAT groups (Additional file 12). Of them, 209 groups belonged to one-to-two type, nine belonged to one-to-three type, four belonged to two-to-two groups, and one belonged to one-to-four group, respectively. Five types were identified in one-to-two cis- NAT groups (Additional file 12). Except 46 groups which were composed of convergent together with divergent cis- NATs, it also included 74 enclosed together with enclosed, 40 enclosed together with divergent, 48 enclosed together with convergent, and one convergent together with convergent. Interestingly, we found only 17.4-34.7% cis- NAT groups could be detected with simultaneous expression evidence and nat-siRNAs in the same sample from normal, cold, salt or drought conditions; another 29.2-48.6% cis- NAT groups were either expressed (Additional file 13). Similar to one-to-one type cis- NATs, protein kinase domain, leucine rich repeats, NB-ARC domain and the PPR repeat family were the major enriched PFAM protein domains.

Materials and growth conditions

Rice seeds (O. sativa ssp. japonica cv. Nipponbare) were used in all experiments. Seeds germinated at 28°C in darkness for 2 d were transferred to a plant growth chamber to grow for 14 d under controlled conditions (12/12 h and 30/24°C of light/dark cycles) to produce seedlings and epidermal cells. For drought and salt stress treatments, 14-d-old seedlings were incubated in solutions containing 20% PEG-6000 and 200 mM NaCl, for 4 h at 30°C, respectively. For cold treatment, seedlings at the same developmental stage were treated at 4°C for 24 h in darkness.

Strand-specific cDNA library construction and sequencing

We prepared the strand-specific cDNA libraries according to a protocol [37]. The ssRNA-seq is a simple modification of the RNA-Seq method that incorporates deoxy-UTP during second-strand cDNA synthesis and subsequent destruction of the uridine-containing strand in the sequencing library. Thus it enables identifying the orientation of transcripts. Total RNA was isolated using the TRIzol reagent (Invitrogen), then total genomic DNA was removed from tissues using DNase (New England Biolabs), which was examined by gel electrophoresis. The OligoTex mRNA midi kit (Qiagen) was used to purify poly(A) mRNA from the total RNA samples. Next, mRNAs were fragmented using the RNA fragmentation kit (Ambion). The first cDNA strand was synthesized using random hexamer primers and second-strand cDNA was synthesized where dUTP was used instead of dTTP. In this step, Actinomycin D was used to increase strand specificity by inhibiting second-strand cDNA synthesis. At 15°C 0.5 μl of actinomycin D solution (120 ng/μl), 0.5 μl of RNase OUT (40 units/μl, Invitrogen) and 0.5 μl of SuperScript III polymerase (200 units/μl,Invitrogen) were added to the reaction. Then EB (20 μl) (10mM Tris–Cl, pH 8.5, Qiagen) was added to the reaction, and the dNTPs were removed by purification of the first strand mixture on a self-made 200 μl G-50 gel filtration spin-column equilibrated with 1mM Tris–Cl, pH 7.0. After second strand synthesis and DNA fragmentation process, the sequencing libraries were further constructed by following the manufacturer’s instructions (Illumina). Fragments of 300-400 bp were recovered and purified, and then enriched by PCR for 15 cycles. Each library was loaded into one lane of the Illumina GA IIX for 2 × 120 bp pair-end sequencing at a concentration of 2 pM, except that library of normal seedlings was loaded into two lanes. Image analysis and base calling were finished using the Illumina GA processing pipeline v1.4.

Laser microdissection (LM)-captured rice seedling leaf epidermal cells and aRNA preparation

Leaves of 15-d-old seedlings of rice variety TP309 grown in a growth chamber at 12/12h and 25/22°C of day/night cycle were used. Seedling leaves were cut into pieces about 5 mm long and immediately processed by microwave-accelerated acetone fixation (BP-111-RS laboratory microwave, Microwave Research & Applications Inc.) and paraffin-embedded as described by Tang et al. [55]. Cross-sections of 10 μm thickness, parallel to leaf vascular bundles, were obtained using a Leica RM2235 rotary microtome. Paraffin-tape transfer system (Instrumedics) and Veritas Microdissection Instrument (Acturus Bioscience) were used for LM-capturing epidermal cells.

We carried out microdissection of epidermal cells, which are a group of cells (a homogenous population cells). There were about 200-500 cells per sample. The total RNAs of epidermal cells were extracted by PicoPure RNA isolation kit (Acturus, CA, USA) with the DNase (RNase-free, Qiagen) treatment. The integrity of the total RNAs was evaluated by Agilent 2100 Bioanalyzer using RNA-6000 Pico LabChips (Agilent Technologies). Because the total RNA quantity of LM-captured epidermal cells was about 10 ng, we used a TargetAmp two-round aminoallyl-aRNA (antisense RNA) amplification kit (Epicentre Biotechnologies, Madison, WI, U.S.A.) with Super-Script III and SuperScript II reverse transcriptases (Invitrogen, Carlsbad, CA U.S.A.) to amplify laser-microdissected RNA [56]. For each amplification, approximately 0.5 ng of total RNA (in a 2-μl volume) was used as starting material and, typically, 5 to 10 μg aminoallyl cRNA was recovered. For evaluating the fidelity of two-round RNA amplification, the TargetAmp one-round aminoallyl-aRNA amplification kit (Epicentre) was used as a one-round amplification control. For each one-round amplification, 400 ng of total RNA was used as starting material and, typically, more than 10 μg of aminoallyl cRNA was recovered. The quality of amplified RNA was also evaluated by Agilent 2100 Bioanalyzer, and only those with a ‘bell-shaped’ curve with peak size > 300 nucleotides were used for RNA-seq.

Small RNA sequencing

The same materials described above were used for small RNA library construction using the Illumina small RNA preparation kit (v1.5) following the manufacturer’s instructions. Small RNAs of 18-34 bp were enriched by polyacrylamide gel electrophoresis, and ligated to 5^′ and 3^′ adapters. The ligation product was reverse-transcribed into cDNA, which was then amplified by 15 PCR cycles and subjected to Illumina’s Solexa proprietary sequencing. Here, each library was loaded into two lanes of the Illumina GA IIX for 35-bp sequencing, except that the library of cold treatment was loaded into one lane.

Northern analysis and qRT-PCR

The PCR primers and Northern probes were designed for cis- NAT validation (Additional file 14). Northern blot analysis was carried out as described [57]. To ensure that each pair of primers had specificity of its transcript, one was selected from the overlapping region of cis-NATs, and the other from the flanking region. Total RNA was extracted using the Trizol reagent (Invitrogen) according to the manufacturer’s instructions. After treated with DNaseI (NEB), 5 μg of total RNA was used to synthesize the oligo (dT) primed first-strand cDNA using SuperScript™ II reverse transcriptase (Invitrogen). For semi-quantitative RT-PCR, rice actin expression was used as a control, cis- NAT PCR products were amplified by 36 cycles, while actin products were amplified by 26 cycles. Real-time RT-PCR was performed on the Applied Biosystems 7500 real-time PCR System. Diluted cDNA was amplified using SYBR Premix Ex Taq™ (TaKaRa). The expression levels of transcripts were normalized by endogenous eEF-1α (AK061464) transcripts. Three technical replicates were taken for each set.

Assembly of rice transcripts

All ssRNA-seq data obtained from normal and cold, salt and drought stress conditions were respectively mapped to the rice reference genome [41] using software TopHat [43]. Tolerances were set to allow at most two mismatches for paired-end reads in each alignment; and reads with multiple alignments were ignored. The corresponding outputs together with RAP-DB genome annotation data [42] were subsequently used for transcript assembly to detect known and unannotated transcripts and isoforms by another software package, Cufflinks [44, 46]. Then, we integrated these transcripts together according to another command, Cuffmerge, from the Cufflinks package. Finally, referring to the renewed gene models, we calculated FPKM (Fragments Per Kilobase of exon Models) of each transcript under normal, abiotic stressed conditions, and leaf epidermal cells as well [45]. Open reading frames (ORFs) of novel transcripts were predicted using the ‘getorf’ program of EMBOSS package [58], with the longest ORF extracted for each transcript. For functional annotation, all transcripts were searched against the PFAM database [47] using HMMER v3.0 (E-value < 0.0001) [59]. The MSU Oryza Repeat Database [60] was used to determine transposable element coordinates on the rice pseudomolecules and all transcribed regions, which were annotated using RepeatMasker [61].

Identification of cis-NATs and nat-siRNAs

The renewed gene annotation results were applied as models. Those assembled genes with uncertain transcriptional orientation were filtered out first. To profile all potential cis-NAT pairs in the rice genome, transcripts that originated from the same locus but from opposite strands and with non-redundant overlapped length > 25 nt were chosen. From each group of sense-antisense pairs extracted, we selected as representative pairs those with the longest exonic overlap. We wrote perl scripts to classify cis- NAT types: convergent, divergent or enclosed. Furthermore, we examined expression levels of all assembled genes using ssRNA-seq data. Thus, we could identify expressed transcripts from putative cis- NATs. According to the number of antisense partners in each cis- NAT pair, two types of cis-NATs were defined. One was one-to-one type, i.e. one transcript in a cis-NAT pair has only one antisense partner. Another is many-to-many cis-NATs, i.e. one transcript in a cis-NAT has more than one antisense partner.

Raw data of small RNA reads from four corresponding small RNA libraries were aligned to the rice genome [41] by SMALT v0.5.7 with default parameters [62]. Mapping scores below the threshold were not reported. After removing the 5^′- and 3^′-adapters, only those continuously and perfectly matched reads with length of 18-34 bp were extracted for further analysis. Reads with multiple alignments or with any mismatches were completely filtered out. Next, we discarded those sequences that could be unambiguously mapped to rRNA, tRNA, sn/snoRNA, mitochondria and chloroplasts. We calculated small RNA densities according to a published method [28]. Small RNAs (at least one unique read) which fully located in overlapped regions of cis-NATs were extracted as putative nat-siRNAs. The cis-NATs were identified from the set of DEGs across normal and three abiotic-stress libraries. The analysis was carried out based on the following criteria: FPKM of transcripts were used for comparison by computing fold changes (with absolute value ≥ 2) and Fisher’s exact test (p < 0.001) according to an ‘R’ statistical package named ‘DEGseq’ [63].

In addition, R² is the square of the Pearson product-moment correlation coefficient relating the regressor and the response variable. R² gives some information about the goodness-of-fit of a model. In regression, R² is a statistical measure of how well the regression line approximates the real data points. ‘R’ software was used to calculate R²[64] of cis- NATs between normal and stressed conditions.

Data release

The raw sequences were deposited in the EBI European Nucleotide Archive with accession number E-MTAB-721 (http://www.ebi.ac.uk/arrayexpress/browse.html?keywords=E-MTAB-721) and ERP001962 (http://www.ebi.ac.uk/ena/data/search?query=ERP001962). The assembled transcripts can be freely downloaded xfrom http://www.ncgr.ac.cn/scientific_data.asp.