Eukaryotic organisms have evolved ingenious mechanisms to increase protein diversity. At the genome level, diversity is amplified through various mechanisms of gene duplication, such as whole-genome (WG) duplication and segment and tandem duplications, during evolution (Freeling, 2009). At the protein level, alternative RNA splicing provides another layer of diversity, whereby transcripts of many multi-exonic genes are processed with the inclusion or exclusion of particular exons in the final mRNA molecule (Berget et al., 1977; Chow et al., 1977). It has been reported that ~95% of multi-exonic genes are alternatively spliced in humans; in contrast, the number is ~33% in rice (Pan et al., 2008; Zhang et al., 2010). Pre-mRNA splicing is regulated by a system of trans-acting RNA-binding proteins that bind to cis-acting sites on the primary transcript, followed by recruitment of spliceosomal proteins (Ser/Arg-rich proteins, SRs) to the target site. By recruiting small nuclear ribonucleoproteins (snRNPs) that form the spliceosome, SRs are well known for their role in alternative splicing (Glisovic et al., 2008; Lunde et al., 2007). Although many RNA isoforms can be generated from a single locus, the majority of loci generate only two or three isoforms in abundance.

Splicing isoforms consist of several groups corresponding to patterns of exon truncation or extension, intron retention or inclusion/exclusion of entire exons in a mature transcript (Cambell et al., 2006; Kriventseva et al., 2003). Such splice modes have been classified by several research groups. For example, Lewis et al. (2003) proposed classifying splice modes based on splice site usage and effects on the coding sequence. These authors reported that splice sites can be introduced or lost according to the mode of splice type, such as an intra-exon splice, an imperfect exon skip, alternate sites, exon inclusion, and a perfect exon skip; changes in coding region were also considered. Campbell et al. (2006) classified alternative splice forms into 9 types, including alternate acceptor/donor sites and retained or skipped exons in rice. Several web tools are also available for detecting alternative splicing (Huang et al., 2003; Koscielny et al., 2009; Zhu et al., 2014). Expression Atlas (http://www.ebi.ac.uk/gxa) is a value-added database that provides information about gene, protein and splice variant expression in different cell types, organs, and developmental stages (Petryszak et al., 2014). Indeed, splicing events are often differentially regulated across tissues and during development, suggesting that individual isoforms may have specific spatial or temporal roles. Such analyses of alternative splice forms might provide clues about which genes are under alternative splicing pressure during development or caused by various biotic or abiotic stresses.

Rice (Oryza sativa) is an important crop that supplies food more than 3 billion people. Based on the entire genome sequence and map-based sequence, the genome size of rice is approximately 389 Mb. In addition, the recent International Rice Genome Sequencing Project (IRGSP) 1.0 version annotated approximately 38,000 protein-coding sequences (Price et al., 2002; Project, International Rice Genome Sequencing 2005; Salunkhe et al., 2011), among which 2,478 non-redundant sets of rice transcription factors (TFs) (Gao et al., 2006; Perez-Rodriguez et al., 2010; Priya and Jain, 2013) involving major families such as AP2-EREBP (165), bHLH (160), NAC (144), and MADS (90) were found. As the genome size of rice is relatively small, it has been used as a model monocot plant to provide information at the genome scale for numerous biology fields, including organogenesis, development, and stress responses. The attributes associated with the small genome of rice have allowed researchers in these fields to perform proteomic analysis (Choudhary et al., 2009; Degenkolbe et al., 2013; Rabello et al., 2008), transcriptomic analysis using microarrays (Degenkolbe et al., 2009; 2013; Rabbani et al., 2003; Rensink and Buell, 2005), and genetic analysis using next-generation sequencing (NGS) (Huang et al., 2012). Moreover, to identify genes related to agronomic characters, various cultivars have recently been sequenced using NGS followed by mapping to reference sequences (Arai-Kichise et al., 2014; Huang et al., 2010; Yamamoto et al., 2010).

Over the last two decades, DNA microarray technology has been extensively applied in gene expression studies in various areas of biology, such as transcriptome changes during shoot, root, and callus development (Che et al., 2006; Jiao et al., 2009; Ma et al., 2005; Su et al., 2007; Wang et al., 2010). Microarrays have also been used to determine gene expression profiles of whole plants or various organs during stress or physiological changes (Kreps et al., 2002; Lenka et al., 2011; Minh-Thu et al., 2013; Rabbani et al., 2003; Sarkar et al., 2014; Seki et al., 2002; Shankar et al., 2014; Shinozaki and Yamaguchi-Shinozaki, 2000). In these analyses, several versions of rice chips have been used, such as the GeneChip (25 bp) from Affymetrix (www.affymetrix.com) and customized chips (60 bp) from Agilent (www.agilent.com) or NimbleGen (sequencing.roche.com). A spotting technology-based NSF45K array (50–70 bp) has also been employed to distinguish alternatively spliced transcripts (Jung et al., 2008). Recently, NGS without prior annotation information was utilized to analyze organs under various physiological conditions (Kyndt et al., 2012; Wakasa et al., 2014; Yang et al., 2015a; 2015b; Zhai et al., 2013). However, due to the lack of issues associated with probe redundancy and annotation, RNA sequencing (RNA-Seq) simplifies data interpretation and is rapidly replacing microarray technology.

In this study, we designed feature probes for an alternatively spliced transcript detection microarray (ASDM) for rice. The 40,139 transcripts/36,176 loci in ASDM can be detected in the current Rice Annotation Project (RAP) database containing 44,552 transcripts/37,868 loci. ASDM was applied to differentiate transcriptomes among the leaf, root, 1-cm panicle and young seed. In a parallel analysis, the results were compared with those obtained by RNA-Seq. Both technologies revealed the significance of alternative splicing in terms of organ enrichment. The several benefits of alternatively splicing with organ bias are discussed.

Plant materials and RNA isolation

All four organs in this study were prepared from Oryza sativa L. cv. Ilmi. Leaves and roots were harvested separately at 14 days after sowing in soil. Flowers before pollination were collected in the range of 0.5–1.5-cm panicles (P1cm). Seeds at 20–30 days after pollination (S21DAP) were collected at the internal ripening stage of the grains. Total RNA was isolated from tissues using Hybrid-R™ (GeneAll, Korea) extraction following the manufacturer’s manual. For experimental reproducibility, all samples were prepared twice from distinct plants and immediately frozen in liquid nitrogen.

Microarray design

An ASDM for rice was designed based on the genome information of IRGSP_1_0_representative_2014_06_25 (http://rapdb.lab.nig.ac.jp). Four 60-nt-long feature probes as representative features (normally ending with ‘-01’) were designed for each gene, starting 60 bp upstream of the end of the stop codon and shifting 30 bp, such that 4 probes covered 60 bp of the coding sequence (CDS) and 90 bp of the 3′ untranslated region (UTR). For loci with alternatively spliced transcripts, additional probes were designed to distinguish alternatively spliced transcripts. Unique exons (UniqExon) or regions (UniqExonU or UniqExonD) were identified by comparing genome coordinates between alternatively splice transcripts. Regions created by joining adjacent exons were designated as ‘UniqExonJ’. Genes from chloroplasts (123) and mitochondria (74) and selection markers such as gfp, gus, hyg, bar, and kan were included. In total, 175,192 probes were designed, with an average probe size of 60 nt and a Tm between 75 and 85°C. The genome version contains 37,868 loci and 44,552 transcripts, and ASDM covers 36,176 loci and 40,139 transcripts. The microarray was manufactured at Agilent (Agilent Technologies, United States, http://www.agilent.com).

Microarray analysis

Probes were labeled using Agilent Low RNA Input Linear Amplification Kit PLUS following the manufacturer’s protocol (http://www.agilent.com). Total RNA (500 ng) was used to generate labeled cRNA by incorporating cyanine CTP. The microarray was scanned using an Agilent DNA microarray scanner (G2505C), and raw intensity data were extracted with the feature extraction software. The data were processed with the limma package in R (http://www.bioconductor.org/), and background adjustment and normalization were performed with the background Correct and normalize Between Arrays functions, respectively, in the program. plot Densities functions were used to analyze the distribution of probe intensities for 8 microarray data. Frequency densities before and after normalization were plotted against log2-based intensities (Supplementary Figs. S1a and S1b). The microarray intensities ranged from a minimum 23 to a maximum 545,883. Consistency between the microarrays was tested using leaves (Supplementary Fig. S1c). Log2-based intensities of 2 leaf samples were compared. The linear model was Y = 0.998x + 0.018, and the Pearson correlation coefficient was 0.996. Multiple analysis was performed after filtering out control probes and low-expression probes (Ritchie et al., 2015). The package adopted the linear modeling approach implemented by lmFit and the empirical Bayesian statistics implemented by eBayes. Genes with adjusted-P-values or false discovery rates (FDRs) below 0.05 were collected, and those with log2 scale values higher than 1 or less than −1 in at least one organ compared to that of the leaf were further selected. The data were processed with the limma package based on R statistical language.

Gene Ontology (GO) analysis

Biological term enrichment among CVI, CVII, CVIII-i, and CVIII-ii genes was assessed using GoMiner (Ashburner et al., 2000; Zeeberg et al., 2003), as described previously (Minh-Thu et al., 2013). To identify a tentative ortholog of a rice gene in the Arabidopsis genome, BLASTP analysis was performed for the two species; genes with a score of 70 or higher were considered tentative counterparts. The microarray contains 27,448 genes, and 21,556 of these match Arabidopsis genes with scores equal to or greater than 70; these genes were used as the total gene set in GoMiner. Among the 21,556 rice genes in this study, significant changes in biological processes, molecular functions and cellular components were found for 1,342, 406 and 308, respectively. As the relationships of GO terms are hierarchical, the top terms (66) in each clade for the 1,342 biological process terms were collected and subjected to hierarchical clustering using limma language in the Bioconductor project (http://master.bioconductor.org/). FDR values were obtained from 100 randomizations. GO terms for which the FDR was less than 0.05 in at least at one group were collected. For presenting GO terms in hierarchical clusters, unique GO terms were first selected. The FDR values ranging from less than 0.05 were scaled 0.5 to 5 such that the more enriched terms have values close to 5 and darker color codes in the heatmap function (http://master.bioconductor.org/). Other non-significant FDR values were transformed to 0 using Perl scripting language. To reduce redundancy, the top terms in each clad were selected, and hierarchical clustering was performed with the Bioconductor R package (http://master.bioconductor.org/).

RNA sequencing

mRNA was purified using the TruSeq RNA sample preparation kit (http://www.illumina.com/) following the company’s protocol. Following purification, the mRNA was fragmented, and first-strand cDNA was generated using reverse transcriptase and random primers. Second-strand cDNA synthesis was then performed using DNA polymerase I and RNase H. These cDNA fragments were subjected to an end-repair process, the addition of a single A base, and adaptor ligation. The products were purified and enriched by polymerase chain reaction (PCR) to create a cDNA library in a bridged amplification reaction that occurred on the surface of the flow cell. A flow cell containing millions of unique clusters was loaded into a HiSeq 2000 for automated cycles of extension and imaging. Each sequencing cycle was carried out in the presence of all four nucleotides, generating a series of images with each representing a single base extension at a specific cluster. Libraries were constructed for the leaf, root, P1cm, and S21DAP. Repetitions were performed for all samples, except roots. In the initial sequencing, 2.7–0 Gb was sequenced at each end. These results generated 16–43 × 10⁶ paired-end reads.

Data processing

Programs such as TopHat, Cufflinks and HTSeq were used to process sequence reads. Differential exon usage was examined with the DEXseq package in Bioconductor. Briefly, the rice genome sequence (IRGSP-1.0_2014-06-25) in the RAP database (http://rapdb.lab.nig.ac.jp) was indexed using Bowtie 2 with an FM Index based on the Burrows-Wheeler transform (Langmead et al., 2009). Reads in FASTQ files were aligned using TopHat (v2.0.4) with default parameters (Kim et al., 2013). TopHat aligned the RNA-Seq reads to genomes using the short read aligner Bowtie 2 and then analyzed the mapping results to identify splice junctions between exons. In the program, samtools (version 1.2) was used for indexing and sorting of bam and sam files. For transcriptome assembly and isoform quantitation from the RNA-Seq reads, Cufflinks was used and tested for comparison (Trapnell et al., 2012). Prior to statistical analysis of exons using the DEXseq package, per-gene read counts were performed with htseq-count in HTSeq-0.6.1 using the annotated gtf file (Anders et al., 2015). To extract exons among transcripts, non-overlapping exonic regions were defined with dexseq_prepare_annotation.py, and then the per-exon read counts were extracted with the script dexseq_count.py in the DEXSeq package (Anders et al., 2012). For each exon of each gene, the data contain the number of reads in each sample that overlap with the exon. If an exon’s boundary is not the same in all transcripts, the exon is separated into two or more parts, referred to as “counting bin”. A read that overlaps with several counting bins of the same gene is counted for each of these. The program returned one file for each biological replicate with the exon counts. Finally, differential exon usage was assessed with the DEXseq package in Bioconductor (http://bioconductor.org). In the package, generalized linear models (GLMs) for model read counts and parameters are calculated from a negative binomial (NB) distribution of these count data (McCullagh and Nelder, 1989).

Reverse transcription PCR (RT-PCR) and real-time PCR

For RT-PCR, 5 μg of total RNA was reverse transcribed using oligo (dT) and a RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Scientific™, USA). PCR was performed with one-hundredth of the first-strand cDNA mixture and gene-specific primers. The PCR conditions used were as follows: 10 min at 94°C and 30 cycles of 30 s at 94°C, 30 s at 55°C, and 30 s at 72°C. The 20 μl contained 10 μl Solg™ 2X Real-Time PCR Smart mix [Solgent, Korea], 2 μl cDNA, 1 μl primers, and 7 μl water, and real-time PCR was performed using the CFX96 touch real-time PCR detection system (Bio-Rad, USA) with three technical replicates. Expression was assessed by evaluating threshold cycle (CT) values. The relative expression level of tested genes was normalized to the tubulin gene and calculated using the 2^−ΔΔCT method (Livak et al., 2001). The sequences of all primers used in this study can be found in Supplementary Table S11.

Target probes were designed based on unique sequences among alternatively spliced transcripts

It has been reported that the 3′ regions of genes in plants are more specific than the 5′ regions (Eveland et al., 2008). To reflect these features, four 60-nt-long feature probes were designed for each gene, starting 60 bp upstream of the end of the stop codon and shifting 30 bp, such that the four probes covered 60 bp of the CDS and 90 bp of the 3′ UTR (arrows collectively indicate R for locus Os08g0482700 in Fig. 1A). To distinguish alternatively spliced transcripts, additional probes were designed based on unique sequence variations among alternatively spliced transcripts. For example, two transcripts, Os08t0482700-01 (encoding cupredoxin domain-containing protein) and Os08t0482700-02 (encoding conserved hypothetical protein), are alternatively spliced from mRNA produced by the Os08g0482700 locus. In addition, probes were designed to represent the unique features of each transcript. The Os08t0482700-01 transcript is composed of two exons and a 546-bp CDS contained within these exons. In contrast, the Os08t0482700-02 transcript is composed of one exon and has a 228-bp CDS located within this exon. The Os08t0482700-01 transcript has unique junction sites that join the spliced transcripts, and a 60-bp feature probe (UniqExonJ) was designed based on two 30-bp-long regions in the CDS for each exon. Thus, the probe represents a part of the unique exon (Os08t0482700-01_UE) of Os08t0482700-01 and is used to represent Os08t 0482700-01 alone. When representative probes cannot distinguish between alternatively spliced transcripts, ‘-01_UE’ can be used for expression of the ‘-01’ transcript. The Os08 t0482700-02 transcript was compared to Os08t0482700-01; a portion of the 5′ UTR, the CDS and a portion of the 3′ UTR are unique, and these regions were designated as UniqExonU, UniqExon, and UniqExonD, respectively. Three 60-bp target probes were designed according to these regions to represent Os08t0482700-02.

A total of 37,868 genes are representatively annotated in the IRGSP and RAP databases (http://rapdb.lab.nig.ac.jp), and 5,254 loci encode 11,938 alternatively spliced transcripts (Table 1). Among these alternatively spliced transcripts, 7,502 transcripts were designed with the unique exon features indicated above. Ultimately, 175,192 feature probes representing 40,139 of the 44,552 transcripts were designed and designated ASDM.

Confirmation of organ preferential genes with alternative splicing transcripts

We applied the rice ASDM to detect alternatively spliced genes among the leaf, root, panicle at 0.5–1.5 cm (P1cm) and seed at 20–30 days after pollination (S21DAP). The 175,193 feature probes for 41,953 transcripts were normalized as described in the Methods section (Supplementary Figs. S1a and S1b). Consistency for intensities between microarrays using total RNA from leaves was higher than 0.99, suggesting that the experiment was very reproducible (Supplementary Fig. S1c). The data have been deposited at National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO); GSE87536. Normalized intensities are presented in Supplementary Table S1.

To confirm the detection of alternatively spliced transcripts by ASDM, several genes were tested. As shown in Fig. 1B (ASDM), among the four organs, the strongest intensity for Os08t0482700-01 were found in the root. The Os08t04827 00-01_UE transcript represented by the unique junction Os08t0482700-01 confirmed this result. The intensities of Os08t0482700-01 and Os08t0482700-01_UE showed a similar pattern among the organs, being highly expressed in the root, followed by the leaf, S21DAP and P1cm. Although Os08t0482700-02 exhibited lower expression compared to Os08t0482700-01, similar to Os08t0482700-01 and Os08t0 482700-01_UE, it also was highly expressed in the root. Semi-quantitative PCR confirmed these patterns of expression (Fig. 1C). Primers were designed to represent the unique regions of each transcript, and as expected, these primers amplified segments presenting a tendency similar to the microarray data shown in Fig. 1B.

We also examined Os04g0612500, Os10g0471100 and Os07g0673400 in the leaf, P1cm and S21DAP samples, respectively (Supplementary Fig. S2). The probe intensities of the microarray showed the highest expression of Os04g06 12500 in the leaf, which was confirmed by the result for Os04t0612500-01_UE (encoding a 36.4-kDa proline-rich protein) (Supplementary Fig. S2A). The probe designed using the 5′ and 3′ ends of Os04t0612500-02 also revealed expression in the leaf. RT-PCR using primers for Os04t 0612500-01, Os04t0612500-01_UE and Os04t0612500-02 amplified products of the expected sizes. Through evaluation of mRNA levels by real-time PCR, approximately 300 times higher expression in the leaf compared to the other tissues was found. Real-time PCR results also indicated strong expression in the leaf.

The microarray data showed that transcript Os10t04711 00-01 (encoding a protein with high similarity to proteins involved in cuticular wax biosynthesis) is specifically expressed in the P1cm stage (Supplementary Fig. S2B). The intensity of Os10t0471100-01_UE also confirmed that this transcript is highly expressed at this stage. The degree of intensity of Os10t0471100-02 (CER1) was lower than that of Os10t0471100-01, and this transcript was also found to be exclusively expressed in the P1cm sample. In addition, the mRNA levels detected by real-time PCR showed enhanced expression of this transcript in P1cm.

The microarray probe intensities showed the highest expression of Os07g0673400 in S21DAP, and the result for Os07t0673400-01 (encoding Rossmann-like alpha/beta/alpha sandwich fold domain) confirmed this result (Supplementary Fig. S2C). The probe designed using the 5′ end of Os07t0673400-02 also revealed expression in the S21DAP sample. RT-PCR with primers representing Os07t0673400-01, Os07t0673400-01_UE and Os07t0673400-02 generated products of the expected sizes. The mRNA levels detected by real-time PCR showed heightened expression of this transcript at the S21DAP stage. These data suggest that ASDM can detect expression of loci in the tested organs as well as alternatively spliced transcripts.

Genome-wide comparison of rice ASDM and RNA-Seq suggest co-linearity between the two technologies

In a parallel analysis, RNA-Seq was performed using the Illumina HiSeq 2000 platform, as indicated in the Methods section. Approximately 17–43 million reads were obtained from 30–80 Gb raw reads (Supplementary Table S2). These data were aligned to the reference genome 15–22 million times with TopHat, a splice-aware aligner (Kim et al., 2013). Unmapped hits varied widely; the lowest number of reads (1.1 M reads) was found in the seed and the highest number (17.6 M) in one of the leaf samples. The distribution of raw RNA-Seq data is shown in Supplementary Fig. S3. The dataset was loaded and normalized with the estimate Size Factors function in the DEXSeq package (Anders et al., 2012). In Supplementary Fig. S3, the exon dispersion estimates, fitted mean-dispersion values function and shrunken values are indicated by black points, red lines and blue points, respectively. Finally, median values for transcripts were obtained. The leaf sample provided the number of transcripts aligning more than once in the genome and showed the highest number of transcripts (27,852), with the panicle at 1 cm sample displaying the lowest (25,507) (Supplementary Table S3). The transcripts with maximum counts for the leaf, root, P1cm, and S21DAP were Os11t0707000-01 (119,487.5, encoding ribulose-bisphosphate carboxylase activase), Os12t0571100-01 (54,236, encoding metallothionein-like protein type 1), Os02t0121300-01 (34,784, encoding cyclophilin 2) and Os01t0762500-00 (224,247, encoding glutelin subunit), respectively. The intensities of these transcripts by ASDM were also very high, at 313,918 (Os11t0707000-01), 416,636.6 (Os12t0571100-01), 244,628.4 (Os02t0121300-01), and 251,958.8 (Os01t0762500-00), respectively, suggesting co-linearity between the two technologies. In terms of transcripts with low numbers, among the 40,139 transcripts in the microarray, 12,287–14,632 transcripts were counted less than once by RNA-Seq. The median value of these transcripts by ASDM was 100.5.

We tested genome-wide co-linearity between the microarray intensities and RNA-Seq counts. The log2-based average intensities of the leaf microarray data and the log2-based counts for the organs were compared (Fig. 2). Pearson’s correlation from two methods ranged 0.30–0.36 for the leaf, root, P1cm and S21DAP, suggesting correlation between the data. We also downloaded publicly available RNA-Seq data from TENOR (http://tenor.dna.affrc.go.jp), and data for the leaf and root at 10 days after germination were compared to those in ASDM. Pearson correlation coefficients ranged 0.25–0.27. The p-values of the F-statistics were also less than 2.2e-16 and significant (Supplementary Fig. S4). The distributions were similar to those for the leaf and root, as described above.

A list of exons of the transcripts was generated from a GTF file in IRGSP-1.0 subversion 2016-03-09 with a python script, as described in the Methods section (Anders et al., 2012). The overlapping exon regions are marked as “counting bins” by the program, and the read counts are distinguished by the statistical model in the program. The bins of transcripts such as Os08t0482700-01 and Os08t0482700-02 for the root (Fig. 1B, RNA-Seq), Os04t0612500-1 and Os04t0612500-2 for the leaf Os10t0471100-1 and Os10t0471100-2 for P1cm and Os07t0673400-1 and Os07t0673400-2 for S21DAP (Supplementary Fig. S2, RNA-Seq in b panel) are listed in Supplementary Table S4. The pattern of the exon counts for these transcripts was compared to the ASDM data, with a very similar result for each organ, even though expression values can vary depending on the technology employed.

Application of rice ASDM and RNA-Seq technologies to distinguish between genes that are enriched in organs and those that are constitutively expressed across organs

Among 45,128 transcripts, the number with normalized counts below 1 ranged from 12,000 to 15,000, and the median intensity of ASDM for these transcripts was 100.5. Based on the microarray data, the transcripts (32,718) with intensity greater than 300 (3* the median intensity of the transcripts, where RNA-Seq counts were greater than 1) in at least one organ were chosen for further analysis. The average values are presented in Supplementary Table S5. Among these transcripts, 21,497 were for a single exon, and 11,129 were for loci with alternatively spliced transcripts. To identify genes differentially expressed in each organ, we compared normalized intensities between distinct organs by adapting a statistical analysis approach using values of coefficient of variation (CV), which is the standard deviation of the mean and is used as a standardized measure of dispersion for probability distribution or frequency distribution (Everitt, 1998). Thus, a higher CV value suggests high variation among samples and is likely to be organ preferential in our data. This technique was also adopted by a rice gene expression database, CREP (Collections of Rice Expression Profiling, http://crep.ncpgr.cn, Wang et al., 2010), and presented in our previous paper (Chae et al., 2016). CREP contains 190 Affymetrix GeneChip Rice Genome Arrays using RNA from 39 tissues collected throughout the life cycle of two indica varieties: Zhenshan 97 and Minghui 63. To delineate CV values in our dataset, we attempted to compare our results with those in CREP. We also estimated the distribution of the 100 genes most constitutively expressed in CREP in our ASDM and RNA-Seq datasets, as described previously (Chae et al., 2016). The mean values of these genes were 0.18 (sd 0.02), 0.11 (sd 0.05) and 0.34 (sd 0.19) for CREP, ASDM and RNA-Seq, respectively (Supplementary Fig. S5).

Based on these CV values in ASDM, genes were tentatively categorized into 3 groups, as shown in Fig. 3A: CVI (higher than 1.1), CVII (higher than 0.7), and CVIII (below 0.7). We classified these groups as an organ-preferential group (CVI, 5,410), an organ-enriched group (CVII, 6,547), and a nonspecific group (CVIII, 20,761). CVIII was sub-divided as -i (highest), -ii (medium) and -iii (lowest) according to the degree of expression according to ASDM. As the CV values of RNA-Seq were comparable to those of ASDM, the CV values from RNA-Seq are also categorized as those from ASDM in Fig. 3B and Supplementary Table S3: CVI (higher than 1.4), CVII (higher than 0.8), and CVIII (below 0.8). The number of transcripts for the CVI, CVII and CVIII categories were 5,721, 8,574 and 18,371, respectively. CVIII was also sub-divided as -i (highest), -ii (medium) and -iii (lowest), at 3,467, 5,858, and 9,046, respectively.

When the intensities and counts were set ≥ 300 and ≥ 1 for ASDM and RNA-Seq, respectively, the effective number of expressed genes was approximately 37,000 for both technologies. Genes commonly detected using these technologies were searched according to CV group (Fig. 3C): approximately 40–50% of transcripts were commonly detected for each of the CVI and CVII groups, whereas approximately 60% of transcripts were detected for CVIII. For each technology, the number of transcripts commonly detected for most CV groups was greater than the number of excluded transcripts.

Hierarchical clustering of genes with a single transcript and alternatively spliced transcripts

As the genes of the CVI group for ASDM showed organ enrichment, those with a single transcript (4,172) and alternatively spliced transcripts (1,169) in this group were delineated by hierarchical clustering (Fig. 4). Highly expressed genes (618) in the leaf with a single transcript were enriched for a metabolic process represented by photosynthesis (Supplementary Table S6). These genes were Os07t0556200-01 (encoding Rieske [2Fe-2S] region domain-containing protein), Os02t0578400-01 (encoding photosystem II oxygen-evolving complex protein, PsbQ family protein), and Os07t0105600-01 (encoding photosystem II oxygen-evolving complex protein, PsbQ family protein). At least 6 C₂C₂-CO-like and putative zinc finger TFs that promote flowering induction in Arabidopsis under long photoperiods, such as Os06t0264200-00, Os09t0240200-01 and Os03t0711100-01, were highly expressed in the leaf. In the root, 665 genes were highly expressed and clustered into 3-4 groups. The root system development-related genes Os07t0128800-00 (encoding PR1a protein), Os05t0135200-01 (encoding plant peroxidase domain-containing protein), and Os02t0662000-01 (encoding Rcc3 protein) were highly expressed. Five AP2-EREBP family members, including Os03t0341000-01, Os01t 0752500-00 and Os03t0191900-01, were also strongly expressed. It is worth noting that approximately 30 non-protein coding transcripts, such as Os02t0791650-00, Os05t0400301-00, and Os02t0178901-00, were similarly highly expressed. For P1cm, the reproductive process was strongly activated, with Os10t0489900-01 (encoding CDT1a protein), Os03t0800200-00 (encoding protein argonaute MEL1), Os03t0137800-01 (encoding cyclin-dependent kinase inhibitor family protein) being highly expressed, as were the MADS box TFs Os05t0423400-01 (encoding PISTILLATA-like MADS box protein), Os02t0682200-01 (OsMADS6), and Os07t0605200-01 (OsMADS18). Highly expressed genes in S21DAP clustered into 3 groups. Processes related to seed storage metabolites such as amylopectin, amino sugar, and lipid localization-related processes were enriched, including Os08t0520900-00 (encoding iso amylase), Os02t0456100-01 (encoding glutelin), and Os03t 0808500-00 (encoding plant lipid transfer protein/Par allergen family protein), and starch biosynthesis-related genes such as Os06t0229800-01 (encoding starch synthase IIA.) and Os06t0160700-01 (encoding starch synthase I chloroplast precursor, EC 2.4.1.21). CCAAT-containing TFs Os10t 0191900-00, Os10t0397900-00, and Os01t0102400-01 and NAC domain TFs such as Os03t0327800-01, Os01t03931 00-01 and Os01t0104500-01 were also highlighted.

Highly expressed genes with alternatively spliced transcripts in the leaf clustered into 4-5 groups (Supplementary Table S7). Metabolic process, cellular component organization and biogenesis were enriched by the genes Os12t 0277500-01 (encoding RuBisCO subunit binding-protein, alpha subunit chloroplast precursor) and Os12t0420200-01 (encoding NAD(P)-binding domain-containing protein). In the leaf, alternatively spliced transcripts for TFs such as Os04t0497700-01 (encoding C₂C₂-CO-like, CONSTANS-like protein) and Os06t0348800-01 (encoding G2-like) were found. Two to three groups were strongly expressed in the root, coding for proteins that are induced under decreased oxygen levels, such as Os11t0210500-01 (encoding alcohol dehydrogenase) and Os03t0183000-01 (encoding AP2 domain-containing protein RAP2.6). Transporters such as Os10t0444700-02 (encoding phosphate transporter 6), Os02t0650300-01 (encoding iron-phytosiderophore transporter for iron homeostasis) were also highly expressed. With regard to P1cm, the reproductive process and developmental process were enriched, with Os08t0538700-01 (encoding retinoblastoma-related protein) and Os08t0442700-01 (encoding similar to SERK1) highly expressed. Os07t 0108900-02 (encoding MADS domain-containing protein for sexual reproduction) was highly expressed compared to Os07t0108900-01; Os02t0258300-01 (encoding zinc finger FYVE/PHD-type domain-containing protein) was also highly expressed. Four C₃H family members, such as Os03t0329200-01 and Os03t0328900-01, were also highly expressed. For S21DAP, genes involved in the generation of precursor metabolites and energy, such as Os01t0663400-01 (encoding aspartic proteinase oryzasin 1 precursor) and Os08t0345800-01 (encoding glucose-1-phosphate adenylyltransferase), were enriched. Two prolamins, Os05t0329100-01 (encoding prolamin) and Os05t0329300-01, were also highly expressed.

To examine biological term enrichment, GO analyses were performed for genes of the CVI, CVII, CVIII-i, and CVIII-ii groups from the ASDM dataset (Ashburner et al., 2000; Zeeberg et al., 2003). This comparison was performed for loci of whole genes and those with single (16,999) or alternatively spliced (5,844) transcripts (Supplementary Figs. S6a and S6b). GO analysis of loci of whole genes by RNA-Seq was also performed (Supplementary Fig. S6c). These data suggested that GO analysis of the genes with single transcripts and alternatively transcripts reveal the representative characteristics of organs, even though subtle differences between the ASDM and RNA-Seq datasets were observed.

The portion of genes with alternatively spliced transcripts is increased among genes that are constitutively expressed across organs

In the rice genome, there are 5,254 loci with more than two alternative transcripts (http://rapdb.lab.nig.ac.jp); among these, 3,539 loci were included in the current microarray. We determined the influence of organ specificity on the number of alternatively spliced transcripts or the influence of the degree of expression on the genome. As shown in Table 2, the ratio of loci with alternatively spliced transcripts was 0.14 (5,254/37,869) in the current genome annotation and 0.15 (5,244/36,141) in ASDM; these ratios were designated as the index of loci with alternatively spliced transcripts (IAST) in the group. IAST of expressed genes was 0.19, slightly higher than those of the genome and ASDM. These values were 0.14, 0.22 and 0.24 for CVI, CVII, and CVIII, respectively. The CVI group was subdivided into the leaf, root, P1cm, and S21DAP, with values of 0.26, 0.08, 0.13, and 0.09, respectively. For the CVII group, these values were 0.18–0.28 and generally higher than those in the CVI group. When CVIII was divided into -i (highest), -ii (medium) and -iii (lowest) according to the degree of expression, IAST values of these groups ranged from 0.31–0.34. The IAST values appeared to increase from CVI to CVIII, suggesting that the ratio of genes with alternative splicing increased as the genes were expressed throughout each organ. We tested this hypothesis by performing a literature review using the web-based database CREP (Wang et al., 2010) and projected the IAST values of the 100 genes that were most constitutively expressed in CREP. Interestingly, IAST was 0.43 (Supplementary Table S8 and Consti100 in Table 2) for these constitutively expressed genes, higher than the values reported above.

We also assessed IAST for CV groups based on RNA-Seq analysis (Table 3). The IAST values for CVI, CVII, and CVIII were 0.09, 0.12 and 0.21, respectively, and comparable to those based by ASDM. The values among the CVI and CVII organs were relatively higher in the leaf and lower in the root, P1cm, and S21DAP samples, but the values were somewhat lower than those based on ASDM. The IAST values among the CVIII subgroups were higher than those among other groups, as shown for ASDM, and were also less than those for 100 constitutive genes (Consti100 in Table 2). These data confirm the observation in the ASDM analysis, whereby IAST of genes constitutively expressed throughout organs was higher than that of organ-enriched genes.

A transcript confers organ preference in a locus, and the alternatively spliced transcript(s) confers the same preference or reduced organ specificity

We tested organ preference among alternatively spliced transcripts of the CVI group enriched in each organ. The number of these transcript variants were 527, 179, 120 and 124 from 367, 137, 91 and 96 loci of the leaf, root, P1cm, S21DAP, respectively (Table 4). We tested CV groups and organ preference between a major transcript and its alternatively spliced transcripts. For example, two transcripts, Os04t0612500-01 and Os04t0612500-02, are alternatively spliced from pre-mRNA of Os04g0612500 (Supplementary Fig. S2A). Os04t0612500-01 is classified as a CVI_Leaf and Os04t0612500-02 is also preferential to the leaf thus these transcripts showed the same bias for the leaf; the gene can be grouped as loci_CVI_same_organ. Os10t0471100-01 and Os10t0471100-02 from Os10g0471100 is an example of a grouping for P1cm (Supplementary Fig. S2B). Similarly, the lists of transcript variants that are classified in the same organ were 305, 81, 57 and 54 from 147, 39, 28 and 26 loci for the leaf, root, P1cm, S21DAP, respectively (Supplementary Table S9). In contrast, a number of transcript variants were found for which the major transcript belongs to CVI but the alternatively spliced transcript(s) belongs to as_ts_Not_CVI; CVII, CVIII or a group expressed in low amount. For example, Os07t0673400-01 from Os07g0673400 (Supplementary Fig. S2C) belongs to the seed-preferential group, and its alternatively spliced transcript, Os07t0673400-02 belongs to CVII and is classified as_ts_Not_CVI. Os08t0482700-02 from Os08g0482700 is expressed in a low amount and may be an additional example of as_ts_Not_CVI (Fig. 1 and Supplementary Table S4). The numbers of as_ts_Not_CVI were 222, 98, 63 and 70 for the leaf, root, P1cm and S21DAP. Interestingly, none of the alternatively spliced transcript variants classified as CVI showed other organ specificity or bias, and these might be grouped as as_ts_CVI_diff_organ (Table 4). These finding suggest that any two alternatively spliced transcripts in the CVI group tended to show the same organ preference or the other alternatively spliced transcript showed less organ preference.

A longer CDS transcript among alternatively spliced transcripts at a locus tends to be more expressed

We tested the extent of expression among spliced transcripts within a locus by comparing transcripts with the highest and lowest expression values. To simplify expression variability, we compared the log2-based ratio between transcripts with the highest and lowest expression according to the ASDM dataset. The log ratio of the highest and lowest values is shown for all tested organs in Supplementary Fig. S7. For example, Os12t0291100-01 exhibited 14.6-fold higher expression than that of the lowest, Os12t0291100-02 (468,018.3/32,040.6), in the leaf at locus Os12t0291100 (encoding ribulose bisphosphate carboxylase small chain). The log2-based values are presented according to organ for the CVI group (Supplementary Fig. S7). For instance, CVI_Leaf showed 323 loci with alternatively spliced sites. The mean log2-based value for transcripts with the highest and lowest expression was 3.3 (sd = 1.9). The other organs exhibited similar mean ranges. We also expanded the evaluation to the CVIII group (Supplementary Fig. S8) using 3,382 of 4,028 loci to determine values. The mean value was 2.7, similar to that of the CVI organ-preferential group, suggesting that reciprocal expression was complementary among transcripts within a locus, regardless of organ specificity. We confirmed similar variation in tests based on RNA-Seq analysis (data not shown).

Campbell et al. (2006) classified types of alternative splicing into 9 groups based on the results of Program to Assemble Spliced Alignments (PASA): alternate acceptor, alternate donor, alternate terminal exon, retained exon and skipped exon, initiation within an intron, termination within an intron, and spliced intron and retained intron. We categorized the type of alternative splicing to assess splice types and variation in expression between transcripts with the highest and lowest expression. PASA was applied first to whole transcripts, detecting 5,028 alternative splices for 2,438 of 4,143 loci in the genome and 3,598 alternative splices for 2,415 loci in the microarray (Supplementary Fig. S9 and Supplementary Table S10).

It has been reported that expression varies among alternatively spliced transcripts. Nonsense-mediated mRNA decay (NMD) has been postulated as the target of RNA degradation (Lewis et al., 2003). When an alternative splice is introduced, the stop codon is classified, and the corresponding mRNA isoform is labeled an NMD candidate, resulting in lower expression compared to the isoform with a longer CDS. If the spliced transcript escapes NMD, the result is a transcript with a shorter CDS. To examine how CDS length, due to alternative splicing, affects the degree of expression, transcripts with alternative splicing are sorted by the degree of expression and the CDS lengths compared according to organ and CV group. As shown in Table 5, among alternatively spliced transcripts, longer CDS transcripts appeared to be more expressed than shorter ones for all organs and CV groups. A paired t-test showed a p-value of 0.029, suggesting that these differences are significant.

Genome-wide gene expression profiling using microarray and RNA-Seq technologies has provided abundant information regarding the basic biological mechanisms underlying development, stress responses, and other processes (Gordon et al., 2014; Ma et al., 2005; Severin et al., 2010; Shinozaki and Yamaguchi-Shinozaki, 2000; Wakasa et al., 2014). In the current version of rice ASDM, 60-bp-long oligomers were used, and 40,953 of the 44,552 transcripts deposited in IRGSP 1.0 were included. We improved the feature probe design to assess alternative splicing. Because long oligomers can be designed for spliced regions composed of two exons, the approach may have advantages over short oligomer-based microarrays. Several loci, such as Os08g0482700, Os04g0612500, Os10g0471100, and Os07g0673400, with alternatively spliced transcripts were discerned in four organs: the leaf, root, early stage of the panicle at 1 cm and seed at 20–30 days after pollination (Fig. 1 and Supplementary Fig. S2). Genome-wide comparison with RNA-Seq data suggested co-linearity of transcript expression between rice ASDM and RNA-Seq, even though the counts for RNA-Seq were somewhat compressed (Fig. 2 and Supplementary Fig. S3), which might be due to saturation problems associated with hybridization-based microarrays. The results for distribution showed a wide range of variation for genes with low expression values according to both technologies. Lower reproducibility has been reported for gene classes that are less strongly expressed, even for RNA-Seq data, despite the excellent overall correlation (Mortazavi et al., 2008). Several factors, such as probe design for ASDM and PCR amplification for RNA-Seq, can affect the number of genes determined to be expressed in low amounts, which appears to be the case with our ASDM and RNA-Seq datasets. Nonetheless, these technologies do exhibit colinearity for highly expressed genes. In a study comparing RNA-Seq with a model organism, Saccharomyces cerevisiae, high consistency between the two platforms for gene expression-level analysis (correlation ≥ 0.91) was reported (Nookaew et al., 2012), and in a study of human T cell activation, Zhao et al. (2014) reported that RNA-Seq performed better in detecting low-abundance transcripts and demonstrated a broader dynamic range than microarray. Although analysis of RNA-Seq data requires extensive bioinformatic skills and computer resources (CPU and RAM), with improvement in cost and analytical complexity, RNA-Seq is being utilized more frequently.

ASDM and RNA-Seq provide similar GO profiles for the transcripts, as based on CV values

CV values for constitutive expression vary according to the number of organ samples. In this report, CV values of 0.7 and 0.8 were tentatively used for ASDM and RNA-Seq, respectively (Fig. 3). Wang et al. (2010) and our previous study (Chae et al., 2016) reported values of 0.1 and 0.3 using 39 tissues and 7 representative organs/tissues, respectively. For both technologies, the variation in and criteria for constitutive expression increased as the number of samples decreased. Because values vary depending on the approach, the medium range of CVII groups was adopted as an intermediate for each technology. These demarcations showed that analysis of genes with alternatively spliced transcripts can be achieved for both in a similar manner. The effective number of expressed transcripts was approximately 33,000 for ASDM and RNA-Seq (Fig. 3C), and the genes commonly detected ranged from 40 to 70%, more than the number transcripts exclusive for each technology in each CV group. Among the organs examined, higher numbers for constitutively expressed genes in both datasets were found for commonly expressed group. In particular, the results of GO analysis of the CV group between ASDM and RNA-Seq were quite similar (Supplementary Fig. S6), suggesting that both technologies are applicable for transcriptome analysis. Hierarchical clustering of the CVI group among the organs showed representative transcripts and TFs. For example, representative processes were observed, such as photosynthesis-related terms for the leaf, metallothionein for the root, cdc2 kinase for P1cm and glutelin and starch-debranching enzyme for S21DAP (Yuan et al., 2008; Yun et al., 2011). GO analysis was conducted on the loci of entire genes, followed by GO analyses of single-transcript and alternatively spliced transcript loci. In general, GO enrichment according to loci with alternatively spliced transcripts was in agreement with those with a single transcript.

Expressed TFs with a single transcript in P1cm included MADS box TFs Os05t0423400-01 (encoding PISTILLATA-like MADS box protein), Os02t0682200-01 (OsMADS6), and Os07t0605200-01 (OsMADS18). In a previous report, prolonged expression of Os02t0682200-01 and Os07t0605200-01 to the stage after pollination was observed (Chae et al., 2016), suggesting that these genes are involved in panicle development. MADS box and AP TFs are crucial for floral morphology (Coen and Meyerowitz 1991; Pelaz et al., 2000; Yun et al., 2013). Highly expressed genes in the S21DAP sample included those coding for CCAAT-containing TFs, such as Os10t0191900-00, Os10t0397900-00, and Os01t 0102400-01, and NAC domain TFs, such as Os03t0327800-01, Os01t0393100-01 and Os01t0104500-01. NAC genes, which constitute one of the largest families of plant-specific TFs and are present in a wide range of species, have been implicated in the regulation of many developmental processes, including seed development, embryo development, and shoot apical meristem formation (Duval et al., 2002; Kim et al., 2007a; Sperotto et al., 2009). These NAC TFs are also enhanced during the seed stage according to the Hana database (http://evolver.psc.riken.jp).

Analysis of the transcriptome using ASDM might be useful for alternatively spliced transcripts that include TFs. Indeed, a TF with an alternatively spliced transcript in P1cm, Os07t 0108900-02 (encoding MADS-domain-containing protein for sexual reproduction), was highly expressed compared to Os07t0108900-01. In rice, florigen and Hd3a, associated with the 14-3-3 protein in the cytoplasm, translocate to the nucleus and regulate expression of the MADS box-containing TF gene Os07g0108900 (OsMADS15, Tsuji et al., 2013), which is a homolog of an Arabidopsis floral identity gene APETALA1. Overexpression of the gene results in precocious phenotypes with regard to early internode elongation, shoot-borne crown root development, reduced plant height and early flowering at the seedling stage (Lu et al., 2012). As both transcripts contain a MADS box, an expression vector may need to be designed to discern between them. Os02t0258300-01 (encoding zinc finger FYVE/PHD-type domain-containing protein) was also highly expressed, as were four C₃H family members, such as Os03t0329200-01 and Os03t0328900-01. In the S21DAP sample, Os12t 0621100-01, Os12t0621100-02 (encoding C₂C₂-YABBY TF) and Os03t0799600-01 (encoding PHD-containing TF) were significantly expressed from loci with alternatively spliced transcripts. Some YABBY genes are expressed in the apical-ventral region during vascular tissue development and early embryogenesis (Itoh et al., 2016; Liu et al., 2007), suggesting that these TFs are involved in early seed development.

Constitutively expressed genes are under greater alternative splicing pressure

Pre-mRNA or alternative splicing is found in eukaryotic organisms, some of which become highly developed through various processes. Many trans- and cis-acting factors have been identified, and one basic question with regard to alternative splicing pressure for common and organ-preferential genes remains. To address this question, we examined the influence of organ preference on the number of alternatively spliced transcripts and the influence of the degree of expression on the genome. Ratios were designated as IAST. The IAST value of expressed genes was 0.19, slightly higher than those of the genome and ASDM. These values were 0.14, 0.22 and 0.24 for CVI, CVII, and CVIII, respectively (Table 2). When CVIII was divided according to the extent of expression, IAST ranged from 0.31–0.34. IAST values also were higher in groups with greater constitutive expression, culminating in 0.43 for the most constitutive 100 genes in the literature (Wang et al., 2010). These observations were also confirmed by RNA-Seq data (Table 3). Although IAST values for CVI, CVII, and CVIII were 0.09, 0.12 and 0.21, respectively, slightly lower than those based on ASDM, the tendency toward increase with constitutively expressed groups and variations among organs were similar in these independent analyses.

To date, it remains unclear why constitutively expressed genes would have a higher IAST. All exon boundaries require canonical splice sites (e.g., GT-AG, GC-AG, AT-AC), and all bases in the genome can be exposed to mutation (Burset et al., 2000; Lorkovic et al., 2005; Mount and Steize 1981). DNA damage events can occur at a high frequency in higher eukaryotes, such that a single mammalian cell can accumulate ~10,000 abasic (apurinic) lesions per day (Lindahl and Nyberg, 1972). For comparison, mutation rates in microbes with DNA-based chromosomes are close to 1/300 per genome per replication (Drake et al., 1998). A high level of transcription has been associated with elevated spontaneous mutations and recombination rates in eukaryotic organisms (Datta and Jinks-Robertson, 1995). For example, the effect of transcription on the rate of spontaneous mutation in the yeast S. cerevisiae was examined using a lys2 frameshift allele under the control of a highly inducible promoter. As the rate of spontaneous reversion was shown to increase when the mutant gene was highly transcribed, transcriptionally active DNA and enhanced spontaneous mutation rates are correlated in yeast. Several hypotheses have been proposed to explain transcription-associated mutagenesis, e.g., transcription impairs replication fork progression in a directional manner, and DNA lesions accumulate under high-transcription conditions (Kim et al., 2007). Genes constitutively expressed across organs might have higher rates of mutation than organ-preferential genes, which may undergo reduced transcription. However, a high rate transcription might not be the sole factor explaining this phenomenon. Subdivision of CVIII by intensity, i.e., CVIII-i, -ii, and -iii, did not appear to result in variability (Tables 2 and 3). Alternatively, the neutral theory of molecular evolution may offer an explanation, such that organ-preferential genes within cells of differentiated organs might be under pressure toward alternative splicing (Kimura, 1968). Specifically, the products of alternative splicing may be detrimental to organogenesis and development. Indeed, positions within introns corresponding to known functional elements involved in pre-mRNA splicing, including the branch site, splice sites, and polypyrimidine tract, show reduced levels of genetic variation in human polymorphism data (Lomelin et al., 2010).

Alternatively spliced transcripts do not show a distinct preference of organ

According to our analysis, a transcript confers organ preference and alternatively spliced transcript(s) confer(s) the same preference or show less organ specificity. However, alternatively spliced transcripts do not exhibit additional organ distinctness or preference (Table 4). In other words, among CVI genes that show strong organ preference, none of alternatively spliced forms belong to the CVI group of the other organ. These finding suggest that multiple specificity may not be allowed for any two alternatively spliced transcripts in the CVI group. However, alternative splicing may be allowed for genes that are less detrimental to organogenesis and development, as mentioned above with regard to the observation of lower IAST in CVI compared to CVIII. Nonetheless, alternatively spliced transcripts that escape an organism’s surveillance evolve additional functionality and may be more abundant in other organs, resulting in CVII or CVIII. These aspects may place more constraint on pre-mRNA splicing of CVI genes, with lower IAST, as shown in Tables 2 and 3.

Numerous alternatively spliced transcripts have been implicated in flowering (Marquardt et al., 2014; Wang et al., 2014), circadian clock (Filichkin et al., 2010; Seo et al., 2011), and abiotic stress (Shang et al., 2017). In these analyses, alternatively spliced transcripts have a similar biological function. In the example of OsDREB2B expressed under drought conditions, two splice variants are differentially expressed in response to heat and drought stresses in rice (Matsukura et al., 2010). The transcript of OsDREB2B1 is more abundant under normal growth conditions; it contains a 53-base pair exon 2 insertion in the mature mRNA that introduces a shift in open reading frame. In plants exposed to high temperature stress, OsDREB2B2, in which exon 2 is spliced out in the mature mRNA, dominates. Genome-wide research in salt stress has revealed several features of pre-mRNA splicing in plants under abiotic stress; the expression level of these stress-induced variants is quite low relative to that of the major splicing variant(s), and most stress-induced splicing variants are associated with intron retention (Cui and Xiong, 2015; Ding et al., 2014) with a premature CDS. These observations might be in line with the allowance of alternative splicing for maximum function by avoiding expression of other alternatively spliced transcript(s) under unwanted circumstances.

Pre-mRNA splicing might be involved in balancing maximization of function and additional functionality to avoid distinct functions. For alternatively spliced variants of mRNAs, there is no way to distinguish functional differences that may be at risk of a waste of resources. Differential regulation through the components of the spliceosome might be possible, which would provide additional complexity in gene regulation. In this regard, pre-mRNA splicing might be involved in unique strategies contrasting those shown in gene families that provide direct functional distinctness by using unique promoters to respond to different signals (Freeling, 2009). As an example of a gene family, two transcripts, Os10t0505700-01 and Os11t0620300-01 (Nonspecific lipid-transfer protein 2) belong to the same ortholog group, OG5_137065 (orthomcl.org). However, Os10t0505700-01 is leaf preferential, while Os11t0620300-01 is seed preferential.

Among alternatively spliced transcripts, those with a longer CDS show higher expression

Our data also suggest that reciprocal expression is complementary among transcripts within a locus, regardless of organ specificity (Liu et al., 2011). The log2-based values of the highest and lowest expressed transcripts were 2.5–2.7 (sd = 1.9) (Supplementary Figs. S7 and S8). For CV classes, we found that longer CDS-containing transcripts largely showed higher expression (Table 5). More frequently, those with a longer CDS had the highest expression, as opposed to those with the lowest expression (1,509 compared with 945, respectively). Among the 9 alternatively spliced groups, spliced_intron was the form most frequently found among transcripts with high expression and a longer CDS (data not shown). In a survey of human transcripts, a stop codon >50 nt upstream of the final exon–exon splice junction of an inferred mRNA isoform is called a premature termination codon (Lewis et al., 2003), rendering them apparent targets for NMD. Coupling of alternative splicing and NMD may indicate that the cell possesses a large number of irrelevant mRNA isoforms that must be eliminated. When an alternative splice is introduced, the stop codon is classified, and the corresponding mRNA isoform is labeled an NMD candidate, resulting in lower expression compared to the isoform with a longer CDS.

There may be additional mechanisms for maintaining alternative splicing. Reciprocal or complementary expression patterns have been reported, whereby a copy among gene duplicates is expressed in a certain organ or tissue type and another copy is expressed to a different extent in other organs or tissue types (Liu et al., 2011). Approximately 30% of whole-genome duplicate pairs and 38% of tandem duplicate pairs show reciprocal expression patterns in Arabidopsis. This result suggests that a balancing pressure exists among these duplicated genes and that the balance is broken when the genes undergo sub-functionalization due to partitioning of the ancestral expression pattern or ‘new functionalization’ due to acquisition of a new expression pattern (Wegmann et al., 2008). Such expression balancing might be subject to alternative splicing pressure. Indeed, commonly expressed genes are under evolutionary pressure toward alternative splicing, whereas organ-preferential genes might be less subjected to alternative splicing. Favorable conditions for the occurrence of alternative splicing might be the result of constant expression and balancing of alternatively spliced isoforms under a developmental process. Based on our data, a large fraction of isoforms for transcripts with low expression and a longer CDS demonstrate reciprocal expression with other alternatively spliced transcripts. In a sense, reciprocal expression might be a common phenomenon for gene families and also for alternatively spliced transcripts of a gene.

In this study, rice ASDM and RNA-Seq were applied to differentiate alternatively spliced transcripts between representative rice organs, including the leaf, root, panicle at 1 cm, and young seed. Transcripts were classified according to organ enrichment, and the results are collectively explained in terms of evolutionary pressure. The portion of loci with multiple transcripts due to alternative splicing was higher among constitutively expressed genes than organ-preferential loci. The allowance of alternative splicing is maximized by avoiding expression of the other alternative splicing transcript under unwanted circumstances. Among loci, transcripts with a longer CDS tended to be more highly expressed than those with a shorter CDS. These genome-wide technologies might be useful for studying transcriptomes and their biological significance with regard to pre-mRNA splicing.