Cell-type-resolved alternative splicing patterns in mouse liver.

Alternative splicing (AS) is an important post-transcriptional regulatory mechanism to generate transcription diversity. However, the functional roles of AS in multiple cell types from one organ have not been reported. Here, we provide the most comprehensive profile for cell-type-resolved AS patterns in mouse liver. A total of 13,637 AS events are detected, representing 81.5% of all known AS events in the database. About 46.2% of multi-exon genes undergo AS from the four cell types of mouse liver: hepatocyte, liver sinusoidal endothelial cell, Kupffer cell and hepatic stellate cell, which regulates cell-specific functions and maintains cell characteristics. We also present a cell-type-specific splicing factors network in these four cell types of mouse liver, allowing data mining and generating knowledge to elucidate the roles of splicing factors in sustaining the cell-type-specialized AS profiles and functions. The splicing switching of Tak1 gene between different cell types is firstly discovered and the specific Tak1 isoform regulates hepatic cell-type-specific functions is verified. Thus, our work constructs a hepatic cell-specific splicing landscape and reveals the considerable contribution of AS to the cell type constitution and organ features.

1. Introduction

Alternative splicing (AS) plays a major role in increasing transcriptome variability and proteome diversity., It is estimated that 95% of multi-exon genes in human undergo AS, and the resulting isoforms from one gene could present identical, similar or even opposite protein functions., AS variants are variably expressed between different tissues and cell types, and most of the recent studies focused on the varying effects of AS in different species or tissue types.,, AS is highly prevalent across distinct species, and the divergences in AS patterns significantly contribute to shaping species-specific differences. Furthermore, it was suggested that proteins with tissue-specific isoforms were likely to occupy key hubs in protein interaction networks.,,, Although the roles of AS in tissue specificity have been well-studied, splicing differences at the cellular level and in multiple cell types from one single organ were almost neglected to date.

Thus, we expected to establish splicing profiling of four cell types from the mouse liver as a paradigm of exploring cell-type-specific splicing patterns. The liver is the largest metabolic organ of the mammals and participates in multiple biological functions, the majority of which are performed by parenchymal cells, also named hepatocytes (HCs), occupying approximately 80% of the liver volume. The remainders are nonparenchymal cells (NPCs), including the liver sinusoidal endothelial cells (LSECs), the Kupffer cells (KCs), the hepatic stellate cells (HSCs) and several rare cell types. Despite the small fraction in total liver volume, NPCs also perform significant functions.

Recently, although some cell-type-specific proteomic analyses had performed in mouse liver, that globally investigated the proteome profiles of individual cell types and possible mechanisms of intercellular crosstalk, we still want to know how this proteomic complexity was generated and how AS contributed in such processes. Given that tissue-specific exons play crucial roles in determining tissue identity, we wondered whether this feature also expands to the cellular level in the liver and how AS regulates the biological functions of isoforms in different cell types. Thus, in our study, we obtained highly purified liver cells and applied RNA-Seq to build the first dataset of cell-type-specific transcriptomic profiling. Next, we comprehensively compared the splicing patterns between the four hepatic cell types and analyzed their biological characteristics. Furthermore, we identified the cell-type-specific splicing factors that dominated the splicing patterns and functional characteristics in respective cell types. Last, we took HCs-specific Tak1 isoform as an example to demonstrate the effect of splicing switching on cell-type specificities. Overall, our results revealed a hepatic cell-type-specific splicing profile that played a considerable role in determining cellular-specific functions.

2. Materials and methods

2.1. Isolation of four hepatic cell types

C57BL/6 male mice (8–10 weeks old) were killed, and their livers were excised by operation according to the ethics committee guidelines of the Beijing Institute of Radiation Medicine. To reduce the biological variation among individuals, the liver samples from three mice were pooled for the sequencing of each cell type. All procedures were reported before.

2.2. RNA-seq data generation and processing

