In silico splicing analysis of the PMS2 gene: exploring alternative molecular mechanisms in PMS2-associated Lynch syndrome

Abstract

Lynch syndrome (LS) is one of the most common hereditary cancer syndrome in human populations, associated with germline variants in MLH1, MSH2/EPCAM, MSH6 and PMS2 genes. The advent of next generation sequencing has proven a significant impact in germline variant detection in the causative genes; however, a large proportion of patients with clinical criteria still receive uncertain or negative results. PMS2 is the least frequent reported gene, associated with up to 15% of LS cases with late-onset disease and low penetrance phenotype; however, the proportion of PMS2-LS cases is considered to be highly underestimated. In this context, our analysis aimed to improve the current diagnostic yield by focusing on missense and intronic PMS2 variants available in public clinical databases (ClinVar, LOVD). We performed an in silico assessment of the wild-type DNA sequence and the reported genetic variants, employing splicing bioinformatics tools known for their effectiveness in other genes. Splicing variants were predicted in silico and using GTEx short-read RNA expression data. Out of the 2384 missense variants discovered, 90% were classified with uncertain significance (VUS). 4.9% of missense variants were shown to have a potential splicing consequence (DS > 0.2) using SpliceAI. As described in the original publication, SpliceAI-visual was proven effective in annotation of short intronic variants (< 50 bp). Four short intronic variants were identified using SpliceAI-visual as potentially splicing disturbing, in spite of using a lower threshold (DS > 0.1). Exons 2, 3, 4, 5, 6, 7, 8, 11, 12 and 14 were consistently predicted in at least three out of eight software with weak canonical splice sites. Additionally, we noted that both Exonic Splicing Enhancers (ESEs) and Exonic Splicing Silencers (ESSs) contribute significantly to alternative splicing and exonic selection in PMS2 gene. Specifically, ESE motifs were consistently more abundant in highly expressed exons 5, 11 and 14, while ESS motifs played a fundamental role in exons 6, 7 and 10. Computational analysis performed in our study serves as a valuable filtering step for guiding further RNA experiments. Additional functional data is necessary to validate our findings.

Introduction

