Analysis of Genes with Alternatively Spliced Transcripts in the Leaf, Root, Panicle and Seed of Rice Using a Long Oligomer Microarray and RNA-Seq

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).

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.

Figure 2

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.

Figure 3

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.

Figure 4

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.

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.