Total RNA of four hepatic cell types was isolated with Trizol Reagent (Invitrogen) according to the manufacturer’s protocol. RNA quality was assessed using Agilent 2100 Bioanalyzer. RNA-seq library was prepared using Illumina’s reagent and protocol. Paried-end 90 bp sequencing was performed on an Illumina HiSeq 2000 genome analyzer. Clean reads were mapped to the mouse GENCODE M4 (GRCm38) using Tophat2 (version 2.0.10) with a mate inner size of 20. FPKM were then calculated as relative abundance value for each gene and transcript by Cufflinks (version 2.2.1) supplied with a protein-coding gene model annotation file in GTF format (GENCODE M4).

2.3. Differential analysis of AS between four cell types

MATS was applied to identify the differential alternative splicing between four cell types with options ‘-t paired –len 90 –a 8 –c 0.0001 –expressionChange 20’. Significant different events were filtered under the threshold of FDR < 0.05 for five common splicing classes. Percent spliced-in (PSI) values for each alternative exon were extracted from MATS outputs.

2.4. Feature analysis of CEEs

Cell-type-enriched exons (CEEs) required statistically significant differences (ΔPSI > 0.25 and FDR < 0.05) between at least two cell types. All exons were mapped to annotated proteins (Ensembl release 78) using the method of Ensembl Perl API. Protein domains of Pfam were extracted from Ensembl annotations for mapped exons. IUPred was applied to predict protein disorder information. The fraction of segments classified as disordered by the IUPred cutoff of 0.4. ANCHOR was used for predicting protein binding regions. PTM sites were collected from the UniProt database after mapping the Ensembl protein ID to UniProt accessions. The statistically significant differences were assessed by χ2 test.

2.5. Splicing validation by RT-PCR

Total RNA was reverse transcribed using the PrimeScriptTM RT Master Mix (TaKaRa Bio Technology), and PCR was performed using Premix TaqTM (TaKaRa TaqTM Version 2.0 plus dye, TaKaRa BioTechnology). The PCR products were separated by electrophoresis in 2% agarose gels and detected using the Gel Doc EZ System (Bio-Rad, USA). Sequencing of the PCR products was performed by BGI.tech (Beijing, China). All the primers were listed in Supplementary File S1.

2.6. GO and pathway enrichment analysis

