Comprehensive Analysis of Alternative Splicing and Functionality in Neuronal Differentiation of P19 Cells

Background Alternative splicing, which produces multiple mRNAs from a single gene, occurs in most human genes and contributes to protein diversity. Many alternative isoforms are expressed in a spatio-temporal manner, and function in diverse processes, including in the neural system. Methodology/Principal Findings The purpose of the present study was to comprehensively investigate neural-splicing using P19 cells. GeneChip Exon Array analysis was performed using total RNAs purified from cells during neuronal cell differentiation. To efficiently and readily extract the alternative exon candidates, 9 filtering conditions were prepared, yielding 262 candidate exons (236 genes). Semiquantitative RT-PCR results in 30 randomly selected candidates suggested that 87% of the candidates were differentially alternatively spliced in neuronal cells compared to undifferentiated cells. Gene ontology and pathway analyses suggested that many of the candidate genes were associated with neural events. Together with 66 genes whose functions in neural cells or organs were reported previously, 47 candidate genes were found to be linked to 189 events in the gene-level profile of neural differentiation. By text-mining for the alternative isoform, distinct functions of the isoforms of 9 candidate genes indicated by the result of Exon Array were confirmed. Conclusions/Significance Alternative exons were successfully extracted. Results from the informatics analyses suggested that neural events were primarily governed by genes whose expression was increased and whose transcripts were differentially alternatively spliced in the neuronal cells. In addition to known functions in neural cells or organs, the uninvestigated alternative splicing events of 11 genes among 47 candidate genes suggested that cell cycle events are also potentially important. These genes may help researchers to differentiate the roles of alternative splicing in cell differentiation and cell proliferation.


Introduction
The number of human genes (approximately 23,000 according to the Human Genome Project [1]) is just slightly higher than the number in lower eukaryotes such as worms. Most human exons are constitutive. However, .90% of human genes undergo in alternative splicing [2], a much higher percentage than that in lower eukaryotes [3]. Alternative splicing changes the use of exons, producing multiple types of transcripts from a single gene. Each transcript is translated into unique protein isoform. Alternative splicing enhances proteomic diversity to support complexity in higher eukaryotes [4]. Many alternative isoforms have distinct or opposing functions modified from another isoform, and play important biological roles [5,6]. Alternative splicings often occur in a spatiotemporal manner, and some are regulated by alternative splicing regulators. In the case of neural tissues, several regulators, such as nPTB, Nova1, Fox-1, and their target pre-mRNAs have been identified [7][8][9]. Understanding alternative splicing networks is very important, because splicing is one of important gene expression processes and affected on gene expression network [10].
There are five basic types of alternative splicing that involve only the differential use of splice sites: exon skipping, 59 or 39 splice site (ss) selection, mutually exclusive, and retained intron [4]. In addition, 39ss selection and/or the selection of poly-A additional signals help regulate alternative termination. Alternative promoters promote the differential use of the first exons adjacent to the transcriptional initiation sites, a process that is sometimes considered as an alternative splicing event [5,11]. Alternative splicings can occur together, causing complicated splicing patterns. The differential use of exon, generated by alternative splicing events, including alternative termination and alternative promoters, generally changes some of amino acid codes of the expressed protein.
Several technologies, such as microarrays and high-throughput sequencing, have been developed to analyze alternative splicing in genome-wide scale [reviewed in [12]]. Junction and exon arrays are the two main microarrays to investigate alternative splicing employing probes specific to the exon(s). Exon array probes are designed for each exon, while probes of junction arrays are designed for exon-exon junctions. The GeneChip Exon Array produced by Affymetrix is a convenient assay in which operations from sample preparation to scanning are performed in a standard gene expression profiling environment [13]. Although various studies using Exon Array have been published, few have efficiently extracted alternative splicings from Exon Array results and validated by RT-PCR [14][15][16][17][18][19]. Even fewer reports have applied and analyzed the validated potential alternative exons using gene ontology (GO) or pathway analyses [18,19].
P19 mouse embryonic carcinoma cells are a multi-potent cell type that can differentiate into either neural or cardio muscular cells depending on the retinoic acid (RA) concentration in the medium [20,21]. We previously reported the regulation of neural-specific splicing of the Mef2c exon b in P19 cells [22]. The transcription factor Mef2 plays important roles in muscle differentiation and synaptogenesis [23,24]. Fox-1 promoted the inclusion of Mef2c exon b via the GCAUG sequence [22]. In addition to Mef2c exon b and Fox-1, many alternative splicing events likely occur in a neuralspecific manner during the neuronal differentiation of P19 cells. These splicing events are presumably regulated by splicingregulators that are also expressed in a neural-specific manner. The alternatively spliced products, such as transcription factors, affect the transcriptional networks and function in neural events.
The goal of the present study was to comprehensively investigate neural-splicing using P19 cells. A GeneChip Exon Array was performed using total RNAs purified from undifferentiated P19 cells (Day 0), neuronal differentiated cells (Day 7), and cells from the early glial stage (Day 10). Nine filtering conditions were used for probe sets based on their annotations, estimated gene expression levels, splicing index (SI), detection above back ground (DABG) values, and predictions of alternative splicing. A total of 262 differentially alternatively spliced (DAS) candidate exons of neural splicing were extracted. Thirty exons among the 262 DAS exons were validated by RT-PCR. The RT-PCR experiments suggested that 87% of these exons were actually changed between undifferentiated and neuronal differentiated cells. This result indicated the method of the candidate exons extraction was very successful. Results from informatics approaches, including GO analysis, text-mining, and pathway analysis, also suggested that many of the candidate exons were related to neural events. In addition to gene-level informatics analyses, text-mining for the alternative isoform were performed. This informatics analysis is a new trial to identify exon-level annotations.

