Alternative splicing of IRF3 plays an important role in the development of hepatocarcinoma

ABSTRACT

Alternative splicing is a process causing mRNA translation to produce different proteins, and it is crucial for the development of tumours. In this study, we constructed a prognostic model related to alternative splicing events in hepatocarcinoma using bioinformatics analysis, including the alternative splicing of CSAD, AFMID, ZDHHC16, and IRF3. The model is an independent prognostic factor and can accurately predict a patient’s prognosis. IRF3 is a transcription factor related to the immune response. Its alternative splicing can affect the expression of various genes related to prognosis and plays an essential role in the tumour microenvironment. We also verified the expression of IRF3 exon skipping isoform in hepatocarcinoma at the mRNA level. In conclusion, we discovered that the alternative splicing of IRF3 is essential for the development of hepatocarcinoma. This study provides new insight into the development of treatments for hepatocarcinoma.

KEYWORDS:

• Alternative splicing
• exon skip
• hepatocellular carcinoma
• IRF3
• immune infiltration

Introduction

Hepatocellular carcinoma (HCC) is a prevalent malignant and invasive tumour. Early-stage tumours can be effectively treated with surgery, interventional chemotherapy, radiofrequency ablation or even liver transplantation. However, following the tumour has spread, these therapies are no longer effective [1]. Therefore, early detection and precision therapy are particularly important for HCC. Several valuable markers for early HCC identification have been discovered [2,3], while the effect remains to be considered. As a result, new strategies for early diagnosis and effective treatment should be further explored.

Alternative splicing (AS) is considered as a key factor in expanding the complexity of cell functions. AS of mRNA precursor molecules can produce a variety of mRNA and protein isoforms that differ in structure, location, and function [4]. The existence of AS provides a molecular basis for the diversity of biological phenotypes, and it also participates in and affects the progression of diseases, as evidenced by tumour research. Th1 cells, for example, can generate the replacement subtype of IRF1, decrease IFN γ secretion, and influence the anti-tumour effect [5]. AS may also alter the protein function and influence the incidence of HCC [6,7]. For instance, muscle blind splicing regulator 3 (MBNL3) can induce the transcript of lncRNA-PXN-AS1 containing exon 4, which can prevent the degradation of paxillin (PXN) and promote the occurrence of HCC [8]. AS has also become a novel target for tumour therapy [9–11], due to its critical function in tumour proliferation and invasion [12].

Approximately 90% of HCC cases are accompanied by chronic liver disease, which is caused by excessive drinking, viral hepatitis, and inflammation caused by non-alcoholic fatty liver disease [13]. The liver, being the most significant metabolic organ, is impacted by a variety of external factors. Chronic hepatitis is frequently caused by immune cell activation, the production of pro-inflammatory factors, the accumulation of exogenous toxins, and an excess of metabolic products. Chronic inflammation may alter immunological modulation and impact the development of liver cancer [14]. The immune microenvironment can induce alterations in alternative splicing types to ensure the survival advantage of cancer cells [15]. Therefore, the immune microenvironment plays a crucial role in the progression of HCC.

In this research, we structured a prognostic model for HCC based on AS events. The exon skip event of IRF3, an immune-related gene, was incorporated in the model Therefore, we thoroughly analysed the IRF3 AS events and explored their relationship with tumour microenvironment. This study provides new insight into the early diagnosis, prognosis assessment and targeted treatment of HCC.

Material and methods

Data acquisition and sorting

The transcriptome and clinical data were downloaded from LIHC patients in The Cancer Genome Atlas (TCGA) database. The TCGA SpliceSeq database was used to download the alternative splice events data. The events, which have a percent-splice-in (PSI) value in 75% of the samples, meet screening criteria. The exon structure of IRF3, the function of the corresponding region, and the isoforms of IRF3 were derived from the Ensembl database and the UniProt protein database. The IRF3 target gene set was from the GSE31477 data set which was selected in the Htftarget database [16].

Screening of survival-related exon skip events

ES (exon skip), AD (alternate donor site), AP (alternate promoter), AA (alternate acceptor site), AT (alternate terminator), ME (mutually exclusive exons) and RI (retained intron) are the seven types of AS events in TCGA SpliceSeq database. The AS events are named by combining the corresponding parent gene, ID number and splicing type. As an illustration, consider the meaning of ‘IRF3 | 50994 | ES,’ where IRF3 represents the corresponding gene 50,994 for the ID number, and ES for the splicing type. The KNN algorithm was used to fill the data gaps, and the events with a standard deviation of PSI < 0.05 were eliminated. ES events were combined with the survival time of clinical data, and a univariate Cox analysis was conducted to screen survival-related ES events. The UpSet plot diagram is depicted through the ‘upsetR’ package.