Lynch syndrome (LS, OMIM #120435, #609310, #614350, #614337, #613244) is one of the most prevalent cancer predisposition syndromes known in humans, with an estimated prevalence in the general population of 1 in 300. LS is characterized by pathogenic and likely pathogenic germline variants in DNA mismatch repair (MMR) genes – MLH1, MSH2/EPCAM, MSH6 and PMS2 – discovered by current molecular approaches [1,2,3,4,5,6]. Despite the latest advancements in high-throughput sequencing technologies, germline testing still continues to pose significant challenges in PMS2-associated Lynch syndrome (PMS2-LS) in current clinical practice. Pathogenic and likely pathogenic PMS2 variants are reported in literature in up to 15% of all individuals diagnosed with Lynch syndrome, a recently acknowledged underestimation [1,2,3,4,5,6]. Supporting observations for this assertion include the high reported frequency of PMS2 defects in individuals with Constitutional Mismatch Repair Deficiency (CMMRD) syndrome (OMIM #276300), a rare pediatric cancer syndrome [1, 3, 7]. In heterozygous state, PMS2 pathogenic variants are associated with Lynch-related cancers, mainly colorectal and endometrial cancers [8]. However, the low penetrance phenotype and late-onset disease associated with heterozygous PMS2 variants are contributing factors that could impede an accurate assessment of the true prevalence of healthy carriers in the population [5, 9, 10].

The PMS1 homolog 2, mismatch repair system component (PMS2, UniProt C9J167) represents an essential part in human DNA mismatch repair process. PMS2 combines with MLH1 in vivo to form the MutLα (MLH1-PMS2) heterodimer, which consist of N-terminal and C-terminal domains, as well as a linker region [11, 12], containing the Nuclear Localisation Signal (NLS) (Fig. 1) [13,14,15]. Currently, there is limited literature to describe how the MutLα complex is regulated. It is well established that MLH1 loss of expression is associated with premature PMS2 degradation. Additionally, several MLH1 missense variants in the C-terminal region have been linked to the loss of PMS2 protein expression and mismatch repair deficiency (MMR-d) [16]. Previous research suggests that phosphorylation of PMS2 in the C-terminal domain might play an important role, although specific amino acids involved in this process have not yet been identified [17].

Fig. 1

The PMS2 functional domains aligned to MANE transcript (NM_000535.7): Histidine kinase-like ATP-ase domain, MutL transducer hPMS2-like domain, Nuclear localization domain, MutL C terminal dimerisation domain

Historically, the only exonic positions regarded as critical for RNA splicing were those corresponding to acceptor and donor canonical splice sites, specifically the first and last three nucleotides of the exon [18, 19]. Further research uncovered the significance of additional factors in exon recognition, such as splicing regulatory elements (SREs), chromatin structure, transcription rate, and the secondary and tertiary structure of the transcript [20,21,22]. With the recent development and implementation of in silico tools and deep learning-based algorithms [23, 24], many exonic variants formerly classified as missense, nonsense, or silent are susceptible to reclassification, with some revealing a predicted impact on splicing [24,25,26,27,28]. In this scenario, the pathogenic mechanism usually involve a mixture of transcripts with abnormal splicing patterns and transcripts carrying the causative variant [29,30,31].

The vast majority of PMS2 variants documented in public databases to date are missense. Only a small proportion of these variants are, however, clinically significant, being classified as pathogenic (class 5) or likely pathogenic (class 4), according to ACMG 2015 guidelines [32,33,34]. Comparing to other LS-related genes, PMS2 truncating variants, including nonsense, frameshift and splicing variants, are less frequently encountered in clinical databases [34, 35]. Nevertheless, the major mechanism of disease is loss of function (LOF), with the latest ClinGen expert group assessment supporting haploinsufficiency as the major driver of disease [36]. In other human diseases, exonic variations are linked to splicing alterations in as much as 25% of cases [37], primarily involving exon skipping, but in PMS2-associated LS this aspect remains largely unknown. Previous research indicates that alternative splicing is a common phenomenon during the natural processing of the PMS2 gene [38]. Limited data is available on splicing role in tissue-specific and developmental stage-dependent PMS2 gene expression, as well as how variants altering splicing influence the protein function. Based on these observations, the main objective of our analysis was to narrow the existing diagnostic gap in PMS2-related LS by identifying potential significant splicing variants that could contribute to the underlying mechanisms of the disease. To assess the impact of splicing in silico, our analysis primarily focused on missense and short intronic germline variants related to Lynch syndrome and documented in public clinical databases. Bioinformatics tools, proven effective for other genes and previously employed in literature, were utilized to analyze both the wild-type DNA sequence and reported variants. To address the limitations of in silico tools and provide strength to our analysis, publicly available gene expression data was complementary employed for quantifying exon expression in a tissue-specific manner.

Results

Acceptor loss, the major mechanism for missense variants with predicted splicing impact in SpliceAI

Within the consulted databases (ClinVar and LOVD), a total of 2384 missense variants were documented (Fig. 2). Notably, the preponderant majority of these variants fell within the category of variants of uncertain significance (VUS), comprising 90.81% of the total. Variants with conflicting interpretations constituted 6.87%, while benign/likely benign variants accounted for 1.63%, and pathogenic/likely pathogenic variants were observed in 0.67%.

Fig. 2

Distribution of reported PMS2 missense variants across exons, stratified depending on ACMG 2015 classification. Red (P/LP – pathogenic/likely pathogenic variants), green (B/LB – benign/likely benign variants), orange (Conflicting variants), yellow (VUS – variants of uncertain significance). Each bar signifies the calculated density of variants, depending on the exon length and total number of reported variants

From all class 1–3 and conflicting missense variants (2384 variants) included in the study, 117 variants (4.90%) were anticipated to exhibit at least a mild effect on RNA splicing (DS > 0.2) and 34 (1.42%) a moderate or high impact (DS > 0.5) (Figure S1). Sixteen variants of uncertain significance were reported in 3′ canonical splice sites (first exonic position), with 7 (43.75%) having a potential impact (5 of acceptor loss type and 2 of acceptor gain type). Thirty-eight variants of uncertain significance and one conflicting variant were reported in 5′ canonical splice sites (last 3 exonic positions). Among these, 14 variants (35.89%) were expected to have an impact of donor loss type. When categorized by predicted consequences, the exons most commonly linked to specific splicing mechanisms were as presented in Fig. 3; Table 1. Additionally, wild-type ESRseq scores (HExoSplice in silico score calculated for enhancer splicing elements predicted in the wild type sequence) and ΔESRseq scores (HExoSplice in silico score calculated for enhancer splicing elements induced by genetic variants) corresponding to splicing variants (DS > 0.2) in exons enriched in acceptor loss (AL) and donor loss (DL) variants (Table 1) were compared to those from variants without such a prediction (Table 2, Figure S3).

Fig. 3

Class 1–3 and conflicting PMS2 missense variants with predicted splicing impact (DS > 0.2) in SpliceAI. Yellow (AG—Acceptor Gain), red (AL—Acceptor Loss), green (DG—Donor Gain), blue (DL—Donor Loss)

Table 1 PMS2 exons harboring missense variants with predicted splicing impact (SpliceAI)

Table 2 Wild-type ESRseq and ΔESRseq scores comparison between predicted splicing variants (DS > 0.2) and variants with no predicted splicing impact (DS < 0.2) in exons enriched in acceptor loss (AL) and donor loss (DL) variants

SpliceAI-visual, a valuable prediction tool for PMS2 complex short intronic variants

Out of 838 intronic variants reported in ClinVar, we discovered 71 non-point short genetic variants (< 50 bp). Upon filtering by clinical significance, 61 (85.91%) variants had conflicting interpretation or were interpreted as class 1–3 (benign, likely benign and VUS) according to ACMG 2015 criteria. Among these, 11 out of 61 variants affected the intronic canonical splice sites (the first 6 and last 3 intronic nucleotides). Notably, 14 variants were predicted to have a potential splicing impact when we used a lower threshold (DS < 0.2) for DS in conjunction with SpliceAI-visual. To the best of our knowledge, at the time when this manuscript was written, with one exception, no clinical or functional data were available for the mentioned variants. Additional details regarding the selected variants are available in the Table S1.

Canonical splice site interpreted as novel splice site by SpliceAI

The variant of uncertain significance NM_000535.7:c.1970_2006 + 9dup is a 46 bp duplication that spans over the 3′ end of exon 11 and exon-intron junction (Fig. 4A). In this particular case, relying solely on the DS provided by SpliceAI might suggest a mild increase in the strength of the canonical donor splice site. When the sequence was inspected using SpliceAI-visual, it was visible that SpliceAI located the original donor site in the duplicated region, showing only a slight decrease in strength (0.07). The canonical donor site was consequently interpreted as a novel splice site with DS score of 0.27. Upon further analysis, the duplicated region was anticipated to induce a frameshift and introduce a premature termination codon in the sequence. This event is expected to trigger transcript degradation via nonsense-mediated decay (NMD).

Variants increasing the strength of a weak canonical splice site

NM_000535.7:c.538-5_538-4del is a likely benign variant located in the proximity of acceptor site of exon 6 (Fig. 4B). As previously described, the 5′ end of exon 6 corresponds to a relatively weak splice site (REF score: 0.67). Within the exonic sequence, there is at least one alternative cryptic acceptor splice site that is stronger than the canonical site (REF score: 0.81). While the variant DS (0.28) suggests a mild increase in strength of the acceptor splice site, the cumulative effect renders the canonical splice site stronger than cryptic sites (ALT score: 0.94). This, in turn, could consequently alter the natural proportion of Δ6 and Δ6p transcripts with potential clinical impact. Similarly, variant NM_000535.7:c.538-12dup with conflicting interpretation of pathogenicity increases the strength of the same acceptor splice site (ALT score: 0.78), despite having only a mild DS (0.12).

Intronic inclusion and premature termination predicted by SpliceAI-visual

The variant NM_000535.7:c.354-18_354-15dup previously classified as likely benign, constitutes a 4 bp intronic duplication near the 3′ splice site of exon 5 (Fig. 4C). SpliceAI indicates a mild acceptor gain effect (DS: 0.17) that could be easily filtered out using a standard DS cutoff value of 0.2. When using SpliceAI-visual, we observed that the variant enhances the strength of an intronic cryptic acceptor site (ALT score: 0.57). This site could compete with the canonical site and lead to intronic inclusion. Using the alternative splice site is important in this context since the variant induces a frameshift and includes a premature termination codon naturally present in the intronic sequence. This, consequently, is predicted to induce NMD and eventual transcript degradation.

Fig. 4

SpliceAI-visual predictions (IGV interface, MobiDetails) for PMS2 short intronic variants: (A) NM_000535.7:c.1970_2006+9dup, (B) NM_000535.7:c.538-5_538-4del, (C) NM_000535.7:c.354-18_354-15dup. REF—reference (wild-type) score, ALT—alternative (variant) score. Vertical blue bars signify donor sites and orange bars signify acceptor sites

Bioinformatics assessment of donor and acceptor splice sites strength

Canonical exon-intron junction motifs were evaluated using eight alternative software tools (Table 3, Figure S4). As was previously indicated [39], splice sites with scores falling below the lower boundary of the 90% confidence interval (90% CI) were deemed weak. Eleven exons displayed at least one putative weak splicing signal, with exons 2, 3, 4, 5, 6, 7, 8, 11, 12 and 14 consistently reported in at least three out of eight predictions. Exons 9, 10, 13, and 15, were lacking identifiable weak canonical splice sites.

Table 3 Summary of PMS2 exons predicted to harbor strong and weak canonical splice sites

Exonic splicing regulatory elements (SREs) predictions: ESEs, ESSs and ESS/ESE ratio

The coding regions of all PMS2 exons, except from the first and last exons, were assessed utilizing motif matrices in ESEfinder to detect potential binding sites for SR proteins: SF2/ASF, SF2/ASF (IgM-BRCA1), SC35, SRp40 and SRp55 (Fig. 5, Table S3). Significance in the analysis was attributed solely to motifs with high scores surpassing the standard thresholds recommended by developers. Exons 3, 7, 10, and 13 exhibited the lowest overall density of ESE motifs, with an isolated overrepresentation of binding sites for SRp40 in exon 10. The findings from the parallel analysis, employing both HOT-SKIP and HExoSplice, are comprehensively outlined in Table 4 and Table S4.

Fig. 5

Exonic density of SR proteins motifs identified using ESEfinder 3.0 software. Red lines represent the 95% CI for means. Exons outside the interval are depleted (below the lower bound of 95% CI) or enriched (above the upper bound of 95% CI) in SR protein motifs and were noted in the graph

Table 4 PMS2 exons predicted to be prone to exon skipping, based on ESEs and ESSs densities and ESS/ESE ratio

High density/ratio values are considered those above the upper bound of 90% CI, whereas low density/ratio values are below the same interval.

Exons displaying consistently high predictions for ESE (exons 5, 11, 14) or ESS (exons 6, 7, 10) densities in all three SREs-specific software (Table S4), were further examined based on their wild-type ESRseq scores returned by HExoSplice, in order to explore motifs distribution across exons (Fig. 6, Table S5). Correlation between exonic position and ESRseq score in exon 5 (R = − 0.34, p = 0.001) and exon 14 (R = 0.30, p = 0.006) was weak-moderate and statistically significant. In exon 5, ESEs were located preferentially in the first 50 bp, while ESSs were found in the latter half of the exon. This pattern was reversed in exon 14.

Fig. 6

Exonic distribution (HExoSplice) of SREs motifs in exons with high density of ESEs or ESSs (exons 5, 6, 7, 10, 11 and 14). ESRseq scores > 0 correspond to ESE motifs and ESRseq scores < 0 represent ESS motifs. The red color signifies the high density of SREs motifs in a genomic region

PMS2 expression data—GTEx database and RefSeq coding transcripts

The available data regarding PMS2 gene expression revealed that exons 11, 13, and 14 had a significantly higher median count per base in all tissues included in the study (Table 5, Table S6). On the other hand, exons 1, 2, and 15 consistently exhibited low expression in almost all tissues. 60 protein coding transcripts were included in the study (Table S7, Figure S5). The most frequent splicing variants in studied transcripts were presented in Fig. 7.

Table 5 Mean RNA expression (median count/base) of PMS2 exons in Lynch syndrome-related tissues

Fig. 7

Splicing variants encountered in PMS2 RefSeq coding transcripts. Δ—complete or partial (adjacent to canonical splice sites) exonic deletion, ▼—intronic sequence inclusion (adjacent to canonical splice sites), p—acceptor site shift, q—donor site shift

Discussion

Low level of exonic splicing variants in PMS2 predicted by SpliceAI

In our study, only 4.9% of PMS2 missense variants were predicted to have a significant splicing consequence (DS > 0.2) using SpliceAI. This value is lower than expected given the data available in literature which states that up to 25% of exonic disease-causing variants may disturb exon definition in mature transcripts [37]. There are several potential explanations for this result: (1) A lower performance of bioinformatics predictions in exonic positions outside the canonical splice sites [19, 40]; (2) Variants situated deep in the exonic regions are less likely to influence exon inclusion [41]; (3) The DS threshold > 0.2 could impact variant classification in PMS2 gene, clinically significant variants being reportedly overseen in other genes when a standard cutoff value was used [24]; (4) Alternative splicing have a secondary role in the biology of PMS2 gene, according to what was previously reported in literature [42].

Exons harboring weak canonical splice sites may require regulatory elements for proper definition

Critical for exon definition, canonical splice sites may demand additional signals, such as ESEs or ESSs, when their sequence is less conserved [18, 23]. Not surprisingly, in the past, other analysis of variants that disrupt SREs supported the same idea of SREs exhibiting activity predominantly in exons having splice sites of weak or moderate consensus [40]. In accordance with these observations, in our study, exons 2, 4, 5, 6, 8, 12, 14 predicted to have weak splice sites exhibited a high frequency of at least one ESE motif in ESEfinder. Exon 3, however, although predicted with a weak 3′ splice site, had a lower than average level of ESEs in ESEfinder and HExoSplice, but a significantly higher one in HOT-SKIP, suggesting a potential for alternative splicing events [43]. Notably, the expression of exon 3 was below average in all tissues, with whole blood and medullary kidney showing significantly low expression levels. In contrast, the median read count per base in bladder tissue was above the average. When taken globally, in all exons with a weak 5′ss or 3′ss, except for exon 4, ESE or ESS motifs were overrepresented in at least one prediction. Nevertheless, in exon 4, ESEfinder detected a high density of SC35 and SF2/ASF (IgM-BRCA1) specific domains, but the overall ESEs density was above the average, reaching the upper bound of the confidence interval.

PMS2 exons with weak splice sites may be prone to exon skipping

As indicated by the low ESEs frequency, high ESSs frequency, or high ESS/ESR ratio [44, 45], exons 2, 3, 6, 7 and 8 were predicted by at least one SRE software as being potentially prone to exon skipping. This association is plausible, given similar in silico observations reported in the literature for other genes [41], which were subsequently confirmed in functional studies [39]. Indeed, when consulting transcripts reported to date, high number of wild type and mutant isoforms aligned with this prediction [42, 46]. Bulk RNA-seq data was concordant over exons 2, 3 and 7, with the mention of variations in exon 3 and 7 that were tissue-dependent. On the other hand, exons 6 and 8 displayed an average median read count per base across all tissues, as indicated by GTEx database. Moreover, exon 8 was shown to harbor higher than average levels of SF2/ASF sites, but the overall density of SR proteins-binding domains remained low. SpliceAI was concordant with canonical splice site software regarding exon 8, indicating a high density of acceptor loss variants, which may serve as an indirect proof of the same biological event. Indeed, Δ8_11 transcript, naturally occurring by a splice acceptor shift mechanism, was reported [42]. Other similar ‘multi-cassette’ RefSeq transcripts were curated, suggesting once again the role of alternative splicing in processing of this exon. Intriguingly, exon 4, also featuring a weak donor site, undergoes alternative splicing in certain transcripts, while exons 5 and 14 utilizes distinct acceptor sites [46], an aspect also foreseen by SpliceAI-visual (Figure S2), which identifies strong cryptic acceptor sites. Recently, several exonic and intronic variants associated with exon 4 skipping were reported in literature, highlighting even more the relevance of alternative splicing in PMS2-related LS [40, 46, 47].

Other exons potentially sensitive to skipping with no weak consensus were 9, 10 and 13. Exon 10 was mentioned by all 3 tools, with a low ESEs incidence in ESEfinder and high ESSs levels in HOT-SKIP and HExoSplice. On the other hand, exon 10 is one of the most highly expressed exons across all tissues. Of note, several protein-coding RefSeq transcripts revealed, however, an isolated exclusion of this exon. Minor Δ10 transcripts were also evident in vitro assays [46]. Exon 9 is sometimes skipped in a ‘multi-cassette’ event with other exons. However, from our knowledge, there are no natural transcripts where exon 9 is skipped individually. In this case, ESRs predictions were conflictual, with ESEfinder indicating a high density of ESE motifs, while the other tools supported a low level. Similarly, exon 13 is a constitutive exon, being included and highly expressed in the majority of transcripts. This comes as no surprise, given that it encodes the MutL C terminal dimerization domain, a region known for its high conservation in the protein [48, 49].

High ESE levels concordantly predicted in exons 5, 11 and 14 – critical for PMS2 function

The vast majority of ESE domains in exon 5 are located next to the acceptor site in HExoSplice. Interestingly, an alternative 3’ss located proximally in the intronic sequence (and predicted by SpliceAI-visual—Figure S2) was observed in at least one naturally occurring transcript. This distribution prompts questions regarding the role of ESEs in the region, particularly considering that 3 out of 8 predictions indicated a 3’ weak splice site [50]. Of note, exons 2–5 codify N-terminal ATPase domain of PMS2 (HATPase_c_3, InterPro, PF08676), which is important for maintaining mismatch repair proficiency. However, it is worth highlighting that this domain may be working in an asymmetric manner with the similar domain in MLH1, with the later appearing to play a more decisive role in this biological process [51,52,53].

By comparison, in exon 11, ESEs and ESSs hexamers are evenly distributed throughout the exon, with a region with high ESSs density and a relative ESEs depletion near the 3’ss. Exon 11, being a long exon, suggests that ESEs may play a fundamental role in facilitating the accurate selection of exonic boundaries by the splicing machinery [43]. Additionally, certain isoforms include a shorter variant of this exon due to a shift in the 5′ splice site, in contrast with the presence of a relatively strong 5′ donor consensus. Similarly, SpliceAI predicted a significant number of variants in exon 11 that could lead to donor site gain or loss. At least one transcript, Δ11q_14p, has been reported to occur via splice donor shift process [42]. An analogous event was noticed in atypical CMMRD related to a benign missense founder variant that created a novel donor splice site in intron 11, generating a 5 bp deletion frameshift at the exon 11–12 junction [54, 55]. Remarkably, exon 11 is rarely spliced out in reported transcripts, which indicates a fundamental biological significance of the regional coding sequence. Upon delving into the GTEx exon expression data, we observed a consistent trend of high expression for exon 11 across all investigated tissues. Upon reviewing the literature, we discovered that PMS2 interacts with MLH1 in a critical area spanning amino acid residues 675 to 850 [11]. This sequence defines the C-terminus of the protein, a region sensitive to phosphorylation and involved in regulating the degradation of PMS2, that overlaps with exons 12–15 [17]. Additionally, missense variants within exon 11, located in close proximity to this hotspot region, may impact MutLα heterodimer formation by altering the binding affinity of the protomers and nuclear localization domain [56].

Within exon 14, ESSs are primarily located at the 5’ end (3’ acceptor site), whereas ESEs are complementarily positioned in the middle and 3’ end (5’ donor site). Since exon 14 was consistently predicted to have a weak acceptor site, and confirmed by wild type and mutant Δ14p transcripts [46], the high ESEs incidence may be relevant for its inclusion in mature transcripts, as it was reported in other exons with similar properties [43]. As mentioned earlier, exon 14 encodes the MutL C terminal dimerisation domain (MutL_C, InterPro, IPR014790) and, thus, critically contributes to MutLα heterodimerization [48, 49]. As anticipated, the RNA expression analysis revealed, once again, a consistent pattern of elevated expression across all considered tissues.

High ESS levels concordantly predicted in exons 6, 7 and 10

Both HOT-SKIP and HExoSplice indicated potential elevated levels of exonic splicing silencers in exons 6, 7 and 10. When we visually inspected exon 6 in HExoSplice, ESSs were strikingly located in the 5’ half of the exon. In combination with ESE depletion in the same region, it could explain, at least partially, the acceptor splice shift, Δ6p and Δ6 transcripts described in literature [42, 46] and present in RefSeq data. SpliceAI-visual identified at least 2 exonic cryptic acceptor sites in the region (Figure S2) that may play a decisive role in the process. In line with this, SpliceAI predicted a significant number of acceptor gain and acceptor loss variants in exon 6 based on the variants collected in this study. Moreover, the median ΔESRseq score of predicted splicing variants was negative and significantly lower from other reported variants in this exon. In exon 7, ESSs are located preferentially at both 5’ and 3’ ends of the exon, middle of the exon being enriched in ESEs instead. Besides the predicted weak donor site, we suspect that this particular distribution of SREs may provide a possible explanation for the observed exon 7 skipping and intron 7 inclusion in some transcripts [42]. In exon 10, ESSs are relatively evenly distributed, interspersed with ESEs, showing no specific pattern.

Limitations

The present study conducted a bioinformatics evaluation of the significance of alternative splicing in the expression of the PMS2 gene. This analysis is relevant given the numerous PMS2 variants of uncertain significance reported, that complicate the translation of DNA sequencing data into clinical practice. However, we acknowledge that the present study has several limitations, primarily related to the fact that the analysis was conducted exclusively in silico. Despite including bioinformatics tools commonly used and verified by other authors, future experimental data are still required to validate the results. Additionally, the included expression data may exhibit certain known biases [57, 58]. The analyzed transcripts are solely protein-coding, excluding non-coding transcripts, which may not fully capture the diversity and complexity of alternative splicing. Furthermore, while the primary focus of the current study was represented by the PMS2 coding variants, the potential importance attributed to noncoding variants and antisense transcripts in exon selection should not be disregarded, warranting further investigation.

Materials and methods

Reference sequence and variant nomenclature

PMS2 variants were described according to Human Genetic Variation Society guidelines (https://hgvs-nomenclature.org/stable/). MANE Select transcript (NM_000535.7, ENST00000265849.12) was considered the reference sequence, position c.1 corresponding to the first coding nucleotide. The variants analyzed were retrieved from the public databases ClinVar (https://www.ncbi.nlm.nih.gov/clinvar) and LOVD (https://databases.lovd.nl/shared/genes/PMS2) on 4 December 2023. These are the most comprehensive clinical genomic databases currently available for free. Variants with no reported classification based on ACMG criteria [32] have been excluded from the analysis. Variant annotation was performed using Ensembl Variant Effect Predictor (VEP) (https://ensembl.org/Homo_sapiens/Tools/VEP/). For exonic regions, we focused on nonsynonymous missense variants. Truncating variants, known to be deleterious, and synonymous variants, very rarely reported in consulted databases, were excluded. Among the selected intronic variants, we specifically included short genetic variants (< 50 base pairs). Complex and large variants were beyond the scope of our study.

Tissue-specific bulk RNA expression data and public available transcripts

Protein-coding PMS2 transcripts under analysis were retrieved from NCBI RNA reference sequences collection (RefSeq), a public database that provides an extensive and carefully curated collection of sequences [59, 60]. RefSeq transcripts available on UCSC genome browser (https://genome.ucsc.edu) were included in the analysis. To further enhance our understanding, the GTEx database was employed for quantifying exon expression using median read per base provided in a tissue-specific manner [58]. The analysis focused on tissues more commonly affected in Lynch syndrome.

Bioinformatics analysis of splicing impact and statistical analysis

Several bioinformatics approaches were employed to assess the impact of PMS2 variants on RNA splicing. The strength analysis of canonical acceptor and donor splice sites was conducted using freely available online tools, including ESEfinder 3.0 for splice sites (https://esefinder.ahc.umn.edu/), FSplice (http://www.softberry.com/), MaxEntScan (http://hollywood.mit.edu/), NetGene2 (https://services.healthtech.dtu.dk/) and NNSplice (https://www.fruitfly.org/seq_tools/splice.html), SpliceAI [23, 24, 61], Spliceator (https://www.lbgi.fr/spliceator/), SpliceRover (http://bioit2.irc.ugent.be/rover/splicerover). The potential splicing impact of reported variants was estimated using SpliceAI, one of the most proficient deep learning-based tool reported to date [23, 24, 61]. Delta scores (DS) greater than 0.2 were used to screen for splice-altering variants, as outlined in the original publication [23]. To enhance the characterization of intronic variants, we utilized the SpliceAI-visual IGV interface [62] accessible on MobiDetails (https://mobidetails.iurc.montp.inserm.fr). Moreover, SpliceAI-visual was used to predict the presence of cryptic splice sites within the wild type sequence. Splicing regulatory elements (SREs) were screened across the coding regions using three alternative strategies: ESEfinder 3.0 for SR proteins (https://esefinder.ahc.umn.edu/), HOT-SKIP (https://hot-skip.img.cas.cz/) and HExoSplice (http://bioinfo.univ-rouen.fr/HExoSplice/). Exon skipping was considered in either case based on the resulting low putative ESEs and high ESSs densities, as well as increased ESS/ESE ratios [44, 45]. All splicing predictions were done using default specifications for each in silico tool (Table S8). A brief flowchart for the methods used in this study is shown in Fig. 8. Satistical analysis was performed using MedCalc^® Statistical Software version 22.019 (MedCalc Software Ltd, Ostend, Belgium; https://www.medcalc.org; 2024). Statistical significance was assigned for p-value < 0.05. Except where otherwise stated, a 95% confidence interval (CI) was employed. The Student’s t-test was used to statistically evaluate numerical and normally distributed values. The Mann-Whitney U test was used as a non-parametric test in non-normally distributed data.

Fig. 8

An overview for the methods used in this study: a) data extraction from clinical databases – ClinVar and LOVD, and variant annotation using Ensembl (VEP – variant effect predictor) b) splicing predictions tools employed for canonical splice sites and SREs predictions

Conclusions

This study highlights the importance of in silico tools to further guide the PMS2 mRNA analysis and potentially increase the diagnostic yield in Lynch syndrome. Bioinformatics software used underlines the frequency of splicing alterations associated to PMS2 gene, providing a possible explanation for the current underdiagnosis in PMS2-associated LS. In this regard, we identified several missense and short intronic variants that were predicted to have a splicing impact, making them good candidates for future functional analysis. Exons 2, 3, 4, 5, 6, 7, 8, 11, 12 and 14 were shown to have weak canonical splice sites in at least three prediction tools. Moreover, we observed the importance of ESEs and ESSs in the PMS2 gene, their density and localization across the sequence, providing valuable input about splice site selection and overall exon expression. ESE motifs were more prevalent in highly expressed exons 5, 11 and 14, while ESS motifs were overrepresented in exons 6, 7 and 10. Some limitations of in silico tools emerged, however; therefore, additional functional data are necessary to evaluate the biological significance of our computational observations.