Exon array
Both cDNA and cRNA syntheses were performed as described in the supplier's instructions. Total RNAs from undifferentiated (Day 0), neuronally differentiated P19 (Day 7), and early glial stage cells (Day 10) were applied to a GeneChip Mouse Exon 1.0 ST Array (Affimetrix). Hybridization, washing and scanning were performed on a fluidics station and a GeneChip scanner (Affimetrix) according to the supplier's instructions. The accession number of the Exon Array data is GSE23710.

Estimation of signals
Signal data of the CEL files were sketched, normalized, and summarized with the PLIER method for exon-level intensities and with the Iter-PLIER for gene-level intensities using the Affymetrix Expression Console. The ''core'' option was used in both cases, which limits analysis to RefSeq-mapped probesets and to the mRNA of full-length CDS. The presence or absence of each probe was determined by the DABG p-values. Signal estimates of 3 time points (Day 0, 7, and 10) with biological duplicates were summarized with various annotations supplied by Affymetrix, such as probeset positions, transcript accessions, and probe sequences.

Extraction of the candidate exons and genes
Although multiple probesets are designed within the same exon, probeset intensities represent the expression levels of exons. The ratio of the probeset intensity to the gene intensity, designated as the normalized intensity (NI), was estimated in 3 time points (NI d0 , NI d7 , and NI d10 ). The splicing index (SI) was calculated by taking the log ratio (base 2) of the averaged NIs of biological duplicates [19,26]. Probesets annotated as having a high potential of crosshybridization were removed from the analysis. The SI d7~l og AvgNI d7 =AvgNI d0 ð Þ , SI d10~l og AvgNI d10 =AvgNI d0 ð Þ splicing patterns of exons were predicted from using the UCSC Blast-like alignment (blat) search tool (http://genome.ucsc.edu/ cgi-bin/hgBlat). An automatic system of alternative splicing prediction was created using the UCSC blat and NCBI BLAST (Basic Local Alignment Tool) searches. A detailed explanation of the candidate exon extraction process is described in the Results section (see also Table 1). A gene was defined as ''expressed'' when $50% of probesets of the same gene had a DABG p-value of #0.05 in both biological duplicate of the same time point [19,26]. Genes that were defined as ''expressed'' on both Day 0 and 7 or on either Day 0 or 7 were analyzed. For gene-level expression analysis, these ''expressed'' genes were categorized as differentially expressed upregulated (DEX-up) or downregulated (DEX-down). A gene was categorized in the DEX-up group if its expression level increased $2-fold and #10-fold from Day 0 to Day 7, and decreased from Day 7 to Day 10. A gene was categorized in the DEX-down group if its expression decreased $2-fold and #10-fold from Day 0 to Day 7, and increased from Day 7 to Day 10.

Semiquantitative RT-PCR
The total RNAs of P19 cells were prepared as described above. Total RNAs from mouse tissues were obtained from a commercial source (Ambion). The cDNAs were made from 2 mg of total RNAs using SuperScript III (Invitrogen) and 0.5 mg of oligo dT primer in 20-mL reaction mixture. The cDNAs are used as templates for the PCR. The PCRs were performed using the GoTaq polymerase (Promega) with specific primers. The DNA primer sequences, number of cycles, and annealing temperatures for the candidates are described in Table S1. Those for b-Actin and GluR1 were described in a previous report [22]. The PCR products were analyzed on a 6% polyacrylamide gel. The gels were stained with SYBR Green I (Takara), and the images were analyzed in LAS-3000 (Fuji Film). The sequences of the PCR products are confirmed in the 3100 DNA sequencer (ABI). Accession numbers of sequences for novel alternative variants are shown in Table 2. Densitometry was performed using MultiGauge v 3.0 software (Fuji Film). Each experiment was performed $3 times to confirm reproducibility.

Gene ontology analysis, pathway analysis and textmining
The GO analyses were performed for genes with candidate exons and for genes in the DEX-up and DEX-down groups [27]. Statistically-significant biological process terms were obtained with Fisher's Exact Test (p#0.01) using PathwayStudioH (Ariadne Genomics, Inc, Rockville, MD, USA) [27,28].

Extraction of the alternative exon candidates
Because probes of ''core'' exons was used, most of the analyzed probes were constitutive exons. Nine filtering conditions were constructed to extract the candidate exons of alternative splicings ( Table 1). The filtering procedures (including our original procedures) were arranged by considering the general features of alternative splicing and previous concepts of exon extraction [26]. Because many outermost probesets are designed at potential exons outside of the genes, the first and last probesets of each gene were removed. Since the length of most human exons is ,300 nucleotides (nt) [29], probesets with sequence of .500 nt were excluded. The long probesets mostly indicate untranslated region (UTR). Probesets with an annotation of a high potential for crosshybridization were also removed. After filtering with these conditions, 221,336 probesets of core exons of 16,661 genes were reduced to 159,552 probesets of 15,238 genes ( Table 1, third procedure).
To identify alternative splicings, genes should be expressed on both Day 0 and Day 7. The reliability of gene expression level was determined by the DABG p-values of the probesets (see Methods) [19,26]. Probesets of genes expressed on both Day 0 and Day 7 were retained on the list (Table 1, fourth procedure). In addition, exon expression level should be detected at Day 0 or Day 7. Probesets with in a DABG p-values of #0.05 on Day 0 or 7 were accepted for the next filtering procedure. These filtering conditions reduced the number of probesets to 94,780 (Table 1, fifth procedure).
The splicing index SI of a probeset represents its change in exon expression, which is adjusted according to its gene-level expression. The A-value is the average intensity between Day 0 and 7. SI d7 values and A-values of 94,780 probesets are plotted in Figure 1A. Most probesets were closely plotted in the middle (SI d7 = 0.00). This result indicated that probesets with a small absolute value of SI d7 (ABS_SI d7 ) were not related to the alternations in exon usage between Day 0 and 7. A threshold 61.35 of SI d7 (2.55-fold) was set, and probesets with $1.35 ABS_SI d7 were retained on the list (Table 1, sixth procedure).
SI of probeset represents a change of exon expression, which is adjusted with its entire gene level expression. However, Because the SI of a probeset cannot indicate the relationship between adjacent exons, the difference in SI d7 values was assessed in a local area of the transcript. We separately extracted the adjacent downstream or upstream probesets that were surely expressed (no cross-hybridization, DABG p-value #0.05). When $1 of the adjacent probesets had an ABS_SI d7 value of #0.667, the probesets were retained on the list (Table 1, seventh procedure). Theoretically the probeset indicates a border between alternative and constitutive exons after this filtering condition. If a probeset Keep a probeset that its ABS_SI d7 value is larger than its ABS_SI d10 value 262 236 * The probeset likely indicate an exon, although multiple probesets may exist in an exon.
1 Signals of these probesets were plotted using the SI d7 value ( Figure 1). Note: In the fourth procedure, genes were retained if $50% of probesets in each gene had p#0.05 in DABG. In the seventh procedure, according to the of probeset ID start position, upstream and downstream neighbor probesets were searched whose expression was detectable (DABG; p#0.05). If ABS_ SI d7 value of the neighboring probeset was $0.667, the neighbor was considered to have changed expression. A probeset was removed if both of its neighbors changed expression. If a probeset with only one neighbor, the probeset was removed if that neighbor changed expression. In the ninth procedure, in addition to comparing of the ABS_SI values, a probeset was retained when its SI d7 and SI d10 were a negative number and $1.35, respectively, or when its SI d7 and SI d10 are a positive number and #21. 35 were kept only when both adjacent probesets had an ABS_SI d7 value of #0.667, then complex splicing events would be accidentally removed.
Since most the probesets were designed to constitutive exons, it is possible that many constitutive exons are involved in the 1,107 exons that were left after the seventh filtering (Table 1). We searched these probesets' sequences in the UCSC blat and found alternative exons. The probeset was categorized as a candidate alternative exon. When a probeset sequence resulted in multiple hits for mRNA fragments in Genbank and EST and also skipped other multiple fragments in the range of each gene on genomic structure, then the probeset was categorized as an alternative exon. An automatic system was used to a double check these sequences. Contradictory exons were manually curated in the UCSC blat search again. Approximately 28% of the probesets (309 probesets) were predicted to be alternative exons (Figure 1B &  Table 1, eighth procedure).
Finally, the SI d7 and SI d10 values were compared. A probeset was retained if its SI d7 was larger than SI d10 when SI d7 was positive (or smaller than its SI d10 value when SI d7 was negative) ( Table 1, ninth procedure). Altogether, 262 probesets containing all alternative splicing types were extracted (Table S2 & Figure 1C). These probesets were identified as DAS exons that potentially had undergone alternative splicing in a neural-specific manner.

Semiquantitative RT-PCR for the candidate exons
To examine whether the potential DAS exons are actually involved in alternative splicing during the neuronal differentiation of P19 cells, semiquantitative RT-PCR was performed. Semiquantitative RT-PCR is more appropriate than qPCR in the analysis of alternatively spliced products. Among 262 DAS exons, 30 exons were randomly selected (Figure 2A & Figure S1A). Typical patterns for the DAS events in P19 cells are shown in Figure 2B. Densitometric analysis showed alternative splicing of 10  Figure S2B). The expressions of these 26 exons increased or decreased consistent with the directions indicated by SI d7 . The changes of these 26 exons on Day 7 were greater than those on Day 10, consistent with the results in Table 1 (ninth procedure). These results show that neuronal cell stage-specific alternative splicing occurred in these 26 of 30 exons, and suggested that 87% of the 262 DAS exons were involved in neural splicing. Seven alternative splicing events had novel splicing patterns (Kif1b, Dlg3, Plekha5, Atf2, Clta, Fam113a, and Zfp326; Table 2) in mouse. The patterns included a novel exon, novel splice sites, and novel combinations. Except for the novel Dlg3 variant, the other 6 variants were detected in adult mouse brain ( Figure S2C). Among 26 alternative exons, the percentages of alternative variants of 17 exons were changed .2-fold in the brain compared to other tissues, which showed expression patterns similar to those in the brain ( Figure 2C & Figure S1C). Although the differentiation of P19 cells is a model system, this result implied that alternative splicings of more than half of DAS exons had similar regulations between neuronal differentiated P19 cells and brain.

Gene ontology analysis of the candidate exons
To test whether the 262 DAS exons of 236 DAS genes were involved in neural events, GO analysis was performed using the Biological Process category. A total of 72 GO terms (p#0.01, Fisher's Exact test) were obtained from the 236 DAS genes. These terms were classified into 10 groups: neural-related process (there were no neural-related terms in other terms), differentiation and development, cytoskeleton and cell adhesion, signaling, post-  translational regulation, transcription, cell cycle and proliferation, cellular transport, apoptosis, and others. The neural-related process group was the most abundant, accounted for 20.8% of all terms, followed by differentiation and development (12.5%), cytoskeleton and cell adhesion (12.5%), and signaling (9.7%) ( Figure 3A).
The 262 DAS exons were categorized into two groups according to whether the SI d7 value was $+1.35 (DAS-up; 153 genes) or #21.35 (DAS-down; 88 genes). The percentage of neural-related GO terms in the DAS-up and DAS-down groups were 17.0% and 21.2%, respectively ( Figure 3B). This result suggests that the protein isoforms translated from these alternatively spliced transcripts, which were either included or excluded exons, were involved in neural events.
Affymetrix Expression Console calculates the expression level of each gene based on the intensities of probesets of the Exon Array. We analyzed the gene expression of the same samples and obtained the differentially expressed (DEX) genes in P19 neuronal cells. The expression of 1666 DEX-up genes increased $2-fold and #10-fold from Day 0 to Day 7, and decreased from Day 7 to Day 10 (Table S3). The DEX-up gene group had 194 GO terms, including neural-related process (20.1%) and differentiation and development terms (18.6%) ( Figure 3B). The DEX-up and DAS gene results were quite similar, suggesting that genes in both groups were related with neural-events ( Figure 3C).
A total of 550 DEX-down genes were obtained, whose expression levels decreased $2-fold and #10-fold from Day 0 to Day 7 and increased from Day 7 to Day 10 (Table S4). The DEXdown genes had 107 cell process GO terms. The percentage of neural-related processes in all 107 terms was only 2.8%, and those of cellular transport and apoptosis were 11.0% and 5.6% ( Figure 3B). The GO terms of DEX-down genes showed unique patterns compared with those of the DAS, DAS-up, DAS-down, or DEX-up genes, suggesting that the DEX-down genes are not responsible for neural-events ( Figure 3C).

Potentially important exons extracted by pathway analysis and text-mining
Text-mining and pathway analysis were performed for the DAS genes using MedScanH text-mining technology in PathwayStudioH (Ariadne Genomics) [28]. Of the 236 DAS genes, the functions (cell process) 151 genes have been reported in multiple publications. The remaining 85 genes were categorized as 'unknown genes'. Of the151 DAS genes with known functions, 66 were categorized as 'well-known genes' in neural processes since their functions in neural cells or organs were reported previously. The remaining 85 were categorized as 'functional DAS genes', having little reported evidence of neural regulatory roles ( Figure 4A & 4C). This suggests that many DAS genes are related with neural events, as already shown by GO analyses. Because 66 well-known genes has been studied in neural cells or organs, further research about functional differences between alternative isoforms will probably clarify their functions in non-neural cells or organs. It will be very difficult to analyze alternative isoforms in the 85 unknown genes, since there are too few clues about their functions. Therefore, we focused on the 85 functional DAS genes.
In order to extract potentially important genes from the 85 functional DAS genes, we focused on DEX-up genes since DEXup genes and DAS genes might affect on their expressions each other. As shown in GO analysis (Figure 3), in contrast to the DEXdown genes, many of the DEX-up genes seemed to regulate neural events. A total of 279 cell process terms closely related to DEX-up genes were extracted from the literature using text-mining and PathwayStudioH database ( Figure 4B). Ninety overlapping terms of DEX-up and DEX-down groups were subtracted from the 279 terms. The remaining 189 cell process terms are presumed to be responsible for events in the neuronal P19 cells.
We next extracted biological relationships between the 189 cell processes and 85 functional DAS genes. Forty-seven DAS genes were linked to these cell processes ( Figure 4A & Table S5). These 47 genes are potentially important, because isoforms created by neural splicing may acquire novel functions that have not been reported in neural cells. Eleven genes out of the 47 genes were involved in cell cycle-related events, processes in the G1 and/or S phase ( Figure 5). This result is reasonable, since cell proliferation and cell differentiation are very closely related.

Informatics and experimental analyses of splicing isoforms
In addition to the above gene-level analyses, we tried to perform exon-level informatics analyses. Among 236 DAS genes, 49 genes have been reported in articles about isoforms or splicing variants of previously investigated transcripts ( Figure 4C & Table S2). Thirtyfour genes were overlapped between the 66 well-known genes and 49 genes with known splicing (Figure 4A & 4C). The exons of 15 of the 34 genes were consistent with previous reported variants (Table 3). Of the 15 genes, distinct functions of the alternative isoforms in each gene were confirmed in 9 genes, including Gnao1 (Table 3). These might mediate neural events as a result of alternative splicing. For instance, results of semiquantitative RT-PCR and the SI value of Gnao1 showed that the Gnao1 B isoform (189 nt) changed to the A isoform (124 nt) in P19 neuronal cells ( Table 2 & 3). This result is consistent with a previous report showing that GNAO1 A isoform, but not the B isoform, is required for the light response of retinal bipolar cells [30].
Among the 49 genes with known splicing, 15 genes overlapped with 85 functional genes ( Figure 4C). In this study, only the splicing of the Abi2 gene in these 15 genes was analyzed by RT-PCR (Table 2). Although a detailed expression analysis of splicing has not been reported, the RT-PCR and SI results indicated that the skipping isoform (149 nt) of Abi2 changed to an inclusion isoform (329 nt) during the neuronal differentiation of P19 cells (Table2). Interestingly, one previous study has reported a neural function of Abi2 gene, showing that Abi2 is involved in 'memory' and 'learning' [31]. To test whether Abi2 splicing is related to this higher brain function, alternative splicing was studied by brain dissection, which revealed that the relative amount of the inclusion isoform was increased in the frontal cortex compared with other tissues ( Figure S3). This preliminary result implies that the neural splicing regulation of Abi2 is quite important for higher brain functions.

Discussion
One important aspect of Exon Array projects in wet laboratories is that the extracted candidates primarily consist of actual alternative exons validated by RT-PCR. However, there have only been a small number of projects in which the DAS exons were efficiently extracted and validated by RT-PCR [12]. In the present study, 87% of the 262 DAS exons were found to change in the neuronal stage of P19 cells by semiquantitative RT-PCR (Table 2). This score indicates that the filtering conditions used efficiently extracted DAS exons from signals using Exon Arrays. Another important aspect in Exon Array projects is the ease of modifying the filtering conditions. To obtain more reliable DAS exons, the threshold of the ABS_SI d7 values in the sixth procedure in Table 1 should be increased (e.g., from $1.35 to $3.0). A higher ABS_SI value likely indicates a remarkable change of alternative splicing in the RT-PCR experiments ( Table 2). Probesets with high ABS_SI d7 values indicated alternative exons in a ratio of approximately 40,50% in the UCSC blat search ( Figure 1B). There are probably many alternative splicings in DAS exons with high ABS_SI values. Indeed, these exons can be efficiently validated by RT-PCR in a neural-specific manner.
Although the ABS_SI d7 threshold can be variable, we think that a threshold of ABS_SI d7 $1.35 in the sixth procedure of Table 1 is appropriate. According to the approximation curve, the percentage of predicted alternative exons decreased in the area of high ABS_SI d7 values. The percentage of predicted nonalternative exons plateaued at an ABS_SI d7 value of approximately 1.35 ( Figure 1B). This is because the potential DAS exons accumulated at an ABS_SI d7 of $1.35. The inflection point between decreasing range and a plateau state (e.g., in in Figure 1B) is probably suitable for use as the ABS_SI threshold. It is also preferable to use accessible software like the Expression Console and Excel. Open access homology search tools such as NCBI BLAST or UCSC blat are useful in the prediction of alternative exons, and reduce the expense of the Exon Array analysis. The extraction method PathwayStudioH was used on 236 DAS genes. The functions of 85 genes were not reported in multiple studies, and these were categorized as unknown genes (*). Sixty-six genes were reported in studies on the functions in neural cells or organs, and were categorized as well-known genes. Another 85 genes were categorized as functional genes. The approach for the 85 functional genes is described below. Among the 236 DAS genes, 49 were reported in articles on 'splicing variant' or 'isoform'. Thirty-four genes overlapped with the 66 well-known genes. Further confirmation was performed (Table 3). (B) Number of cell processes in pathway analysis of DEX-up and/or DEX-down genes (see Figure 3C). Orange, blue, and yellow indicate the number of biological processes from DEX-up genes, DEX-down genes, or both, respectively. Since DEX-up genes are primarily important for neural-events and DEXdown genes are not (see Figure 3), it is predicted that 189 biological processes take place in the neural differentiation. Forty-seven of the 85 functional DAS genes from the remaining 170 genes, whose functions matched with the 189 specific biological processes of the DEX-up group, were extracted. (C) Number of genes in Figure 4A. Upper box shows 66 well-known (purple), 85 unknown (*), and 85 functional genes. Among the 85 functional genes, 38 genes were not linked (green) and 47 genes were linked with the 189 processes and were categorized as potential genes for neural regulation (orange). Among the 236 DAS genes, 49 genes (red) were found in articles concerning 'splicing variant' or 'isoform' (lower box) and 34 genes were found in articles of neural functions. doi:10.1371/journal.pone.0016880.g004 presented here could become the standard method for extracting DAS exons in wet laboratories, because it is useful in projects that employ a small number of samples and promotes the comprehensive analyses of alternative splicing events.
We detected novel splicing patterns in 7 alternative splicing events among 30 randomly selected exons from 262 DAS exons ( Figure 2 & Table 2). Therefore, there are still many undetermined alternative splicing events. It is possible that probesets that were removed in the prediction of alternative splicing (Table 1; eighth procedure) also contain novel alternative exons. The latest software or the scoring of exons using exonic splicing enhancer (ESE) and exonic splicing silencer (ESS) may accurately predict DAS exons [32,33]. Such scoring is especially important for aberrant splicing events, which are frequently observed in cancer cells [34].
In addition to these sequence data, the results of informatics analyses showed that only 49 of the 236 DAS genes have been reported in previous articles of 'alternative splicing variants' (Figure 4). Although this informatics analysis is a new trial to identify exon-level annotations, it is likely that insufficient researches has been performed on splicing isoforms. Moreover, there is no assembled database that shows the relationship between a function and each distinct splicing isoform, even though databases for functional analysis are substantial in annotations of genes and proteins. Expression profiles containing alternative splicings are accumulating in databases. Therefore, exon-level informatics analyses will only increase in importance.
Forty-nine of the DAS genes were splicing 'isoform' or 'variant'. Of these, 34 genes had well-studied functions in neural cells or organs. Nine of the 34 genes have a reported distinct function for their alternative isoforms (Table 3). In particular, one Gnao1 isoform, GNAO1 A, is required for light response [30]. Gnao1 B (189 nt) was a major product in undifferentiated P19 cells and changed to Gnao1 A (124 nt) in neural differentiated P19 cells (Table 2 & Figure S1). Among adult mouse tissues, Gnao1 A was abundant in brain or brain dissections ( Figure S1 and S3). Therefore, it seems that the splicing regulation of Gnao1 is primarily involved in neural differentiation. Although Abi2 is one of the 49 genes with known splicings, it was categorized with the 85 functional genes. Distinct functions for the Abi2 isoforms are unknown. The skipping isoform (149 nt) of Abi2 was replaced with the inclusion isoform (329 nt) in neuronal differentiated P19 cells (Table2 & Figure S1), and this isoform was largely detected in the frontal cortex of the adult mouse ( Figure S1 & S3). Alternative splicing regulation of Abi2 may be important for higher brain functions, based on the present results and the previous finding that Abi2 mediates 'memory' and 'learning' [31].
Along with Abi2, the 85 functional genes may be the most important group to investigate alternative splicings, since they are substantially unknown in neural cells or organs. To extract potentially important DAS genes from the 85 functional genes, both alternative splicings and gene expressions were examined. Informatics analyses showed that many DAS genes and DEX-up genes were involved in neural-related events ( Figure 3C). Based on the GO biological process terms, we assumed that the DEX-up genes, but not the DEX-down genes, were responsible for events in P19 neuronal cells ( Figure 3C). The protein expression levels of the DEX-down genes are probably reduced in P19 neuronal cells. Reductions or disruptions of DEX-down proteins may prevent them from functioning in the differentiated cells.
Some of the DAS genes likely directly or indirectly affected on expressions of the DEX-up genes, and some of the DEX-up genes likely affected on alternative splicings of DAS genes [10]. Therefore, from the pathway analysis, the links between 189 cell processes of DEX-up genes and 85 functional genes out of 236 DAS genes were examined ( Figure 4A). Eleven of the 47 presumably important DAS genes mediated the G 1 /S transition in the neuronal differentiation of P19 cells ( Figure 5). The regulations of alternative splicing in the 11 DAS genes may play important roles in the cell cycle regulation. This informatics result is very reasonable, since cell proliferation and cell differentiation are highly associated. Some alternative events have been also investigated, but there are still unknown frontiers in alternative splicing research (Table 3). PathwayStudioH was used to find linkages between the 85 functional genes and 189 biological processes (see Figure 4A). Forty-seven DAS genes and their linked to biological processes are shown. Cell cycle-related functions were found (orange & red), especially those relate to the G 1 and/or S phase (red Curation confirmed that 34 DAS genes (40 DAS exons) were reported to be variants or isoforms (see Figure 4A). The probeset ID sequences of 40 DAS exons and previous literature on the variants of these genes were compared to confirm whether the sequence and the variants indicated the same isoforms. The literature was also searched for the distinct functions of the isoforms. doi:10.1371/journal.pone.0016880.t003 This study obtained the expression profiles of genes and alternatively spliced transcripts. These profiles included key factors and their potential targets in the transcriptional and alternative splicing networks. Among the DEX-up and DAS genes, we observed 20 and 14 RNA processing-related genes, respectively. These genes may play key roles in neural splicing regulation. For instance, RBM9 was included in both the DAS and DEX-up groups. The splicing activity of the RBM9 isoform, which is indicated by the DAS exon, is reduced compared to that of another isoform [35]. This study also focused on the cell cyclerelated genes ( Figure 5). SRp38/NSSR1, which we previously detected during neuronal differentiation in P19 cells [22], is a known splicing repressor via dephosphorylation in M-phase specific manner [36]. Nucleophosmin (NPM) is another cell cycle-dependent splicing regulator [37] that is phosphorylated during the G 1 /S transition. These splicing regulators and their potential targets may be very important in future studies clarifying the alternative splicing networks involved in the neuronal differentiation of P19 cells.  Figure 2A, 2B, S1) followed by densitometric analysis are shown in the graphical representation. The amount of each PCR product was divided by total amount in each lane. Percentages are shown in the histogram. Bars indicate the standard error. Numbers on the right side of the graph indicate the length of PCR products. Numbers at the bottom of the graph show the time course: Day 0 (undifferentiated stage), Day 7 (neuronal stage), and Day 10 (early glial stage). (TIF) Figure S3 Alternative splicings in brain dissection. Semiquantitative RT-PCR experiments were performed using total RNAs of Skeletal muscle, brain, frontal cortex, hippocampus, cerebellum, spinal code (left). Total RNAs of brain dissection are commercially available (Takara). Gnao1 and Abi2, but not Kif21b, are categorized among the 49 genes with known splicings. Distinct functions of the Gnao1 isoforms, but not of the Abi2 isoforms, were previously reported, but not Abi2. b-Actin was used as a control. Schematic representations of alternative splicings are shown (right). Boxes and middle lines indicate exons and introns, respectively. Gray indicates a possible alternative exon. Arrows indicate locations of the primer annealing sites. Numbers indicate the length of PCR products. (TIF)