Construction of the prognostic model and verification of its independence

LASSO regression analysis was used to screen candidate ES events to prevent over fitting of the model. The risk score of each prognostic predictor was calculated by multivariate Cox regression analysis, which was also utilized to screen prognostic predictors. The formula was as follows: $\mathrm{r}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(\mathrm{E}\mathrm{S}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{t})={\sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}(\mathrm{E}\mathrm{S}\mathrm{i})\times \mathrm{P}\mathrm{S}\mathrm{I}(\mathrm{E}\mathrm{S}\mathrm{i})$. The ESi represents the ith prognostic ES events, PSI (ESi) is the PSI value of the ESi, and coef (ESi) is obtained by the regression coefficient of the multivariate Cox analysis, which represents the contribution of ESi. This method served as the foundation for our prognostic. Kaplan-Meier survival curves were analysed with the ‘survival’ R package. The ROC curve was created with the ‘timeROC’ R package. Cox regression analyses and stratified survival analyses were performed to verify the independence of this model.

Differentially expressed genes were screened for functional enrichment and prognostic analysis

Two groups of patients were established based on IRF3 expression. 25 patients were chosen as the high ES group and 26 patients as the low ES group based on the proportion of IRF3 ES events. Differentially expressed genes were screened with a cut-off of |logFC|>1 and FDR < 0.05 through the ‘limma’ R package. If logFC > 1, the gene was upregulated in the high ES group, and if logFC < −1, the gene was downregulated. The functional enrichment analysis was performed on the differentially expressed genes through the Metascape website (http://metascape.org). With the use of the GEPIA database and the Kaplan-Meier plotter database, the expression and prognosis of the intersection genes between differentially expressed genes and IRF3 target genes were examined.

Correlation among IRF3 expression, ES events and immune checkpoint genes

The Timer2.0 database was used to evaluate the correlation between IRF3 and tumour immune microenvironment. The correlation coefficient among IRF3 expression, IRF3 ES event proportion and immune checkpoints genes expression were calculated by ‘ggcorrplot’ R package.

cDNA microarray and real-time PCR analysis

The liver cancer cDNA microarray, containing 66 samples, was purchased from Shanghai Outdo Biotech co., Ltd (Cat. No.: MecDNA-HLivH087Su02; Shanghai, China). The clinical information of these patients was displayed in supplementary table S1. We designed characteristic primers according to the exon ENSE0002224211 missing or not. The specific primers were as follows: IRF3-IN: forward 5’- CGACA ATCCC ACTCC CTTCC-3“ and reverse 5”-CGACC CCACC AGCCG CAG-3;“ IRF3-EX: forward 5”-AAGCG GGGAA GATCT GATTAC-3“ and reverse 5”-GAGAG TGGGT GGCTG TGGGA-3.’ According to the instructions, the ChamQ Universal SYBR qPCR Master Mix (Vazyme, Q711–02) and Roche LightCycler 480-II instrument (Roche) were used for quantitative polymerase chain reaction (qPCR). In all of the samples, we measured the mRNA expression levels of IRF3-IN and EX isoforms relative to ACTIN. The relative mRNA expression was determined using the cycle threshold (CT) formula: $2\text{^}-\mathrm{\Delta }\mathrm{C}\mathrm{T}$, where $\mathrm{\Delta }\mathrm{C}\mathrm{T}=\left[\mathrm{C}\mathrm{T}\left(\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{t}\right)-\mathrm{C}\mathrm{T}\left(\mathrm{A}\mathrm{C}\mathrm{T}\mathrm{I}\mathrm{N}\right)\right]$. The formula for calculating the proportion of IRF3-EX isoform in IRF3 expression was calculated by $\mathrm{I}\mathrm{R}\mathrm{F}3-\mathrm{E}\mathrm{X}/(\mathrm{I}\mathrm{R}\mathrm{F}3-\mathrm{E}\mathrm{X}+\mathrm{I}\mathrm{R}\mathrm{F}3-\mathrm{I}\mathrm{N})$.

Statistics analysis

The differences between the two groups were evaluated using the Wilcoxon test. The period from the date of diagnosis to the date of death is known as overall survival. The survival curve was drawn with the Kaplan-Meier log-rank test. The relations were evaluated with the Pearson correlation test. Univariate and multivariate analyses were performed by the Cox regression model. The prognostic value of the prognostic model was evaluated by ROC curve analysis. The p < 0.05 indicates that the difference is statistically significant. All statistical calculations were performed using R software (version 4.0.2).

Results

Exon skip events in HCC are the main component of alternative splicing

We downloaded TCGA-LIHC AS events from the TCGA SpliceSeq and visualized them with an UpSet plot (Fig. S1A). The AS data of 352 patients were retained after being combined with the clinical data and eliminating patient information with seriously incomplete PSI values. Using a univariate Cox analysis, we selected AS events related to prognosis and drew an Upset plot (Figure 1a). The findings indicate that ES events accounted for the central part.

Figure 1. Construction of ES event prediction model.

a. The upset plot of AS events associated with prognosis.

b. Heatmap of ES events in the model. We ranked patients based on risk scores and assigned gradient colours to PSI values based on the range of PSI values for each ES event. Red represents high PSI values, and green represents low PSI values.

c. Distribution of survival status and risk score of HCC patients. In the survival status part, each point represents a patient. Red represents death and green represents alive. The risk score of patients was increased from left to right.

d. Kaplan-Meier curves of overall survival of patients. The patients in the low-risk group had a better prognosis than those in the high-risk group. The result was statistically significant (p < 0.05).

Constructing a prognostic model with ES events

Four ES events were included in the model to assess the prognosis using LASSO regression analysis and multivariate Cox analysis (Table 1). PSI value represents the proportion of the exon retained. Therefore, the exon retention of ZDHHC16 is a risk factor, whereas the exon retention of CSAD, AFMID and IRF3 are protective factors for the prognosis of HCC. We calculated the risk score of each patient and divided them into two groups according to the median value. The PSI value of ZDHHC16 in the high-risk group was higher, while the PSI values of the other three genes were the opposite (Figure 1b). The survival time was decreased with increasing scores (Figure 1c). There was a substantial difference in survival between the two groups, with a poor prognosis in the high-risk group (Figure 1d).

Table 1. The prognostic model related to exon skip events.

ES events	coef^1	HR^2	95%CI^3	pvalue
CSAD|21952|ES	−1.796	0.166	0.037–0.740	0.019
ZDHHC16|12708|ES	2.524	12.480	3.580–43.499	P < 0.001
AFMID|94694|ES	−1.595	0.203	0.061–0.670	0.009
IRF3|50994|ES	−2.438	0.087	0.019–0.402	0.002

Independent validation of the model

We subsequently examined the ROC curve for the purpose to assess the effectiveness of this prognostic model. The AUCs of 1-, 2- and 3-year survival were 0.779, 0.762, and 0.753, respectively, indicating that the model was effective at predicting prognosis (Figure 2a). In addition, univariate and multivariate Cox regression analyses explained that this model might act as an independent prognostic indicator in HCC (Figure 2b and 2c). The prognosis is often impacted by stage and age. Therefore, we stratified the patients according to age and stage and verified the model’s predictive accuracy. The prognosis in the high-risk group was worse than in the low-risk group (Figure 2d and 2e). Then, we investigated into the association between the model and clinical features. The risk score was gradually increased with increasing tumour grade and pathological stage, (Figures 3a and 3b). This model has the greatest accuracy for prediction compared with the ROC curve of other clinical features (Figure 3c).

Figure 2. The prediction ability and independent verification of the model.

a. ROC curve analysis of the model over time. The AUC was calculated for ROC curves, and sensitivity and specificity were calculated to assess score performance. This model showed good prediction ability.

b. Univariate Cox analysis. Analyzed the correlation between relevant factors and prognosis.

c. Multivariate Cox analysis. Comprehensive analysis of the correlation between relevant factors and prognosis.

d. Stratified patients by 65 years old, and compared the prognosis of patients with high and low ages according to risk assessment score. The results were statistically significant (p < 0.05).

e. Stratified patients by stage, and compared the prognosis of patients with early and late stages according to risk assessment score. The results were statistically significant (p < 0.05).

Figure 3. Correlation of model with clinical features.

a. The boxplot showed the correlation between risk scores and tumour grades.

b. The boxplot showed the correlation between risk scores and tumour stages.

c. Comparison of ROC curves of different clinical features and ES event model. ES event model has the best prediction ability.

d. The expression of IRF3 in HCC and liver tissues.

e. The boxplot showed the correlation between IRF3 expression and tumour grades.

f. Kaplan-Meier curves of patients between high and low IRF3 expression groups. The patients with low IRF3 expression had a better prognosis. The result was statistically significant (p < 0.05).

The function and expression of genes in the prognostic model

Three of the four genes studied in this model are CSAD, which is involved in the metabolism of synthesizing L-cysteine into hypotaurine, ZDHHC16, which is engaged in protein palmitoylation, and AFMID, which is involved in the subpathway of converting L-tryptophan into L-kynurenine. These three genes are all related to the protein metabolism pathway, whereas IRF3 participates in innate immune response and induces IFN expression. In addition, IRF3 is a transcription factor that can activate or repress protein expression. When we compared the expression of these genes between tumour and normal tissues, we found that only IRF3 was upregulated in tumours while the expression of other genes remained the same according to TCGA data (Figure 3d). The prognosis for patients with high IRF3 expression was poor (Figure 3e). Moreover, the expression of IRF3 was related to tumour grade (Figure 3f), but not stage (Fig. S1B). We speculate that the increased proportion of IRF3 ES events affects the function of IRF3.

The type of ES events in IRF3

Seven ES events related to IRF3 were discovered when we reviewed the ES events associated with prognosis from the previous analysis (Table 2). We revealed that the proportion of ES events arose as IRF3 expression increased, and some ES events were strongly correlated with one another by analysing the correlation between IRF3 expression and the proportion of ES events (Figure 4a). We validated the nucleotide and amino acid sequences of these events in the TCGA SpliceSeq database. The amino acid sequence mainly involved three exons: ENSE00003686782, ENSE00003656814, and ENSE00002224211 by searching the Ensembl database (Figure 4b). Two IRF3 isoforms, searching the UniProt protein database, were identified as meeting the ES conditions (Figure 4c). Both isoforms lack exon ENSE00002224211, which binds to HERC5. It is hypothesized that the deletion of this exon may have an impact on ISG15 binding and protein function since it contains ISG15 binding sites at both ends. The exons of ENSE00003686782 and ENSE00003656814 form the DNA binding region, and deletion of the region may disrupt IRF3 regulatory function.

Figure 4. The structure and isoforms of IRF3.

a. Association analysis between IRF3 expression and IRF3 ES event proportion (1-PSI).

b. Structure composition of IRF3

a. The arrangement of IRF3 exons displayed in TCGA SpliceSeq database. The thin exon represents the untranslated region (UTR), and the thick exon represents the coding region. The splice sequence exon number is obtained based on the TCGA SpliceSeq database.

b. The IRF3 encoding exons displayed in the Ensembl database. The number above represents the starting and ending nucleotide positions of the exon, while the number below represents the name of the exon in the Ensembl database.

c. The protein structure of IRF3 displayed in the UniProt database. The orange region is the DNA binding region, the blue region is the HERC5 binding region, and the red line is the ISG15 binding site. The number above represents the starting and ending nucleotide positions of the region.

c. Amino acid sequences of IRF3 and its two isoforms in the UniProt protein database.

Table 2. The IRF3 ES events associated with prognosis.

ES events	Skip exons^4	From	To	HR	95%CI	pvalue
IRF3|50995|ES	5.1:5.2	4	6.2	0.092	0.037–0.226	<0.001
IRF3|50994|ES	5.1:6.1	4	6.2	0.036	0.009–0.154	<0.001
IRF3|51007|ES	1.5:1.6:2	1.1	3	0.265	0.131–0.535	<0.001
IRF3|50998|ES	1.5:1.6:2	1.4	3	0.092	0.023–0.374	0.001
IRF3|51010|ES	1.4:1.5:2	1.1	3	0.222	0.087–0.568	0.002
IRF3|51012|ES	1.2:1.3:1.4:2	1.1	3	0.119	0.027–0.524	0.005
IRF3|51026|ES	1.5	1.1	2	0.151	0.032–0.707	0.016

The proportion of ES events related to IRF3 affects gene expression

We further explored whether the proportion of ES events related to IRF3 would affect gene expression. Based on the correlation between IRF3 and its ES events, the patients with PSI values of IRF3|50995|ES, IRF3|50994|ES, IRF3|51007|ES, IRF3|51010|ES, and IRF3|51012|ES are 100% were considered as the low ES group (26 patients). Due to both two IRF3 isoforms missing the exon ENSE0002224211, the patients with PSI values of IRF3|50995|ES and IRF3|50994|ES less than 75% were regarded as the high ES group (25 patients). The analysis of differentially expressed genes showed that 255 genes were upregulated and 175 genes were downregulated in the high ES group (Figure 5a and Fig. S1C). And the overall survival indicated that the high ES group had a worse prognosis (Figure 5b). These genes were connected to the cell cycle, carbon metabolism, apoptosis, and DNA replication according to GO analysis and KEGG analysis (Figure 5c). Then we searched the target genes of IRF3 in HepG2 cells through the Htftarget database. Seventeen genes that intersect in two sets were examined (Figure 5d). Additionally, the expression and prognosis of these 17 genes were queried through the GEPIA and Kaplan-Meier plotter databases (Figure 6a, b). Surprisingly, 15 genes upregulated in the high ES group were highly expressed in tumour tissues, and the majority of them were linked to poor prognosis. The two genes upregulated in the low ES group were highly expressed in normal tissues and correlated with a better prognosis. This suggests that the proportion of IRF3 ES events may significantly impact the gene expression of HCC.

Figure 5. Screening and grouping of patients based on the proportion of IRF3 ES events.

a. Screening differentially expressed genes between the high ES group and the low ES group. |LogFC|>1 and FDR < 0.05 were used as the criteria. Red indicates upregulation of gene relative expression, while blue indicates downregulation. The top 20 genes are shown in the figure.

b. Kaplan-Meier curves of patients in high ES and low ES groups. The patients with a low proportion of IRF3 ES events had a better prognosis. The result was statistically significant (p < 0.05).

c. Functional enrichment analysis of upregulated genes in the high ES group.

d. Venn diagrams of differentially expressed genes and IRF3 target genes. The blue circle represents the differentially expressed gene in the high ES group, and the red circle represents the IRF3 target genes. The genes in the middle intersection area are shown on the right. IRF3 target genes set was from the GSE31477 data set selected in the Htftarget database.

Figure 6. Expression and prognosis of differentially expressed target genes in HCC.

a. The expression of target genes between high and low ES groups. The left represents liver tumour tissue, and the right represents normal liver tissue. Data was obtained through the GEPIA database.

b. The prognosis of target genes in the cancer RNA seq database of Kaplan Meier plotter. The genes shown in the figure were statistically significant.

Correlation of IRF3 ES events with tumor immune characterization

During tumorigenesis, DNA damage can promote the activation of an inflammatory response by activating the cGAS-STING-IRF3 pathway [17]. Through the Timer2 database, we discovered that IRF3 was positively correlated with various immune cells, including inhibitory immune cells such as M2 and Treg cells (Figure 7a). In addition, we observed a connection between IRF3 expression and the expression of multiple immune checkpoint genes (Figure 7b and Fig. S2). We screened the immune checkpoint genes with a correlation coefficient above 0.2 and found that these gene expressions were also associated with the proportion of ES events. Moreover, DNA binding region deletion-related ES events were relatively more strongly correlated (Figure 7c).

Figure 7. Correlation of IRF3 and its ES events with tumour immune characterization.

a. The correlation between IRF3 expression and immune cell infiltration in the Timer2.0 database, including CD8⁺T cells, M1 macrophage, M2 macrophage, and Treg cells.

b. The boxplot shows the correlation between immune checkpoint gene expression and IRF3 expression level.

c. Association analysis between the proportion of IRF3 ES events and immune checkpoint gene expression.

Detection of IRF3 exon skipping by cDNA microarray

We named exon ENSE0002224211 skipping isoform as IRF3-EX (represents two splicing variants of IRF-3 in Figure 4c, Q14653–2 and Q14653–3), and exon ENSE0002224211 included isoform as IRF3-IN (represents the normal form of IRF3, Q14653–1). A cDNA microarray containing the clinical data of 66 liver cancer patients was used to identify the relative expression levels of the two IRF3 isoforms (Figure 8a, Fig. S3A). The IRF3-EX was prevalent in 66 patients, and its relative expression level was significantly positively correlated with the total relative expression level of IRF3 (Fig. S3B). Figure 8b depicted the percentage of IRF3-EX, and Figure 8c showed a positive relationship between it and the overall relative expression level of IRF3. A univariate Cox analysis displayed that the proportion of IRF3-EX was a risk factor for prognosis (Figure 8d). Although there was no statistically significant difference in the proportion of IRF3-EX below 10% (20 patients) and beyond 15% (14 patients), the overall survival of patients with high IRF3-EX proportion was inferior (Fig. S3C).

Figure 8. Detection of IRF3 exon skipping by cDNA microarray.

a. IRF3-EX refers to the transcriptional isoforms with missing exon ENSE0002224211.

b. The relative expression level of the IRF3-EX isoform in 66 liver cancer patients.

c. The proportion of IRF3-EX in 66 liver cancer patients. It was calculated by $\mathrm{I}\mathrm{R}\mathrm{F}3-\mathrm{E}\mathrm{X}/(\mathrm{I}\mathrm{R}\mathrm{F}3-\mathrm{E}\mathrm{X}+\mathrm{I}\mathrm{R}\mathrm{F}3-\mathrm{I}\mathrm{N})$.

d. Univariate Cox analysis. Analyzed the correlation between the proportion of IRF3-EX and prognosis.

Discussion

HCC is one of the most widespread clinical malignancies, and it is impacted by the genes, transcription, and epigenetic modification of tumour cells. Its high heterogeneity leads to a poor prognosis. Early detection, diagnosis, and treatment have been identified as therapeutic principles for treating liver cancer. From the standpoint of epigenetics, there are currently numerous prognostic models, such as the lncRNA [18] and N6-methyladenosine (m6A) modification [19] models. These models are crucial for the early diagnosis and prognosis evaluation of HCC.

In contrast to epigenetics, AS is a process following gene transcription. It is regulated by splice factors, which produce a diversity of proteins. Mutations in some splice factors lead to changes in splicing patterns and further induce the occurrence and development of HCC [20]. Poly(rC) binding protein 1 can regulate the AS of various cancer-related genes in HCC [21], and ATP-dependent RNA helicase MTR4 can drive liver cancer metabolism by controlling the AS of essential glycolytic genes [22]. Therefore, we analysed the AS events in HCC and constructed a prognostic model based on ES events. The model shows excellent prediction capacity.

Changes in the AS pattern can affect the function and even guide treatment. The lncRNA CRNDE can reduce the chemoresistance of gastric cancer by inducing SRSF6 to regulate the skipping of exon 14 of PICALM [23]. Targeted skipping of exon 17 of NF1 can treat type I neurofibromatosis [24]. Therefore, we investigated how these four genes functioned in the model. CSAD is considered the rate-limiting enzyme for taurine synthesis. In HepG2 cells, taurine can inhibit the expression of key glycolytic enzymes to suppress proliferation [25], and enhance the effect of tumour chemotherapy drugs [26]. ZDHHC16 regulates the palmitoylation of proteins, which is essential to signal transduction of cancer cells [27]. AFMID is involved in the metabolism of kynurenine, which can promote the invasion, metastasis, and drug resistance of tumour cells [28]. These three genes participate in protein metabolism. IRF3 is a transcription factor and plays a central role in innate immunity [29,30], and has been linked to the anti-cancer immune response [31]. Therefore, we focused on the role of IRF3 ES events in HCC.

The sequence of amino acids constitutes the primary structure of the protein, and the subsequent interactions between these amino acids form the protein domains that can determine the function of the protein. IRF3 has a DNA binding region, a disordered region, a mediating interaction with ZDHHC11 region, an interaction with HERC5 region and a nuclear export signal region according to the UniProt. The skipping of exons will affect the composition of domains, and may potentially result in proteins losing their corresponding functions. There are seven types of IRF3 ES events related to the prognosis of HCC, involving the deletion of three exons. The exons ENSE00003686782 and ENSE00003656814 form the DNA binding region, which can bind to the interferon-stimulated response element (ISRE) in the target promoter region and regulate the expression of related genes [32]; Exon ENSE00002224211 can bind to the E3 ligase HERC5, in which ISG15 binding regions exist at both ends. ISG15 can result in ISGlyation of the target protein, which can influence the activity and localization of the target protein, whereas HERC5 can catalyse the coupling of ISG15 and IRF3 and induce the continuous activation of IRF3 [33]. The deletion of these exons may consequently result in a reduction in the transcriptional activity of IRF3 and affect the expression of downstream target proteins. Currently, it has confirmed that the IRF-3 alternatively spliced isoform can regulate the function of IRF3 [34]. Additionally, the splicing isoforms of IRF3 have an obvious cancer-specific pattern in HCC, and can compete with the standard IRF3 protein to bind to the target promoter, inhibiting the transcriptional activation ability [35].

The pathway enrichment of differentially expressed genes revealed that the group with a high proportion of IRF3 ES events had higher proliferative and metabolic capacities. The effects of these genes on HCC have been confirmed by research. RCL1 and ALDH6A1 were downregulated in the high ES group. RCL1 can regulate the cell cycle and restrain the growth and metastasis of HCC cells [36]. ALDH6A1 is a mitochondrial protein that participates in mitochondrial respiration. Hepatic neoplastic transformation suppresses the expression of ALDH6A1, which influences mitochondrial function [37]. DDX39 is one of the upregulated genes that promote the development of HCC by activating the Wnt/β-Catenin pathway [38]. ASF1B, CDKN3, and CENPA can enhance cell proliferation and affect cancer development [39–41]. RRM2 has the ability to control nucleotide metabolism and promote tumour proliferation [42]. SPC24 and TPX2 can affect the PI3K/AKT pathway to support in tumour progression [43,44]. STIP1 can promote HCC development through the PI3K-Akt-dependent anti-apoptotic pathway [45]. TCF19 can control the progression of HCC by binding to histone H3K4me3 [46], interacting with p53 to mediate metabolic regulation [47]. TOP2A may trigger EMT to accelerate the metastasis of HCC by activating the p-ERK1/2/p-Smad2/Snail pathway [48]. These results suggest that the ES events of IRF3 May affect the development of HCC.

The immunological effects of IRF3 have been the primary subject of existing research. IRF3 regulates the transcription of type I IFN genes and IFN-stimulated genes (ISG) by combining the ISRE in their promoters [49,50]. We discovered that the expression of IRF3 was positively interrelated with the infiltration of cells that inhibited tumour immunity, such as Treg cells and M2 macrophages. Moreover, the IRF3 ES events were also interrelated with the expression of immune checkpoint genes. Therefore, in addition to interferon related pathways, IRF3 May affect the tumour immune response in other ways.

In addition, we found that IRF3-EX subtypes influenced the prognosis of patients, and the patients with a proportion over 25% had a worse prognosis than the patients with a proportion is 0 in TCGA-LIHC patients. Furthermore, the IRF3-EX was prevalent in the cDNA microarray (MecDNA-HLivH087Su02), and the patients with a proportion over 15% indicated a trend of poorer prognosis than the patients with a proportion below 10%, although there is no statistical significance. Therefore, a larger sample size is required for the analysis of the relationship between IRF3-EX proportion and prognosis in a larger sample size.

In this study, we constructed a model with ES events in HCC, which can effectively evaluate the prognosis. We conducted an in-depth analysis of the IRF3 ES events, providing broader insights into the progress and treatment of HCC. Nevertheless, there are some deficiencies in this application as well. For instance, the samples we obtained were the cDNA microarray from liver cancer patients, but not tissue. Therefore, we are unable to verify the protein expression level of IRF3 splicing isoforms. Besides, we explore the gene expression difference based on the transcriptome data of TCGA, and speculate the effect of IRF3 isoforms on hepatocarcinoma, without further exploration through experiments. Therefore, it is necessary to verify in clinical samples or cell lines. Additionally, calculating the PSI value in clinical practice may be challenging and the reasons for the generation of IRF3 splicing isoforms are also worth further research.

Ethical approval

The cDNA microarray experiments were conducted according to the guidelines of the Declaration of Helsinki, and approved by the Medical Ethics Committee of the Shanghai Outdo Biotech Company (ID: YBM-05-02).

Author contributions

Z.L. wrote the manuscript. Z.L. performed the bioinformatics analysis. W.G. analysed the results of cDNA microarray, Z.L. and Z.D. conducted the statistical analysis and data interpretation. W.G., T.L. and Y.Z. performed the data review. D.X. and C.H. participated in the study conception and supervision, data analysis, and manuscript editing. All authors read and approved the final manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2023.2276371

Additional information

Funding

This work was supported by National Natural Science Foundation of China (No. 82100675), National Natural Science Foundation of China (No. 82170654) and Excellent Youth Foundation of the first Affiliated Hospital of Harbin Medical University (2021Y12).

Notes

1. Coefficient.

2. Hazard Ratio.

3. Confidence interval.

4. The exon number is obtained based on the TCGA SpliceSeq database.