All enrichment analysis was performed by the DAVID online tool (http://david.abcc.ncifcrf.gov/). Enriched functional terms mainly included three ontology terms (BP_FAT, CC_FAT and MF_FAT) and the KEGG pathway. Significantly enriched terms were determined under FDR < 0.05. The GO and pathway networks for splicing factor-regulated target genes were built using Enrichment Map in Cytoscape 3.2.0 with the parameters of P < 0.001 and FDR < 0.1.

2.7. Oligonucleotide synthesis

The splice switching oligonucleotides (SSOs) were synthesised by Sangon Biotech (Shanghai) Co., Ltd. The sequences were listed in Supplementary File S1.

2.8. Cell culture

The AML12 (alpha mouse liver 12) cell line was purchased from ATCC (CRL-2254), which was established from HCs from a mouse (CD1 strain, line MT42) transgenic for human TGFα. These cells exhibit typical HC features, being a suitable alternative option for the mouse primary HCs. The cells were cultured in DMEM (GIBCO, USA) with 10% of fetal bovine serum (GIBCO, USA) and supplied with 1% of insulin-transferrin-selenium (GIBCO, USA) and 40 ng/ml dexamethasone (Sigma, USA).

2.9. The knockdown of the splicing factors

The siRNAs used for the knockdown of the indicated splicing factors were designed and synthesised by Sangon Biotech (Shanghai) Co., Ltd. After 72-hs transfection, the total RNA of AML12 cells was extracted and the AS of Tak1 gene was analyzed. The siRNA sequences were listed in Supplementary File S1.

2.10. Enrichment analysis of RNA-binding sites

RNA-binding sites of splicing factors were mapped to pre-mRNA using the RBPmap online tool (http://rbpmap.technion.ac.il/) with the parameters of a high stringency level and conservation filter. The input sequences included alternative exon regions and their neighboring intron regions (300 bp). Enrichment of RNA-binding sites was performed using the hypergeometric test, and FDR < 0.05 was used to assess significance.

2.11. Data presentation

Unsupervised hierarchical clustering was performed by Perseus software with default options (distance = ‘Euclidean’). Cumulative distribution curves were visualized using the ggplot2 package in R. Other statistical tests were also performed by R software.

3. Data availability

The RNA-Seq data from this study have been deposited into the NCBI Short Read Archive database under study accession number SRP033468.

4. Results

4.1. A global view of gene expression and splicing patterns

In total, we identified 12,947 genes with FPKM ≥ 1 (Fig. 1A and Supplementary Table S1), each cell type had a series of cell-specific genes, while the majority (65.6%) were detected in all. The detected gene number of HCs, LSECs, KCs and HSCs was 9,967, 10,573, 10,964 and 11,767, respectively. The parenchymal cells, which participate in the liver’s main biological functions, had the least number of total genes but more high-abundance (FPKM ≥ 100) genes than NPCs (Fig. 1A). RT-qPCR of 20 randomly selected genes was performed in additional biological samples, the expression patterns of which were similar to the patterns obtained from the RNA-Seq, indicating the high correlation and good reproducibility between multiple samples (Supplementary Fig. S1A).

Figure 1

Cell-type-resolved transcriptome expression and AS profiling in the mouse liver. (A) The overlap of identified genes and the number of detected genes at different FPKM levels. (B) The overlap of AS events and the number of five splicing types. (C) Left: Percentage of alternatively spliced genes in all identified genes. Right: Percentage of alternatively spliced genes with a different number of splicing variants. (D) The correlation analysis of PSI values calculated from RNA-Seq and RT-PCR.

The genome-wide extent of AS was analyzed by searching against known splicing junctions. We focused on five common types of AS (Supplementary Fig. S1B). The AS events identified in our data achieved an average coverage of 81.5% of all known events in the Ensembl database (Supplementary Fig. S1C), consistent with the fact of that the liver was frequently regulated by AS (Supplementary Fig. S1D and Table S2). Most of the AS events existed in all four cell types (Fig. 1B), and in our data, 46.2% of multi-exon genes underwent AS, half of which had three or more AS variants, indicating the distinguishing contribution of AS to the proteome diversity in mouse liver cells (Fig. 1C). Consistent with previous studies in mammals, the skipped exon (SE) was the main AS type in mouse liver, followed by A3SS (Alternative 3′ splice site) and A5SS (Alternative 5′ splice site), and a low frequency was observed for retained intron and mutually exclusive exons (Fig. 1B and Supplementary Fig. S1C). All types tended to occur equally in different cell types, with a rank from 69.8% to 72.5% for SE (Fig. 1B). Although LSECs didn’t have the largest gene number, it had more AS events than the other three cell types, which indicated a precise regulation on the splicing level beyond the transcriptional regulation in LSECs.

To confirm the identification of AS from RNA-Seq data, randomly selected AS events were validated by RT-PCR (Fig. 1D and Supplementary Fig. S1E), and the results displayed a high Pearson correlation (R2 = 0.808) of PSI values from RT-PCR and RNA-Seq.

4.2. Association of gene functions with AS

We classified the AS-regulated genes detected in all cell types into several function categories using MetaCore software. In general, the genes that most frequently regulated by AS were those encoding ‘regulators (GDI, GAP and GEF)’ (62.2%), ‘transcription factor’ (56.7%) and ‘generic kinase’ (56.2%). While for ‘receptor ligand’, ‘GPCR’ and ‘ligand-gated ion channel’, only 33.2, 29.7 and 24.1% of genes were regulated by AS (Fig. 2A), indicating a putative mechanism that the distinct regulation of AS on genes of different function categories. Besides, AS did not happen equally across different cell types even for the genes in the same function categories. We found that the genes of ligand-gated ion channels did not undergo AS equally across cell types, the genes were more often regulated by AS in LSECs (45.8%) than in HCs (16.7%) (Fig. 2A). The statistical significance of spliced genes proportion for each function category between different cell types was shown in Supplementary Fig. S2A. In the liver, several diseases were caused by ion channel dysfunction. Different ion channels might have specialized biological significances in one particular cell type, yet no relevant investigation has been reported so far.

Figure 2

Functional categories of spliced genes across cell types. (A) Altered percentages of alternatively spliced genes for different functional categories. ‘General’ column represents the combination of spliced genes in all four hepatic cell types. (B) Differentially enriched function terms revealed by enrichment analysis in four hepatic cell types. (C) Unsupervised hierarchical clustering for PSI values of AS events. Enrichment analysis shows the corresponding functional annotation for spliced gene clusters.

We next analyzed the enriched functional terms of AS-regulated genes in distinct cell contexts. The AS-regulated genes in HCs were highly enriched in multiple metabolic processes, such as the organic substance metabolic process and the primary metabolic process. Similarly, in NPCs, the AS-regulated genes also tended to play roles in cell-type-specific functions, just as cytoskeleton organization and immune system process (Fig. 2B and Supplementary Fig. S2B). To our surprise, plenty of AS-regulated genes functioning in regulating biological processes, such as ‘positive regulation of metabolic process’, were significantly enriched in NPCs instead of in HCs (Fig. 2B). To eliminate the potential effects of low abundance genes on enrichment analysis, we separately chose the 1000 and 2000 most abundant AS genes in each cell type to reanalyze the enriched biological processes and obtained the consistent trend (Supplementary Fig. S2C and D). Besides, all the four cell types showed an almost uniform cumulative distribution and the bulk of transcription was dominated by a few high abundance genes (Supplementary Fig. S3A), demonstrating they would be identically affected by detection sensitivity for low abundance genes. These findings confirmed the important roles involving regulation of metabolic process in NPCs and provided new evidence of the cooperation between HCs and NPCs on the isoform level.

Furthermore, we explored the association of AS with gene functions using the quantitative data of splicing in different cell types. We firstly defined the CEEs as described in the Materials and Methods section. The unsupervised hierarchical clustering revealed that these exons exhibited obvious splicing switches between different cell types (Fig. 2C). Functionally, the genes with HC-enriched exons were linked to ‘regulation of kinase cascade’ and ‘phosphoprotein’, largely reflecting the identities of the HCs and differing them from NPC-enriched exons. In NPCs, each cluster also revealed significant enrichment in specific GO terms, including ‘cell surface receptor linked signal transduction’ (LSEC-enriched), ‘regulation of immune effector process’ (KC-enriched) and ‘membrane organization’ (HSC-enriched). These CEE-contained genes reflected the core functions of respective cell types, as the visualized global diversity of all genes with FPKM ≥ 1 and the functional enrichment analysis (Supplementary Fig. S3B and C).

4.3. Functional properties of CEEs

Tissue-enriched exons had exhibited multiple properties that distinguished them from constitutive exons., Thus, we wondered if the CEEs in mouse liver also had similar properties and how AS affecting protein functions. In our data, although some of the CEEs destroyed the functional protein domain (e.g. the immunoglobulin domain 1 of Fgfr1) when excluded from the pre-mRNA, the CEEs had a significantly smaller fraction of overlapping to protein domains compared with other alternative exons and constitutive exons (Fig. 3A;P<2.07E-15). In contrast, an increased fraction overlapping to predicted disordered regions was observed (Fig. 3B;P<4.56E-11). The disordered regions had important regulatory effects on protein–protein interaction (PPI) networks,, and consistently, the CEEs were significantly enriched in regions that involved in PPI regulation (Fig. 3C;P<7.16E-04). In addition, these CEEs contained more post-translational modification (PTM) sites than others (Fig. 3D, adjusted to the length of the exons; P<0.022). Together, in mouse liver, our results suggested that the CEEs might regulate isoform functions through affecting PPIs and modifying PTM process, consistent with the previous reports on the tissue-specific exons.,

Figure 3

Characterization of CEEs. (A–D) Percentages of three groups of exons encoding (A) protein domains, (B) disordered regions, (C) binding regions and (D) PTM sites. Significance was calculated using the χ2 test. (E) Different PSI values of alternative exons between HCs and NPCs. Exons in circles are highlighted as examples. (F) Variants-switching validation between HCs and NPCs for Tak1, Sidt2, Lrch3 and Mbnl1.

A scatter diagram illustrated all the alternative exons that exhibited markedly different PSIs between HCs and NPCs (Fig. 3E). The circles shown in Fig. 3E highlighted the flexibility of the splicing regulatory machinery, with a large number of exons predominantly included in one cell type but excluded in others (just like exon 12 of Tak1, with a PSI of 0.886 in HCs, but 0.051, 0.059 and 0.117 in LSECs, KCs and HSCs, respectively; Supplementary File S3). We used RT-PCR assays to validate individual AS events detected by RNA-Seq (Fig. 3F). The validated examples included (but not limited to) the genes that functioning in regulating the MAPK pathway (Tak1), splicing regulation (Mbnl1), glucose homeostasis (Sidt2) and a gene with putative immune functions (Lrch3). Most of the exons correlated well with our RNA-Seq data, exhibiting a strong cell-type-regulated AS pattern.

4.4. Regulation of as by splicing factors

We were interested in how these CEEs generated and which particular splicing factor (SF) might dominate the splicing patterns. First, we collected 94 splicing factors from RBPmap, 78 of which (with 114 experimentally confirmed binding motifs) were expressed in our data with FPKM ≥ 1 and displayed variable cell-type specificities (Fig. 4A). The most significant diversity was observed between HCs and NPCs, and each cell type had a set of highly expressed splicing factors, as listed in Supplementary File S4.

Figure 4

Identification of cell-type-specific splicing factors. (A) Differential expression patterns of splicing factors revealed by unsupervised hierarchical clustering. (B) Differential expression abundances of four cell-type-specific splicing factors in both mRNA and protein level.

In both mRNA and protein levels, we identified the high-expressed splicing factors in each cell type, such as Esrp2 in HCs, Rbfox2 in LSECs, Srsf9 in KCs and Rbms3 in HSCs (Fig. 4B, proteomics data came from the literature 17). We used RBPmap to search the target genes of these splicing factors, followed by the functional enrichment analysis. These splicing factors controlled a large proportion of the cell-type-specific AS events, and the target genes were highly enriched in cell specialized functions. For example, Esrp2 was highly expressed in HCs but not in three other NPCs, and the targets of Esrp2 were significantly enriched in terms of liver specialized functions, such as lipoprotein and coenzyme metabolism (Fig. 5A). The target genes of the LSEC-specific splicing factor Rbfox2 were enriched in cytoskeleton organization and cell adhesion (Fig. 5B). Rbfox2 was reported to splice the Tak1 pre-mRNA to the short variant Tak1-A instead of the full-length Tak1-B, consistent with the distribution of Tak1 variants in HCs and LSECs. Similarly, in Kupffer Cells, we found that the target genes of Srsf9 were highly enriched in multiple immunity-related biological processes and pathways (Fig. 5C and Supplementary Fig. S4A). Furthermore, although little was reported about the functions of Rbms3 in the liver, our data showed that Rbms3 were exclusively high-expressed in HSCs, regulating a range of biological processes, especially those HSC-specific functions like cell-matrix adhesion (Fig. 5D).

Figure 5

Functional regulation by splicing factors. Enriched biological processes of spliced genes regulated by (A) Esrp2 in HCs, (B) Rbfox2 in LSECs, (C) Srsf9 in KCs and (D) Rbms3 in HSCs.

Besides analyzing the differential expressed splicing factors, we combined the trans-acting splicing factors and cis-element (binding motifs) surrounding CEEs and built a splicing regulatory map, in order to further identify the important splicing factors in each cell type (Supplementary Fig. S4B). Specifically, LSECs exhibited the most significant enrichment of RNA-binding motifs of splicing factors across all cell types, consistent with the largest AS event number in LSECs (Supplementary Fig. S4B and 1B), suggesting the functional diversity of LSECs and potential interplay with other cell types in various aspects.

4.5. Specific Tak1 isoform regulates hepatic cell-type-specific functions

Growing evidence suggested that protein isoforms generated from AS can perform distinct or even opposite functions., We attempted to reveal the functions of cell-specific isoforms of Tak1 in distinct cell types. Tak1 has four main variants generated by AS (Supplementary Fig. S5A), depending on the spliced-in or -out of exon 12 and exon 16. We provided a detailed splicing map of exon 12 of Tak1 (TGF-activated kinase 1) gene. Rbfox2, which highly expressed in LSECs, was reported to splice the Tak1 pre-mRNA to the short variant Tak1-A, consistent with the distribution of Tak1 variants in HCs and LSECs. Besides, we used RBPmap to predict the splicing factors that considered to be involved in the formation of Tak1’s splicing pattern, and several candidates with distinct expression patterns were selected, as shown in Supplementary Fig. S5B. We considered that splicing factors with different expression patterns were more likely to contribute in the formation of the cell-type-specific splicing patterns, while in the exonic regions, the predicted splicing factors did not exhibit significant expression differences between HCs and NPCs, thus, the intronic area might be more important in the AS regulation of Tak1. Several splicing factors were proved to be involved in this regulation, but single splicing factor did not play a decisive role (Supplementary Fig. S5C), which indicated that the specific splicing pattern of Tak1 resulted from the possible cooperation and competition of multiple splicing regulators in mouse liver cells.

We collected 561 public RNA-Seq datasets from the Sequence Read Archive database with RPBmap tools to obtain a comprehensive expression profile of Tak1 variants.Tak1-B was highly expressed in tissues with strong proliferation and metabolism properties, whereas Tak1-A was usually existed in tissues with immune functions, such as the spleen and several immune cells (Supplementary Fig. S5D). Coordinately, a multiple instance learning based analysis revealed that Tak1-B was more concentrated in some basic metabolic processes, while Tak1-A was primarily responsible for immune functions (Fig. 6A), consistent with the distribution of Tak1 variants in HCs and NPCs.

Figure 6

The functional differences between Tak1 isoforms. (A) Two variants of Tak1 were responsible for diverse biological processes. (B) The differences in protein structure between Tak1-A and Tak1-B isoforms. (C) SSOs converted Tak1-B to Tak1-A in AML12 cells. (D) Tak1-A bound to Tab2 (not Tab1) more easily than Tak1-B did. The ratio of Tak1/Tab1 or Tak1/Tab2 was calculated from image intensity. (E) The MAPK and NF-κB pathways were activated when Tak1-B converted to Tak1-A. (F) Several lipids and cholesterol synthesis related genes were up-regulated after the conversion. Significance was calculated using the t-test. *P<0.05.

We downloaded the PDB files of the two isoforms to predict the protein structure differences between Tak1-A and Tak1-B by searching Protein Data Bank (PDB) archive. Although the serine-threonine/tyrosine protein kinase catalytic domain of the two isoforms were unaffected by AS (from Pfam database), they had differences in spatial structures around the peptides that translated from exon 11 to exon 13 (Fig. 6B). The peptide segment altered by the spliced exon 12 fell in the protein-binding region (marked in Fig. 6B), suggesting that gain or loss of Tak1 isoform function might be caused by the affected interaction protein partner.

We then used SSOs to explore the exact functional differences between Tak1 isoforms. We converted Tak1-B to Tak1-A in AML12 cells, in which Tak1-B was the main variant, just like in the HCs (Fig. 6C). As we analyzed previously, the CEEs could change the PPIs by affecting the protein binding domains, here, Co-IP showed that Tak1-A bound to Tab2 (but not to Tab1) more easily than Tak1-B did (Fig. 6D). It was reported that the binding to Tab1 and Tab2 was necessary for the activation of Tak1 and the downstream pathways. Our results showed that as the conversion of Tak1-B to Tak1-A, the downstream JNK, p38 and NF-κB pathways were significantly activated (Fig. 6E). Activated JNK and p38 pathways had important roles in cellular inflammatory response, which consistent with the distribution and the functionality of Tak1-A in NPCs.

We also found that the conversion of Tak1-B to Tak1-A triggered a significant up-regulation of genes related with lipids and cholesterol synthesis in AML12 cells, whereas lipolysis-related genes (Lipc) were down-regulated (Fig. 6F). These results were also consistent with the potential metabolic regulation function of Tak1-B in HCs (Fig. 6A). In addition, Lkb1, the gene associated with maintenance of cell polarity in HCs, was significantly down-regulated (Fig. 6F), suggesting that Tak1-B also played important roles in maintaining cell polarity (Fig. 6A).

Thus, for the Tak1 gene, isoforms distributed in a cell-type-specific manner adequately fulfilled its functional plasticity. The full-length Tak1-B was the dominant isoform in HCs, regulating the cell-type specific functions such as lipid metabolism. While in NPCs, Rbfox2, as well as some other splicing factors, made Tak1-A the main isoform. Tak1-A bound to Tab2 more easily than Tak1-B did, activated the downstream JNK and p38 pathways. Thus, Tak1-A exhibited the immune process regulation properties in NPCs. These observations indicated that the regulation of CEEs by AS was closely associated with the cell-type-specific functions and had important roles in the formation of liver cell-type specificities.

5. Discussion

Although the cell-type-specific proteomic research in the mouse liver provided a better understanding of the identities of the individual cell type and how these cell types worked together, we still wanted to know how this complexity of proteomics was generated and how AS variants distributed and functioned across different cell types. In this article, we identified 13,637 AS events from all four cell types, with an average coverage of 81.5% of all known events. We revealed the first profile of cell-type-specific AS events defined in mouse liver, and this profile was associated with a broad range of functional processes of the liver.

AS variants had differential distribution patterns across various cell types. In HCs, the AS-regulated genes were highly enriched in functions that include ‘regulation of kinase cascade’, ‘phosphoprotein’ and ‘metabolism’, whereas the genes involved in ’regulation of immune effector process’ in KCs were most frequently regulated by AS. This finding indicated that the genes with distinct cell-type-specific splicing manners were relevant to the core biological functions of respective cell types.

Increasing studies reported that some AS events had crucial impacts on multiple biological processes and developmental stages, especially in the nervous system., In our research, the genes coding the kinase and transcription factors were more likely to be regulated by AS, which might be the important mechanism that AS caused wide impacts on multiple biological processes. More widely, the CEEs were significantly enriched in regions that involved in PPI regulation, which might also be the approach that AS affecting isoform functions.

AS is regulated by complicated interplays of cis-elements and trans-acting factors (splicing factors). The splicing factors recognize the ‘binding motifs’ to promote or suppress the splicing of the exon. According to Professor Phillip A. Sharp, some splicing factors are required for dominating the specialized splicing patterns and have critical roles in maintaining cellular homeostasis, namely the ‘master splicing factors’. Here, we found several relatively high-expressed splicing factors in each cell type, like Esrp2 in HCs and Rbfox2 in LSECs. These splicing factors dominated the formation of the cell-type-specific splicing patterns (like the splicing pattern of Tak1 in HCs and NPCs), demonstrating that the terminal of an alternative exon was depended on the splicing factors, cell types and the motif locations, revealing a precise regulation on the pre-mRNA level in specific cellular contexts. Besides, the functional analysis also showed that they regulated the specialized biological functions in respective cell types, which made them strong candidates for the ‘master splicing factors’ in mouse liver cells.

We were most interested in how a single AS event functioned diversely in distinct cell types. The conversion of Tak1-B to Tak1-A activated the downstream JNK and p38 pathways because of the enhanced Tak1-Tab2 interaction. The activated signal pathways had important roles in cellular inflammatory response, consistent with the distribution and the functionality of Tak1-A variant in NPCs. Besides, the conversion triggered a significant up-regulation of genes related to lipids and cholesterol synthesis, whereas lipolysis-related genes were down-regulated, coordinating with the functional predictions of Tak1-B. This example showed how AS regulated the isoform functions to shape the liver cell-type specificities.

In summary, our research provided the first profile of cell-type-specific transcriptome and splicing patterns in mouse liver and confirmed that the cell identities were largely defined by cell-type-specific AS events. These events altered isoform functions through remodeling PPI networks and modifying PTM sites and were controlled by a series of splicing factors. Analysis of cell-type-specific splicing enabled us to better demonstrating the functional landscape of hepatic cell types. Spliced genes played diverse roles in the corresponding cell types by differential exon usage even if they exhibited similar expression profiling at the gene level. This information improved our knowledge of additional regulation and appropriate composition of hepatic cell types.

Click here for additional data file.