Alternative splicing signatures of congenital heart disease and induced pluripotent stem cell-derived cardiomyocytes from congenital heart disease patients.

Congenital heart disease (CHD) is the most serious congenital defect in newborns with higher mortality. Alternative splicing (AS) plays an essential role in numerous heart diseases. However, our understanding of the link between mRNA splicing and CHD in humans is limited. Here, we try to investigate the genome-wide AS events in CHD using bioinformatics methods. We collected available RNA-seq datasets of CHD-induced pluripotent stem cell-cardiomyocytes (iPSC-CMs) (including single ventricle disease [SVD] and tetralogy of Fallot [TOF]) and non-CHD from the Gene Expression Omnibus database. Then, we unprecedentedly performed AS profiles in CHD-iPSC-CMs and non-CHD-iPSC-CMs. The rMAPS was used to generate RNA-maps for the analysis of RNA-binding proteins' (RBPs) binding sites. We used StringTie to identify and quantify the transcripts from aligned RNA-Seq reads. A quantification matrix was generated with respect to different groups by extracting the transcripts per million values from StringTie outputs. Then, this matrix was used for correlation analysis between the expression level of RBP and AS level. Finally, we validated our AS results using RNA-seq data from CHD and non-CHD patient tissue samples. We identified CHD-related AS events using CHD-iPSC-CMs and CHD samples from patients. The results showed that functional enrichment of abnormal AS in SVD and TOF was transcription factor-related. Using rMAPS, RNA-binding proteins which regulated these AS were also determined, and RBP-AS regulatory network was constructed. Overall, we identified abnormal AS in CHD-iPSC-CMs and CHD samples from patients. We predicted AS regulators in SVD and TOF, respectively. At last, we concluded that AS played a key role in the pathogenesis of CHD.

1. Introduction

Congenital heart disease (CHD) is the most serious congenital defect in newborns with higher mortality.[ Our understanding of the genetic regulation of CHD has been improved recently, but despite this effort, major knowledge gaps of genetic regulation are still the key barriers in the effective prevention and treatment of CHD. As an important genetic regulation mechanism, alternative splicing (AS) plays an important role in numerous diseases including coronavirus disease 2019, cancers, and heart diseases.[ For example, AS could act as potential therapeutic targets and biomarkers in coronavirus disease 2019 patients.[ Abnormal splicing is regulating the cancer-specific splicing patterns, which may be a source of neoepitopes in cancer.[ Besides, mutations in splice sites of NR2F2 and GATA4 which are known transcription factors were reported to be associated with CHD.[ Furthermore, 12 small Cajal body-associated RNAs may influence splicing during heart development and regulate heart development,[ suggesting that AS is associated with CHD. To further understand the AS mechanism in CHD, whole genome-wide profiling and analysis of RNA splicing in CHDs need to be explored.

Single ventricle defect is a type of heart defect that occurs in 1 of the 2 pumping chambers in the heart.[ Tetralogy of Fallot (TOF) is the most common cyanotic heart condition in children that have survived beyond the neonatal period.[ As described by Joseph et al,[ they generated induced pluripotent stem cell (iPSC) lines from 5 patients with single ventricle disease (SVD), 5 patients with TOF, and 5 healthy individuals (non-CHD). Then, they differentiated all iPSC lines into cardiomyocytes (iPSC-CMs) and used highly enriched iPSC-CMs for RNA sequencing.[ Those iPSC-CMs from CHD patients are increasingly used in the research of cardiovascular disease modeling, drug discovery, cardiotoxicity screening, and cardiac development[; these are important resources for investigating AS mechanism in CHD. We also collected RNA-seq data from CHD patients to validate AS results.

Based on the above information, we aimed to reveal the overall dynamic changes of AS in iPSC-CMs and CHD patients. In this study, we try to identify all potential abnormal AS in CHD-iPSC-CMs and CHD patients using replicate multivariate analysis of transcript splicing (rMATS), a computational tool that infers differential AS events based on RNA-seq data. We aimed to investigate AS mechanism in CHD-iPSC-CMs and CHD patients.

2. Materials and Methods

2.1. Datasets employed and sequence alignment of RNA-seq reads of CHD and non-CHD

To obtain an understanding of the splicing alterations in CHD, we collected available RNA-seq datasets corresponding to the raw RNA sequence reads of CHD-iPSC-CMs and non-CHD-iPSC-CMs from Gene Expression Omnibus (GSE132401).[ To validate our AS results, independent RNA-Seq data from CHD patients and non-CHD groups were downloaded from the Gene Expression Omnibus (GSE36761). We selected RNA-Seq datasets of CHD when N ≥ 3, because small runs of at least 3 individuals from each condition for rMATS are acceptable.[ Briefly, we downloaded the datasets in FASTQ format using the SRA Toolkit. The quality of the aligned reads to a minimum quality score of 20 for each sample was ensured using FASTX-Toolkit. HISAT was used for aligning short reads from RNA sequencing onto the human reference genome.[ We used HISAT with default parameters against human reference genome hg38 annotation files.[ We obtained SAM files as outputs from HISAT, and SAM files were postprocessed using SAMtools to convert SAM to BAM format.

2.2. Transcript identification and quantification from the aligned RNA-seq datasets

To identify the expression value in CHD-iPSC-CMs and CHD tissue samples, StringTie was used for identification and quantification of transcripts from the aligned RNA-seq reads.[ Then, a quantification matrix was generated with respect to different groups by extracting the TPM (transcripts per million) values from StringTie outputs. Then, this matrix was used for downstream analysis.[

2.3. Analysis of differential AS

RNA-seq data from iPSC-CMs and CHD patients was used to detect differential AS events across conditions. We applied rMATS to identify differential AS events. We used BAM files, which were obtained from the alignment of RNA-seq datasets as discussed above, as input to rMATS.[ The software rMATS quantified AS event based on the difference in the level of inclusion of an exon which is defined as the splice index or percentage splicing index (PSI, ψ score) between 2 conditions. PSI value was based on short (S) and long (L) forms of all splicing events present using the equation as follows:

PSI represents the inclusion/exclusion of an exon for a transcript isoform considering all alternate possible isoforms.[ We ran rMATS for all pairs of CHD and non-CHD as follows: python rmats.py --b1 b1.txt --b2 b2.txt --gtf Homo_sapiens.GRCh38.101.gtf -t paired --nthread 8 --od output (b1 and b2 are bam files generated from SAMtools), generated a summary table with significantly different AS events (with false discover rate [FDR] < 0.1, P < .01). Functional enrichment analysis of genes belonging to these splicing events was performed using gene set enrichment analysis (GSEA). We performed GSEA using WEB-based Gene Set Analysis Toolkit[ with FDR < 0.1. The gene set used in enrichment operations was the gene ontology (GO) gene sets. The analytical parameters were set as follows: the t-statistic mean of the genes was computed in each GO terms using a permutation test with 1000 replications and gene sets containing <3 genes or >500 genes were excluded. Quantitative realtime polymerase chain reaction (qPCR) was also performed by splicing-specific primers using 2× Green PCR Mix using CHD-iPSC-CMs and non-CHD-iPSC-CMs. CHD/non-CHD-iPSC-CMs were cultured by the Kunming Institute of Zoology. The study was carried out in accordance with the recommendations of the Ethics Committee of the Second Affiliated Hospital of Kunming Medical University (no. PJ-2022-14) and was approved by said committee.

2.4. Identifying RNA-binding protein (RBP) regulating splicing alterations in CHD-iPSC-CMs and non-CHD-iPSC-CMs

We converted bam files to Wiggle files. Wiggle files were generated and normalized in terms of reads per million reads. Next, these files were converted into bigwig files using the wigToBigWig tool (http://hgdownload.cse.ucsc.edu/admin/exe/) and visualized in Integrative Genomics Viewer.[ The rMAPS (http://rmaps.cecsresearch.org/) is a web server that generates RNA-maps for the analysis of RBPs binding sites. We uploaded rMATS output and performed the analysis of binding sites around differential AS events for over 100 known RBPs.[

3. Results

3.1. Dynamic AS events in SVD-iPSC-CMs

We first performed an AS analysis in SVD by rMATS. As shown in Figure 1, there are 74 significant AS events in SVD (Fig. 1A) (FDR < 0.1, P < .01), including 37 skipped exon, 4 alternative 5′ splice sites (A5SS), 7 alternative 3′ splice sites (A3SS), 9 mutually exclusive exons (MXEs), and 17 retained introns (RI) (Table S1, Supplemental Digital Content 1, http://links.lww.com/MD/H51). We illustrated PSI value difference for each gene and each group (Fig. 1B) to show the AS in every significant splicing site.

Figure 1.

Significantly different splicing events between non-CHD-iPSC-CMs and SVD-iPSC-CMs groups. (A) Volcano plot summarizing the different splicing events between non-CHD-iPSC-CMs and SVD-iPSC-CMs (likelihood ratio test). False discover rate < 0.1 are colored. (B) Heat map of PSI value showing the 74 significantly different splicing events between non-CHD-iPSC-CMs and SVD-iPSC-CMs; PSI value difference for each gene and each group is illustrated. (C) GO term distribution of alternative splice events between non-CHD-iPSC-CMs and SVD-iPSC-CMs groups. (D) Network of interactions among different splicing events related genes, as illustrated by the GeneMANIA database. A set of known different splicing events related genes were provided as a query. Circles represent genes and the gray box represents enriched pathways. Connecting lines represent interactions between genes. Purple lines indicate co-expression, and light pink lines indicate experimentally determined physical interactions. Light yellow lines indicate predicted relationships, and light brown lines indicate shared protein domains. CHD = congenital heart disease, CMs = cardiomyocytes, GO = gene ontology, iPSC = induced pluripotent stem cell, PSI = percentage splicing index, SVD = single ventricle disease.

As shown in Figure 1C, leukocyte differentiation, regulation of cell cycle phase transition, negative regulation of cell cycle process, and regulation of immune system process were enriched in GO (Biology Progress). We also conducted a gene-gene functional interaction network from AS genes list using GeneMANIA, which was a real-time multiple association network integration algorithm. We found interferon regulatory factor played a key role in AS of SVD (Fig. 1D). The functional and pathway enrichment of AS in CHD-iPSC-CMs was similar with functional enrichment of transcriptomic profiles in CHD-iPSC-CMs (Fig. 1D). Because transcription factor-binding motifs including transforming growth factor-beta-induced factor and interferon regulatory factor were reported to be enriched in SVD-iPSC-CMs,[ AS may be involved in the regulation of transcription factors. Thus, AS might regulate gene expression by regulating transcription factors in CHD.

3.2. Dynamic AS events in TOF-iPSC-CMs

To further explore the relationship between AS and TOF, we performed an AS analysis in TOF. As shown in Figure 2A, there are only 23 significant AS events in TOF-iPSC-CMs (Fig. 2A) (FDR < 0.1, P < .01), including 15 skipped exon 1 A5SS, 1 A3SS, 2 MXEs, and 4 RI (Table S2, Supplemental Digital Content 2, http://links.lww.com/MD/H52). We also illustrated PSI value difference for each gene and each group (Fig. 2B).

Figure 2.

Significantly different splicing events between non-CHD-iPSC-CMs and TOF-iPSC-CMs group. (A) Volcano plot summarizing the different splicing events between non-CHD-iPSC-CMs and TOF-iPSC-CMs (likelihood ratio test). False discover rate < 0.1 are colored. (B) Heat map of PSI value showing the 23 significantly different splicing events between non-CHD-iPSC-CMs and TOF-iPSC-CMs; PSI value difference for each gene and each group is illustrated. (C) GO term distribution of alternative splice events between non-CHD-iPSC-CMs and TOF-iPSC-CMs group. (D) Network of interactions among different splicing events related genes, as illustrated by the GeneMANIA database. A set of known different splicing events related genes were provided as a query. Circles represent genes and the gray box represents enriched pathways. Connecting lines represent interactions between genes. Purple lines indicate co-expression, and light gray lines indicate InterPro predicted interactions. Light brown lines indicate shared protein domains, and gray lines indicate miRNA target predictions. CHD = congenital heart disease, CMs = cardiomyocytes, GO = gene ontology, iPSC = induced pluripotent stem cell, PSI = Percentage Splicing Index, TOF = tetralogy of Fallot.

As shown in Figure 2C, the enrichment of GSEA function in AS is mainly related to epigenetic modification, cell cycle, cell adhesion, and immune cell function. Besides, dephosphorylation, positive regulation of cell cycle, lymphocyte mediated immunity, and cell-substrate adhesion were enriched in the gene-gene functional interaction network from AS gene list (Fig. 2D). Obviously, minor changes in AS have been detected in TOF-iPSC-CMs compared with that in SVD-iPSC-CMs.

3.3. Examples of dynamic AS events in CHD-iPSC-CMs

Here, we also showed examples of AS in CHD-iPSC-CMs. Interferon regulatory transcription factor (IRF) 7 encodes interferon regulatory factor 7, which is a member of the IRF family. IRF7 is a novel negative regulator of pathological cardiac hypertrophy.[ We found a significant reduction in the level of the IRF7 isoforms of SVD-iPSC-CMs (P = 4 × 10−5) (Fig. 3A), suggesting the expression level of IRF7 is regulated by AS. Similar results were observed for another gene KANK1 (KN Motif and Ankyrin Repeat Domains 1) in TOF-iPSC-CMs (Fig. 3B) (P = .001). Then, we incorporated the qPCR data to verify the different transcript isoforms of IRF7 and KANK1 (Fig. 3C), and the qPCR results were consistent with the analysis results. KANK1 is affiliated with the protein-coding class, and diseases associated with KANK1 is inherited congenital disease.[ The mechanism of KANK1 and IRF7 regulating AS was worthy to be studied further.

Figure 3.

The AS examples results of IRF7 and KANK1 genes in Tracks and qPCR. (A and B) Tracks displaying representative AS examples of the read coverage on IRF7 and KANK1 genes. The abnormal AS events are highlighted with box. (C) The shorter isoforms of IRF7 and KANK1 are differentially expressed in non-CHD and CHD by qPCR analysis. AS = alternative splicing, CHD = congenital heart disease, IRF = interferon regulatory transcription factor, KANK1 = KN Motif and Ankyrin Repeat Domains 1, qPCR = quantitative realtime polymerase chain reaction.

3.4. RBP regulating dynamic AS events in CHD-iPSC-CMs

To further identify the role of RBPs in the regulation of the differential AS events of SVD and TOF, we analyzed known RBP motifs and binding sites around the differentially regulated alternatively spliced exons using rMAPS.[ The results depicted the spatial distribution of the FXR2 motif [AGT]GAC[AG][AG][AG] around differential exon skipping events and showed FXR2 motif [AGT]GAC[AG][AG][AG] enrich downstream of the silenced and enhanced exons in SVD (Fig. 4A). Besides, PCBP2 showed the same results. For TOF-iPSC-CMs, the result depicted the spatial distribution of the HuR and TIA1 around differential exon skipping events in TOF-iPSC-CMs (Fig. 4B), indicating HuR and TIA1 regulated AS in TOF-iPSC-CMs. The relationship between splicing factors and AS events were calculated by rMAPS, and RBP-AS regulation networks were constructed using Cytoscape.[ To validate the relationship calculated by rMAPS, the Spearman test was used to evaluate the correlations between splicing factor PCBP2 and AS events in IRF7 (Fig. 4C and D). Spearman correlation coefficient showed that the results of the theoretical calculation of RBP-AS network were reliable.

Figure 4.

Results of RNA-maps. (A) RNA-maps depicting the spatial distribution of RBP motifs in the vicinity of differential alternative splicing events of SVD-iPSC-CMs. The red line represents the enriched motif for enhanced exons, the blue line represents the enriched motif for silenced exons, and the black line represents the motif density for background (nonregulated) exons. Solid lines represent the peak quality Motif score (peak height) as scaled on the left. Dotted lines represent the negative log10 (P value) as scaled on the right. The green box indicates the cassette exon. (B) RNA-maps depicting the spatial distribution of RBP motifs in the vicinity of differential alternative splicing events of TOF-iPSC-CMs. (C) splicing network in CHD-iPSC-CMs constructed by Cytoscape. (D) Correlation between IRF7 expression and the PSI value of PCBP2. CMs = cardiomyocytes, iPSC = induced pluripotent stem cell, PSI = Percentage Splicing Index, RBP = RNA-binding proteins, SVD = single ventricle disease, TOF = tetralogy of Fallot.

3.5. Dynamic AS events in CHD patients

To validate the results of dynamic AS in CHD-iPSC-CMs, we analyzed RNA-seq data from CHD patients. We observed the 1176 significant AS events (905 genes) in CHD (Fig. 1A) (FDR < 0.1, P < .01), including 942 skipped exon, 56 A5SS, 92 A3SS, 37 MXEs, and 85 RI (Fig. 5A), suggesting the change of AS is associated with CHD patients. Interestingly, most of these AS events in CHD patients are downregulated (Fig. 5A), and AS events in CHD-iPSC-CMs are downregulated compared with non-CHD-iPSC-CMs (Figs. 1A and 2A). Moreover, there is significant overlapping among AS events in CHD patients and AS events in CHD-iPSC-CMs were calculated by testing whether the random genes had an equal fraction of AS events in CHD patients overlapping with the AS events in CHD-iPSC-CMs using χ2 tests (Fig. 5B), suggesting AS event results are credible. Furthermore, we also found a reduction in the levels of the IRF7 isoforms of CHD patients (P =.008) (Fig. 5C).

Figure 5.

Significantly different splicing events between non-CHD and SVD patients group. (A) Volcano plot summarizing the different splicing events between non-CHD and TOF (likelihood ratio test). Significant hits at an false discover rate = 0.1 are colored. (B) Venn diagram shows the overlapping between abnormal AS events in CHD-iPSC-CMs and CHD patients. (C) Tracks displaying representative AS examples of the read coverage on IRF7 genes. The abnormal AS events are highlighted with box. AS = alternative splicing, CHD = congenital heart disease, CMs = cardiomyocytes, iPSC = induced pluripotent stem cell, IRF = interferon regulatory transcription factor, SVD = single ventricle disease, TOF = tetralogy of Fallot.

4. Discussion

Using RNA-seq data from patient-specific iPSC-CMs, we innovatively demonstrated AS in CHD-iPSC-CMs. The PSI value in CHD-iPSC-CMs was much lower compared with that in non-CHD-iPSC-CMs, indicating that AS might play an important role in CHD-iPSC-CMs. Compared with TOF-iPSC-CMs, SVD-iPSC-CMs showed major changes in AS. It manifested that SVD-iPSC-CMs and TOF-iPSC-CMs had different biological mechanisms of AS. We proposed rMAPS to identify RBPs regulating AS in CHD-iPSC-CMs. We identified key RBPs, FXR2 and PCPB2, might cause the AS in CHD-iPSC-CMs.[ Moreover, we collected RNA-seq data to analyze AS in CHD patient tissue samples to validate our AS results in CHD-iPSC-CMs. Furthermore, RBP is the main molecular determinant of AS, and perturbations in RBP-AS network activity are related to disease development. We proposed rMAPS to identify RBP regulating AS in CHD. These RBP may cause AS in CHD.[ Additionally, we used qPCR from CHD-iPSC-CMs and non-CHD-iPSC-CMs to validate our results. The experiment results are consistent with our bioinformatics results. Moreover, we also found 1176 AS events in CHD patients and 74 significant AS events in SVD. It was reported that exon 8 of IRF7 translated to an autoinhibitory domain.[ This autoinhibitory domain is capable of silencing the active transactivation domains. Furthermore, induction of endogenous IFN gene expression required relief of IRF7 autoinhibition. Lacking the autoinhibitory induced abnormal interferon gene expression. Interferon plays a key role in the development of the autoimmune CHD.[ Therefore, exon 8 skipping of IRF7 may result in a loss of autoinhibitory domain and cause dysfunction of the immune system. Our study was the first genome-wide analysis that suggested AS was related to CHD.

Given the AS in generating protein diversity and its link to many human diseases, AS has become a key target of therapeutic intervention.[ For instance, AS profiling in identifying biomarkers for the prognosis of endometrial cancer provided comprehensive insights into the molecular mechanisms and therapeutic targets involved in endometrial cancer processes.[ AS is a potential therapeutic option for Pelizaeus-Merzbacher disease, spastic paraplegia 2, and hypomyelination of early myelinating structures.[ Even though it is widely accepted that AS plays a key role in heart development,[ for the lack of genome-wide AS analysis, fewer studies focused on the AS mechanism in CHD. Chloe et al[ found 12 small Cajal body-associated RNAs involved in AS during heart development and regulate heart development; however, AS analysis in genome-wide was not performed. Chen et al[ reported that the regulation of RNA splicing was an important player in transcriptome reprogramming during heart failure. With the development of genome-wide sequencing technologies for the past few years, it is now possible to identify genome-wide AS related to CHD. We expected that this study could provide a novel insight into the underlying mechanisms of RNA processing in CHD and play a good foundation for targeted therapy of CHD patients.

In addition, further researches need to be conducted in the future. Concerning AS, the wet experiment should be used to verify the results of AS calculated by bioinformatic methods.[ For instance, qPCR will be used to verify different transcript isoforms calculated through bioinformatics analysis. However, it was widely accepted that AS calculated by rMATS is reliable.[

5. Conclusion

In summary, our study firstly reported the AS and predicted AS regulators in CHD based on the data of RNA-seq from CHD-iPSC-CMs and CHD patients. We concluded that AS played a key role in the pathogenesis of SVD and TOF.

Author contributions

Conceptualization: Qianqian Su, Xiang Xu.

Data curation: Renchao Zou, Xiang Xu, Xiaoyong Liu.

Formal analysis: Renchao Zou, Xiang Xu.

Methodology: Qianqian Su.

Project administration: Qianqian Su, Xiang Xu, Xiaoyong Liu.

Software: Renchao Zou.

Supervision: Xiaoyong Liu.

Validation: Qianqian Su.

Writing – original draft: Xiang Xu.

Writing – review & editing: Qianqian Su, Xiang Xu.