High-throughput analysis revealed mutations’ diverging effects on SMN1 exon 7 splicing

ABSTRACT

Splicing-affecting mutations can disrupt gene function by altering the transcript assembly. To ascertain splicing dysregulation principles, we modified a minigene assay for the parallel high-throughput evaluation of different mutations by next-generation sequencing. In our model system, all exonic and six intronic positions of the SMN1 gene’s exon 7 were mutated to all possible nucleotide variants, which amounted to 180 unique single-nucleotide mutants and 470 double mutants. The mutations resulted in a wide range of splicing aberrations. Exonic splicing-affecting mutations resulted either in substantial exon skipping, supposedly driven by predicted exonic splicing silencer or cryptic donor splice site (5′ss) and de novo 5′ss strengthening and use. On the other hand, a single disruption of exonic splicing enhancer was not sufficient to cause major exon skipping, suggesting these elements can be substituted during exon recognition. While disrupting the acceptor splice site led only to exon skipping, some 5′ss mutations potentiated the use of three different cryptic 5′ss. Generally, single mutations supporting cryptic 5′ss use displayed better pre-mRNA/U1 snRNA duplex stability and increased splicing regulatory element strength across the original 5′ss. Analyzing double mutants supported the predominating splicing regulatory elements’ effect, but U1 snRNA binding could contribute to the global balance of splicing isoforms. Based on these findings, we suggest that creating a new splicing enhancer across the mutated 5′ss can be one of the main factors driving cryptic 5′ss use.

KEYWORDS:

• SMN1
• cryptic splice sites
• U1 snRNA
• splicing-affecting mutation
• 5′ss

Introduction

RNA splicing is an essential mRNA processing step, producing one or more transcript isoforms, and is tightly regulated to guarantee correct protein expression. Specific RNA sequence motifs, especially donor and acceptor splice sites (5′ss and 3′ss, respectively) and branch point, are required for correct exon recognition and subsequent spliceosome assembly. Recognizing an exon is modulated by cis-acting splicing regulatory elements (SREs). These elements are classified as exonic (ESEs) or intronic (ISEs) splicing enhancers and exonic (ESSs) or intronic (ISSs) splicing silencers. These sequences recruit diverse splicing factors that modify exon recognition and spliceosome assembly [1].

In silico evaluation of the mutations’ effect in SREs is a challenge. Despite an obvious improvement in many algorithms, their real use for predicting mutation impact on RNA splicing efficiency is limited because of exon recognition complexity or a higher degeneracy of the splicing factors’ binding sequences [2–5]. The algorithms are thus dependent on extensive datasets [6–8] including those based on experimental work. Expanding the possibilities of high-throughput technologies can thus contribute to improving these algorithms. Moreover, both direct and non-direct proof of mutations’ impact on RNA splicing can also be helpful for clinical geneticists.

Exon recognition itself is dependent on 5′ss interaction with U1 snRNA [9,10], which is typically mediated by Watson-Crick (WC) base-pairing. Vertebrate 5′ss display a relatively high 5′ss sequence diversity. While consensus sequence CAG/GTAAGT can form WC base-pairs, other 5′ss U1 snRNA duplexes contain mismatched base pairs and the 5′ss selection is related to the thermodynamic stability of these duplexes [11]. Alternatively, the duplex formation can occur in shifted or bulged register [12–14]. Later in the splicing process, U1 snRNP must be replaced by U4/U6-U5 tri-snRNP, which exactly defines the cleavage site. However, a U4/U6-U5 tri-snRNP binding site does not need to overlap the U1 binding site [15,16], which enables the use of a relatively weak 5′ss.

To better understand the principles of exon recognition, we chose SMN1 exon 7 as a model system and performed a systematic analysis of all possible exonic and six intronic variants using next-generation sequencing (NGS) technology. Inactivating mutations in the SMN1 lead to spinal muscular atrophy (SMA) development and, among others, truncation of the last 16 amino acids caused by exon 7 skipping, inactivating the protein [17]. Thus, even the mutations that do not influence the amino acid composition, such as e-3C>T, can result in disease outbreak [18]. Promisingly, the impact of abnormal splicing patterns causing mutations can be rescued by exogenous intervention inducing cryptic 5′ss usage and thereby maintaining exon 7 inclusion [19].

SMN1 exon 7 inclusion itself is regulated by a set of SREs. Among the first detected variants included in exon 7 recognition was a nucleotide exchange C/T in the sixth exonic position, which is a natural polymorphism differentiating SMN1 from SMN2, a paralogous gene with decreased exon 7 inclusion [20]. The SMN2 variant induces exon skipping (ES) through hnRNP A1/A2 binding sequence creation [21], although it seems that the SMN1 variant can bind splicing enhancer SRSF1 [22]. The middle part of SMN1 consists of an A/G rich region binding hTra2β1 [23–25], PSF [26] and indirectly SRp30c [25], hnRNP-G [27] or TDP-43 [28]. At the 3′ end of SMN1 exon 7, the RNA sequence can form a terminal stem-loop structure, stabilization of which promotes ES [29]. Other important elements are located in introns 6 and 7. Miyajima et al. [30] identified a cis-element within intron 6, whose deletion led to increased exon inclusion. Similarly, SMN1 intron seven contains additional elements promoting/repressing exon 7 inclusion. Thus, TIA1/TIAL1 binding to U-rich elements located downstream of authentic exon 7 5′ss supports exon 7 inclusion [31], while hnRNPA1/A2 binding to intronic motif ISS-N1 represses it [32].

In this work, we applied a high-throughput approach to elucidate the impact of all possible single nucleotide mutations on SMN1 exon 7 splicing in the artificial context of the surrounding exons in the minigene. This approach revealed different mutations’ diverse effect in the 5′ss region on ES and the use of cryptic 5′ss. To explain such a wide variety between different mutations, we tested two hypotheses. First, we hypothesized that U1 snRNP preserved the ability to bind mutated 5′ss and helped to recognize distant cryptic 5′ss due to a shifted U1 snRNP binding mechanism. The second hypothesis assumed that a mutation in 5′ss resulted in a new ESE creation, subsequently supporting the U1 snRNP binding to an extremely weak 5′ss.

Results

Design and properties of mutagenesis library

To systematically map splicing affecting mutations in SMN1 exon 7, a three-exonic minigene with the middle SMN1 exon 7 was prepared and subsequently mutated in the whole coding region and neighbouring three intronic nucleotides. To enable the exact coupling of detected splicing isoform with the corresponding causal mutation, first exon was barcoded with 14-nt barcode sequence (Figure 1; Supplementary Figure S1).

Figure 1. SMN1 exon 7 in the experimental system. (A) Scheme of minigene containing exon 7 of the SMN1 gene. Boxes represent exons, lines introns. Part of the first exon is shown as a nucleotide sequence, where pl corresponds to exonic sequence before barcoding and bc after barcoding. Underlined nucleotides represent restriction sites used during plasmid library construction. BglII and Hind III restriction sites are used for exon 1 replacement. During this procedure, original exon 1 is replaced with the barcoded one and NheI and AscI restriction sites are lost, which can be used to control successful barcoding. Arrows represent primers used for amplification preceding NGS. Primers fpet and rint are designed for mutant characterization (after preceding intron excision by PacI and AsiSI) and primers fpet and rpet for splicing analysis. Cryptic 5′ss are depicted including their positions. (B) Schematic representation of detected transcripts. (C) Frequency of particular transcript isoforms in wt minigenes detected using NGS. A frequency is depicted as a bar-plot, where boxes correspond to quartiles and whiskers to 0.1 and 0.9 percentiles, respectively. (D) Wt transcript profiling using agarose electrophoresis for endogenous SMN1 and a minigene. SMN1 exon 7 skipping is visible only as a very blurred band

The final plasmid library consisted of about 3,800 differently barcoded minigene variants sufficiently expressed in transfected HeLa cells (Supplementary Figure S1). About 1,660 plasmids represented wild-type (wt) minigenes and the others represented mutated minigenes with at least one mutation inside the mutated region. In the case of single mutations, 180 different mutations were expected. A typical mutation was represented by 6 (median) different plasmids and the highest number of plasmids with an identical mutation was 25 (the number of plasmids per single mutation are depicted in Supplementary Figure S2). Additionally, single mutations outside the targeted region were also detected (Supplementary Figure S3). From these ones, only e+4A>G was added to the subsequent analysis because of its interesting transcript profile (Figures 2 and 3).

Figure 2. The effect of single nucleotide substitutions on splicing. A diagram depicts frequencies of ES (at the upper panel) and authentic 5′ss usage (at the lower panel) in transcripts containing SMN1 gene exon 7 of all analysed mutations as a heat map, where red and blue colours represent aberrant transcripts (ES at the upper panel and usage of cryptic and de novo 5′ss at the lower panel). Mutations are positioned in the panel according to their position inside the mutated region (depicted along the x-axis with marked wt nucleotide) and mutation (y-axis). Black colour in the panel represents wt nucleotides, green colour non-analysed mutations. Checkmark indicates a statistically significant decrease in the frequency of appropriate transcript variants (P < 0.01) and asterisks mutations with a statistically significant decrease of frequency above the threshold of 1%

Figure 3. Proportion of transcripts in mutants causing de novo and cryptic 5′ss usage. Exon skipping (ES), exon inclusion (EI; normal splicing isoform only), cryptic 5′ss (Cr1, 2, 3) and de novo 5′ss (Dn) splicing variants frequencies are depicted. Corresponding statistics, incl. standard deviations and p-values for differences among mutations, is available in Supplementary Table S1

High-throughput analysis captured a wide variety of splicing isoforms

Owing to the design of our method, we were able to detect 9 different transcripts and analyse complex splicing changes. In the wt minigene, the normally spliced (exon inclusion) transcript prevailed, but ES was also detected, which corresponded to endogenously expressed SMN1 (Figure 1). Among other detected transcripts, using cryptic (Cr) and de novo splice sites occurred as well. Thus, the detected transcripts were: normally spliced isoform that used both authentic splice sites, ES, transcripts using cryptic 5′ss Cr1 (located in the intronic position e+24; Figure 1), Cr2 (position e+52) or Cr3 (exonic position 30, i.e. c.768), transcripts using de novo 5′ss (exonic position 26 [c.764] due to the mutation 27A>T; or in the exonic position 48 [c.786] due to the mutation 48T>G). The least frequent transcripts were those using the intronic cryptic 3′ss (intronic position e-36), or exonic cryptic 3′ss (exonic position 7, [c. 745]). Cryptic 3′ss use was extremely rare (totally about 2/100,000 and 13/100,000 reads, respectively) and therefore these were not considered in the subsequent analyses.

The impact of single mutations on SMN1 exon 7 RNA splicing

One of the main advantages of this approach was the possibility to test the effect of single mutations, which can help to identify regions carrying functional SREs. Out of 181 evaluated single mutations, 58 showed significantly increased ES proportion (P < 0.01; Figure 2, upper part), but only 37 of them exceeded a threshold of 1%. Nine of these 37 mutations were located in the authentic 3′ss. The mildest effect among 3′ss mutations was observed in e-3C>T, where ES reached 11%, which perfectly correlated with a moderate 3′ss strength decline as predicted by MaxEnt (Table 1). Excluding splice site mutations, the strongest effect was detected for 51A>G (ES 31%) and 23A>T (ES 24%). For most of the splicing-affecting mutations, we predicted the cause of this disruption using various prediction tools (SpliceAid2 [33], EX-SKIP [34], ESRseq [6], MaxEnt [35], HEXplorer [36]) and combined it with data from the literature. All these findings, including the corresponding ESRseq and Hexplorer scores, are summarized in Table 1 and discussed in depth later.

Table 1. Summary of splicing affecting mutations revealed by NGS. For the purposes of this table, the percentage of aberrant splicing was estimated from NGS data, considering a total number of aberrant reads and subtracting the average number of chimeric (normally spliced) reads. Probable causes of mutations were predicted based on the profile of created minigene transcripts, described splicing regulations and in silico predictions using MaxEnt, Ex-SKIP and others. MaxEnt scores for individual splice sites are presented in the notes column

Mutations	Percentage of aberrant transcripts	The most probable cause of splicing defects	Notes
e-3C>R | c.739-3C>R	~ 100%	Decreased 3′ss strength	MaxEnt 8.26 [C>A], −2,5 [C>G] vs. 10.92 [wt]
e-3C>T | c.739-3C>T	~ 10%	Decreased 3′ss strength	MaxEnt 10.42 vs. 10.92 [wt]; rs772466166;^18
e-2A>B | c.739-2A>B	~ 100%	Disruption of 3′ss AG dinucleotide	rs141760116
e-1G>H | c.739-1G>H	~ 100%	Disruption of 3′ss AG dinucleotide	^40
1G>T | c.739G>T	~ 9%	Decreased 3′ss strength	MaxEnt 9.81 vs. 10.92 [wt]; rs77969175
6C>T | c. 744C>T	~ 4%	ESS creation (hnRNPA1)	^Citation20,Citation70
23A>T | c.761A>T	~ 24%	ESS creation	HEXPLORER score difference: −159.4; ESRseq score difference: −4.19
27A>C | c.765A>C	~ 9%	ESS creation	HEXPLORER score difference: −99.89; ESRseq score difference: −3.94
27A>T | c.765A>T	~ 94%	de novo 5′ss usage	MaxEnt 10.29
29G>T | c.767G>T	~ 17%	ESE (hTra2β1) disruption/ESS creation	HEXPLORER score difference: −19.45; ESRseq score difference: 0.57;^41
32G>A | c.770G>A	~ 6%	Increased Cr3 strength	MaxEnt 7.46 vs 3.68 [wt]
34T>G | c.772T>G	~ 7%	Increased Cr3 strength	MaxEnt 9.05 vs 3.68 [wt]
39T>C | c.777T>C	~ 12%	ESS creation/ESE disruption	HEXPLORER score difference: 23.68; ESRseq score difference: 0.94
48T>G | c.786T>G	~ 46%	de novo 5′ss usage	MaxEnt 9.72;^29
51A>G | c.789A>G	~ 32%	ESS creation	HEXPLORER score difference: −62.33; ESRseq score difference: −2.53;^25
e+1G>H | c.792+1G>H	~ 100%	Disruption of 5′ss GT dinucleotide	^37
e+2T>V | c.792+2T>V	~ 100%	Disruption of 5′ss GT dinucleotide
e+3A>B | c.792+3A>B	~ 100%	Decreased 5′ss strength	MaxEnt 6.05 [A>G], −1.36 [A>T], −2.51 [A>C] vs. 8.57 [wt]
e+4A>G | c.792+4A>G	~ 90%	Decreased 5′ss strength	MaxEnt 3.34 vs. 8.57 [wt]

Fifteen mutants showed significantly increased cryptic or de novo splice site use (Figures 2 and 3). There were two single mutations (32G>A and 34T>G) characterized by increased Cr3 use (Figure 3), both showing increased Cr3 strength (MaxEnt 7.46, and 9.05, respectively vs. 3.68 of wt). The strongest effect on splicing was observed in the 27A>T variant that forms a strong de novo 5’ss (MaxEnt 10.29) used in 88%. Another de novo 5′ss, created by a 48T>G mutation, showed MaxEnt score 6.6 and use in 21%.

The most remarkable changes were detected in mutants affecting the authentic 5′ss (Figure 3). Surprisingly, the mutations’ effects differed extensively. While e+3A>G and e+3A>T mutations caused massive ES, all e+2 mutations showed only slight changes in the ES frequency, but the use of both intronic cryptic sites Cr1 and Cr2 was greatly increased. All mutations in position e+1 and mutation e+3A>C had an intermediate pattern of these splicing isoforms.

NGS results verification

To prove NGS data reliability, the impact of selected single mutations was verified using classical minigene assay. Capillary electrophoresis was used to quantify minigene transcripts. Splicing profiles for all tested minigenes (wt, 51A>G, 23A>T, 6C>T and e+2T>C) (Supplementary Figure S5) corresponded to the profile obtained in the NGS experiment with a minor shift towards a higher abundance of normally spliced minigene in the NGS experiment. This phenomenon was observed in all tested mutations and the wt minigene, including mutations in AG and GT dinucleotides of the 3′ss and 5′ss, respectively. These nucleotides’ mutations were expected to lead to complete splicing disruption, as was demonstrated in the verification experiment (Supplementary Figure S5). We hypothesized that these observed ‘contaminations’ were caused by chimeric sequences formed from barcoded exon 1 of the minigene and the middle exon belonging to a different transcript. Therefore, we compared the incidence of non-mutated SMN1 read in the population of transcripts with barcodes belonging to single nucleotide mutations (Supplementary Figure S6). This type of chimeric sequence mostly accounted for 3–10% of the reads, which correlated well with the number of normal transcripts observed in AG and GT mutants. These findings show that chimeric sequence creation could explain the slight shift in transcript isoform frequency.

Since the most abundant transcript is normal exon inclusion, its frequency was slightly overestimated at the expense of all other isoforms. However, all previously detected splicing affecting mutations (incl. 6C>T; Table 1) were confirmed in the NGS experiment and despite this slight inaccuracy, changes in the aberrant transcript frequency remained statistically significant.

Comparing aberrant transcript frequencies with data from literature showed that the effect of mutations on RNA splicing efficiency is less conspicuous in our experiments. To test whether this effect can be caused by the selected cell line used for NGS, we compared the wt transcript profile, 23A>T and 51A>G both in HeLa and HEK-293 cells (Figure 4). The latter are more often used to evaluate splicing defects in SMN1 exon 7. Indeed, ES was more frequent in HEK-293 cells for all these minigenes, but the difference reached 10% at most.

Figure 4. Comparison of transcript levels in different cells. The results were acquired with standard minigene assay and exon skipping (ES) and exon inclusion (EI) rates were measured using capillary electrophoresis

Altogether, it seems that cells used for this experiment, in combination with the specific pET derived minigene context impacted on the rate of a transcript isoform. Yet our experiment appeared to be very sensitive, detecting all splicing aberrations published so far.

The impact of double mutations

There were 461 different double mutants inside a mutated region with 1–4 different plasmids per mutant (Supplementary Figure S3, S7 and S8). Similar to single mutations, ES was a predominant aberrant transcript. Most double mutants affecting splicing contained one mutation shown to affect splicing as a single nucleotide variant. The second mutation deepened or compensated for the effect of the first one (Figure 5a). Additionally, a few cases where the mutations affected splicing only in combination were detected, but their number was surprisingly low (Figure 5b) and the effect rather minor. To find a possible ESE hit by these mutations, the number of mutations resulting in enhanced ES was calculated for each position (Figure 5c). The highest numbers were acquired for positions 20–27, which represent a region probably binding PSF and hTra2β1 [26]. This region was also affected by single mutations, but typically the ES counted less than 1% (Figure 2).

Figure 5. The effect of chosen double mutations on splicing. A significant decrease in exon inclusion was observed in variants containing one mutation revealed among single mutations affecting splicing (A) or in variants where no mutation affected splicing (B). For comparison, appropriate single mutations were also included (grey column in A). To predict which region is responsible for predominant ES, the number of mutations in exon positions contributing to enhanced ES were counted (C). Error bars represent SD

Similar to single mutants, cryptic and de novo 5′ss use was observed. For example, double mutant 24A>G 26G>T created a strong de novo 5′ss (MaxEnt 7.92) used in 90% (median) reads or similarly, mutant 48T>G 53G>T (MaxEnt 9.72) led to complete de novo 5′ss use.

Mutations of the authentic 5′ss

One of the most important observations from the mutation set represents a diverse effect of mutations in 5′ss region on exon skipping (ES) and cryptic 5′ss use. Therefore, we investigated these mutations in more details. 5′ss mutations typically cause ES, as demonstrated also for SMN1 exon 7 e+1G>C [37]. However, Singh et al. [19] showed that this mutation can also lead to Cr1 and Cr2 use if artificial U1 snRNA binds to specific loci near authentic/cryptic 5′ss. In our minigene context, Cr1 (MaxEnt 0.76), Cr2 (4.88) and Cr3 (3.68) were used more often, which allowed us to discriminate between mutations predominantly causing ES and mutations resulting in cryptic 5′ss use. Individual mutations’ effect was very diverse, both for single mutations and double mutants (Figure 3, Supplementary Figure S9). Interestingly, cryptic 5′ss use correlated neither with the strength of the mutated authentic 5′ss, as we demonstrated for mutations non-affecting GT using MaxEnt prediction, nor with stem-loop structure formation across the 5′ss [29] (data not shown).

Based on the findings of Singh et al. [19], we speculated that the capacity of mutated 5′ss to bind U1 snRNA may be the major discriminant between ES and cryptic 5’ss use. Therefore, we simulated whether mutated 5′ss bind U1 snRNA and consequently promote distant cryptic 5′ss splicing. Analysing the built RNA duplexes along simulations revealed that the U1 snRNA could form stable duplexes with e+2T>V variant sequences of the authentic 5′ss with a similar efficiency to wt 5′ss, while binding to other mutated sequences is less probable (Figure 6, Supplementary Figure S10). Additionally, the RNA duplex stability of two double mutants (e+1G>A e+2T>C and e+1G>T e+2T>C) with U1 snRNA was tested in MD simulations. While cryptic 5′ss use is nearly equal for both mutants, only e+1G>A e+2T>C mutant shows sufficient stability in MD simulations. Thus cryptic 5′ss use cannot be fully explained just by U1 binding to the authentic splice site.

Figure 6. RNA structures of studied RNA duplexes from molecular dynamics simulations after 1 microsecond (A). The structures contain one strand from 5ʹ-mRNA and a second one from 5ʹ-U1 snRNA. Time development of Rmsd values for these systems (B). The plots reveal RNA conformation stability in the course of the simulations, i.e. small Rmsd values indicate good agreement with the initial A-RNA geometry (e.g. systems e+2U>C, A, G), larger Rmsd values indicate RNA conformation instability (e.g. e+1G>C,A,U). Instabilities in these models can be interpreted as unsuccessful duplex formation events

We hypothesised that any mutation disrupting U1 snRNA binding should result in ES, unless it is compensated by ESE creation supporting a weaker, distant 5′ss. Therefore, we analysed the possibility that an SRE was created by mutation across the original 5′ss. The most informative results were obtained using HEXplorer and ESRseq. Scores were counted for each mutation in the centre of the appropriate 21-mer, or 11-mer using HEXplorer and ESRseq, respectively. Although a single score change does not provide information about binding a specific splicing factor, the increase in both algorithms’ scores were associated with increased exon inclusion. Both algorithms showed that mutations with a substantial score increase tend to have a significantly higher cryptic 5′ss use and vice versa. Compared to HEXplorer algorithm, the ESRseq scores’ relationship with cryptic 5′ss use revealed higher linear dependency, as demonstrated by a higher coefficient of determination r [2] counted for linear regression (0.35 vs. 0.7, respectively; Figure 7). Therefore, we speculate that ESE creation over the 5′ss might have promoted the cryptic 5′ss usage.

Figure 7. Predicting SRE motifs in SMN1 exon7 variants containing intronic 5′ss mutations. HEXplorer score and ESRseq score were calculated for all single and double mutations identified in the NGS experiment and plotted against the frequency of all transcripts using cryptic 5′ss to disclose the potential dependence of cryptic 5′ss usage on splicing regulators binding across former authentic 5′ss and 3ʹ end of exon 7. HEXplorer scores and ESRseq scores were counted for these mutations’ positions for every mutation or mutation combination separately and in the case of double mutants, the average values of both scores for the position of these mutations were counted. HEXplorer score (A) and ESRseq score (B) of wt variant was subtracted from corresponding mutations scores (mut-wt) for all available mutations to reflect local SRE changes

To confirm a predominant role of SRE influence over the 5′ss in driving the Cr 5′ss use, we performed additional multiple mutant analysis. Combining specific mutations in positions e+1/2/3 with exonic mutations, the typical splicing profile was modified in both directions towards an increase or decrease in cryptic 5′ss use (Supplementary Figure S9). To elucidate SRE participation at the end of exon 7 on these shifts, we analysed the differences in HEXplorer and ESRseq score for double mutants (exonic mutation combined with e+1/2/3 mutation) and compared them to the score of the corresponding single intronic mutation and correlated them with the corresponding differences in cryptic 5′ss use (Figure 8). If the impact of e+1/2/3 mutations on RNA splicing was modulated through the exonic SREs, changes in cryptic 5′ss use would positively correlate with the differences in HEXplorer and ESRseq score. Indeed, this was what we observed. Surprisingly, opposed to the previous SRE predictions across the 5′ss, HEXplorer provided better predictive value. In fact, nearly all positive (negative) score shifts corresponded to higher (lower) cryptic 5′ss use. Sorting double mutants according to the exonic mutation position (Supplementary Figure S11) indicated that mutations in the last three nucleotides of the exon lead to lower cryptic 5′ss use and mutations in positions 48–51 to higher cryptic 5′ss usage. This data indicates the SRE change has a major influence on cryptic 5′ss use, at least in our sequence context. Thus, the transformation of 5′ss to ESE can be a major factor driving cryptic 5′ss use.

Figure 8. Modifying effect of exonic mutations on HEXplorer and ESRseq scores of intronic mutants. To elucidate how exonic mutations modify the effect of intronic mutations in e+1, e+2 and e+3 positions, HEXplorer scores (A) and ESRseq (B) scores of double mutants were compared with the appropriate single intronic mutants (double mut – single mut) and plotted against changes in cryptic 5′ss frequency in comparison with this frequency determined for appropriate single intronic mutant. The correlation shows the diversifying effect of exonic mutations on intronic mutations caused by SRE motif creation

Nucleotide exchanges as SMA determinants

To check the relevance of our results for selected single mutations (Table 1), we searched the literature for splicing-affecting mutations. Mutations causing 3′ss decreased strength typically results in ES, which is coupled with SMA outbreak [38–40]. ES abolishes the last 16 amino acids, which disrupts SMN1 oligomerization and subsequently its function [17]. Mutation e-3C>T causes slight decrease in 3′ss strength (MaxEnt score 10.42 compared to 10.92 for wt) and only mild SMA phenotype, which nicely correlates with the detected ES frequency (20%) [15]. Similarly, mutations in the 5′ss induce ES [20,37,38]. In our experiments, e+1G>C mutation caused both ES and cryptic 5′ss use, although cryptic 5′ss usage frequency was higher compared to the published values [19]. From the perspective of SMA development, any mutation inducing exclusively cryptic 5′ss use is expected to be benign, preserving an unimpaired reading frame because of the presence of natural stop codon inside exon 7. Unfortunately, no mutation in e+2 position have been described yet.

Additionally, the exonic mutation 29G>T was found to cause SMA through ES induction when using patient’s RNA [41]. It was predicted to disrupt hTra2β1 ESE. In our experiment, this mutation also showed a high level of ES, on the contrary to any other single variant in this region. Therefore, we suppose that ESS creation rather than ESE disruption is responsible for this effect. Generally, all mutations described as disease-causing in patients were detected in our experiments, while some other mutations were not described in literature. We predict that especially 23A>T and 51A>G should induce ES also invivo. The impact of mutations such as 27A>T, 48T>G or even e+2T>V is uncertain, because only ES is associated with SMA. Exon 7 inclusion regardless of the particular 5′ss used can provide sufficient SMN1 function.

Discussion

Methodological concept of parallel evaluation of splicing-affecting mutations

In this study, we present a high-throughput analysis of splicing-affecting mutations in SMN1 exon 7 and surrounding intronic nucleotides. In recent years, it was shown that such data help to improve knowledge about splicing regulation and can contribute to developing more efficient in silico prediction tools [6,7,36,42–44]. To this point, there are at least two different approaches, either a relatively short sequence in a precisely defined region is extensively varied [7,43], or the whole region is mutated in order to evaluate the effect of local mutations in the context of the wt surroundings [42,44]. We applied the second approach and specifically adopted it to analyse de novo and cryptic splice sites and even reveal mutations having a very minor effect on splicing.

For these purposes, we constructed a minigene containing ‘default’ first and last exons derived from pET plasmid and cloned SMN1 exon 7 to its multiple cloning site. This approach enables any minigene of interest to be created with no necessity for multiple mutagenesis and subcloning. The disadvantage of this approach is the absence of natural neighbouring exons, which can modify the output and thus it is suitable especially for research on a general level. The major advantage is that with several hundreds to thousands of reads per mutant, we were able to detect the differences in splicing pattern at a percentage level, as described for 6C>T. Despite this precision, the detected abundance of transcripts were, in most cases (including donor/acceptor dinucleotides) shifted toward an overrated level of normal splicing isoform. Such a shift was detected in other experiments as well [42]. We showed that these changes are predominantly caused by the presence of chimeras and suggested three possible sources for chimeric sequences. First, the products go through several PCR amplification steps in NGS library preparation and multi-template PCR is susceptible to chimera production [45]. Second, thousands of plasmids can occupy the same nucleus [46] and be subjected to homologous recombination [47] during standard transfection. And third, these chimeric sequences could originate from a trans-splicing event [48]. Regardless of the chimeric sequences’ origin, accurate transcript frequency measurement requires their exclusion from appropriate datasets. One of the possible solutions for this problem could be an additional third exon minigene barcoding.

For our experiment, we selected an SMN1 exon 7 as a well-described splicing system. Although differing in aberrant splicing isoform quantity and type from published levels, our data proved the splicing disruption for all previously described mutations. In this regard, our data offer three possible explanations for these differences: First, the above-mentioned chimeras lead to overestimated normal splicing, second, SMN1 exon 7 skipping is less frequent in HeLa cells compared to other more commonly used cell types like HEK (Figure 4), and third, using a different minigene might play a role. As an example, the detected lower effect of the 6T>C mutation (when compared to the major ES described for the SMN2 exon 7 bearing exactly this mutation) could result from the different minigene context. It was shown that the increased SMN2 exon 7 skipping rate is supported by relatively distal elements promoting ES, such as hnRNP A1 binding motifs located in the intron 7 [21], or another ISS inside intron 6 [30]. More specifically, our minigene contains the SMN1 intronic context, while it was shown that besides the 6C>T change, other SMN2-specific nucleotides in intron 7 notably participate in exon 7 exclusion [21]. In addition, it is generally accepted that the minigene context (such as the promoter, surrounding exons, etc.) influences the aberrant splicing rate.

Diverse effect of splicing-affecting mutations

Overall, only a limited number of mutations causing normal splicing disruption was detected in our experimental system and in most of them, the effect could be clearly coupled with a probable cause. Among them, 5′ss or 3′ss disruption and de novo 5′ss creation dominated. Occasionally, several single nucleotide mutations led to ES by influencing exonic SREs. ESS creation seems to be a more probable reason for changing the splicing pattern than ESE disruption, as indicated by the solitary incidence of these mutations; other substitutions in the same or any neighbouring nucleotide position did not affect the splicing pattern. For example, mutation 29G>T was predicted to disrupt hTra2β1 ESE [41]. However, other mutations in the hTra2β1 ESE, including mutations in the same position, did not result in such substantial aberrant splicing (this work, [41]). In an AG-rich region (nt 18–30), only a few single and double mutations resulted in slight ES, which showed that only more complex changes in target sequences could decrease exon recognition. This region is recognized at least by hTra2β1 [23–25], PSF [26] and their interaction partners hnRNP G [27] and SRSF9 [25]. Binding each of these splicing regulators can promote exon inclusion. Besides, other splicing regulatory proteins such as SRSF1 [22], hnRNP Q [49] or TIA-1 [31] can promote exon inclusion recruiting distinct elements. Each element in an exon containing multiple splicing enhancer elements can activate correct splicing [50] and the splicing efficiency may not be dependent on the number of ESEs [51]. Considering the partially degenerate nature of ESEs, including those binding hTra2β1 [52,53], high exon 7 occupancy with ESEs and an exons’ requirement not to have all ESEs occupied with splicing factors [54] suggests the hypothesis that a decline in binding one splicing enhancer can be compensated by binding another one or being beaten with splicing silencer.

Various mutations in the intronic part of the 5′ss resulted in a highly variable transcript profile. While mutations e+3A>T and e+3A>G caused nearly total ES, mutations in the e+2 position substantially promoted cryptic 5′ss use and mutations in e+1 and e+3A>C showed both these effects. From previous experiments, it seems that activating the distal 5′ss is negatively correlated with the strength of the proximal 5′ss [55]. However, SMN1 exon7 is not this case. We tested two major hypotheses. Based on an analysis of RNA duplex stability in molecular dynamics (MD) simulations, we tested whether U1 snRNA can bind to the mutated region of the authentic 5′ss and promote cryptic 5′ss recognition. It was shown that artificial U1 snRNAs can promote the use of both downstream and upstream 5′ss from a distance [16,19,56] and similarly, shifted endogenous U1 snRNA base-pairing enables efficient exon recognition [13]. Competition assay with mutant RNA oligos [11] showed a significantly decreased duplex stability with both e+1 and e+2 mutated positions. Among all of them, e+2U>C is the strongest competitor, but still weaker than mutants in e+3 or e+4 positions. In SMN1 5′ss the consensual exonic CAG sequence is replaced with GGA. Thus, standard WC base pair C=G in the last exonic position must be replaced with mismatched A+/C pair, which efficiently stabilizes the duplex. The stability of duplex with mismatches is probably increased by other factors including U1 snRNP components, such as U1-C [11], or additional splicing factors such as SRSF1 [57] or TIA [58]. Indeed, in our models, e+2T>V mutants showed the best duplex stability and the highest rate of cryptic 5′ss use.

In general, cryptic 5′ss activation itself was detected in many genes, but the exact mechanism has remained mostly unclear. For example, an e+1G>T mutation in the Fibrinogen β-chain gene leads to cryptic 5′ss use only if SRSF1 binds to an upstream region [59], which nicely corresponds to the ESE function in promoting exon recognition. In SMN1 exon 7, there is a very similar situation, because a mutation in the authentic 5′ss region lies upstream of cryptic 5′ss and can be accompanied with ESE creation, as supported by HEXplorer and ESRseq score alterations for e+2T>V mutants. In other mutations such as e+1G>H, U1 snRNA binding is disrupted, which releases an RNA molecule around the authentic 5′ss to the splicing factor binding and subsequently enhances the efficiency of cryptic 5′ss recognition. However, it is not clear if possible competition between the U1 snRNP and splicing factor binding can influence cryptic 5′ss recognition. Since the differences in single mutations’ splicing patterns did not enable us to distinguish between the impact of ESEs and U1 snRNA binding, we used a higher plasmid library variability generated in our NGS experiment, and analysed a series of double mutants with one mutation between positions e+1 and e+3 and a second mutation upstream. While modelling the RNA duplex’s stability of selected double mutants with U1 snRNA did not support a correlation between duplex stability and cryptic 5′ss use, HEXplorer and ESRseq score alterations, whose values are supposed to correspond with SRE creation/disruption [6,36], reflected cryptic 5′ss use nicely. Nearly all mutations in exonic nucleotides 48–51 increased both HEXplorer score and cryptic 5′ss usage, whereas mutations in exonic nucleotides 52–54 decreased both of them. The increased number of ESEs increases the probability of splicing enhancer molecule binding to the exon [54] and helps the exon to be recognized. Therefore, cryptic 5′ss use emerges from the fact that the sequence downstream of original exon 7 does not contain any normally recognizable 5′ss, and appropriate 5′ss mutants carry new additional ESEs overlapping the original 5′ss in the context of a newly defined and prolonged exon. Although this global view indicates that ESEs’ influence prevails, in single cases, U1 snRNA binding to mutated 5′ss could contribute to changes in splicing pattern. U1 snRNA binding to non-5′ss sequences has been earlier proven [60] and additionally, this binding can potentiate the use of a near 5′ss [19].

Conclusions

The use of the new high-throughput approach for detection of splicing-affecting mutations enabled us to evaluate the effect of all possible single nucleotide mutations and a lot of multiple mutations on splicing. The data indicated that the major factor preventing the exon to be correctly recognized is the creation of strong ESSs and strong de novo 5′ss. On the other hand, cryptic 5′ss activation is driven by mutations that destroy authentic 5′ss and simultaneously increase ESE strength across the mutated authentic 5′ss. U1 snRNA binding properties tune the abundance of transcripts in 5′ss mutants.

Material and methods

Gene annotation

The SMN1 sequence was acquired from the National Center for Biotechnology Information and all positions were numbered according to the transcript NM_022874.2 and protein NP_075012.1. For the purpose of this text, the numbering was simplified to reflect the position within exon 7. Thus, position 1 corresponds to the first exonic position, intronic nucleotides upstream of intron6-exon7 junction are marked ‘e-’ and intronic nucleotides downstream of exon 7 are marked ‘e+’. Nucleotide abbreviations follow IUPAC recommendations. The original exon 7 and neighbouring introns, and all corresponding nucleotides are referred to as wild type (wt).

Principles of the NGS-based method to detect splicing-affecting mutations

To develop an NGS-based method for parallel detection of different transcript isoforms in a pool of different minigene variants, a three-exonic minigene containing specific sequence elements was designed (Figure 1).

In the first step of the procedure, a plasmid library with different variants of the middle exon was created. It included two substeps. In the first one, the middle exon was mutagenised and in the second, the whole plasmid library from the first substep was subjected to intragenic barcode insertion into the first minigene exon to produce the final plasmid library. Thus, every middle exon variant was uniquely characterized by one or more specific barcodes in the first exon. In the second step, this plasmid library was sequenced by NGS to determine which barcode was assigned to the mutation. For this purpose, the whole plasmid library was cut by PacI and AsiSI and the longer fragment was re-ligated in concentrations favouring re-circularization. These shortened plasmids were subjected to PCR amplification of the minigene and subsequently sequenced. In the third step, the plasmid library was transfected into the human cells and analysed using RT-PCR followed by NGS. Therefore, RNA from these cells was extracted, reversely transcribed and amplified using primers fpet and rpet hybridized to the first and last exon of the minigene construct. Both PCR products from steps two and three were analysed in the same NGS run. This approach enabled us to directly assign any transcript (incl. ES) to the corresponding mutation in the exonic or intronic region based on a barcode in the first exon (Supplementary Figure S1).

Minigene construction

The minigene containing the SMN1 exon 7 was derived from the Exontrap plasmid pET01 (MoBiTec). The final expression cassette is depicted in Figure 1. The upstream phosphatase-derived part of the first exon was removed from the pET01 expression cassette. Additionally, AscI, NheI and BglII restriction sites were introduced into the first exon, a PacI site into the first intron 15 nucleotides downstream of the first exon, and a PciI site into the second intron 18 nt upstream of the second exon. This plasmid was used to insert SMN1 exon 7 (including 254 nt of upstream intron and 324 nt of the downstream intron) into ApaI and SmaI restriction sites within the multiple cloning site. Additionally, an AsiSI restriction site was inserted 53 nucleotides upstream of SMN1 exon 7. The nucleotide sequence of this plasmid was verified by Sanger sequencing and the potential influence of these modifications on splicing effectivity was excluded by standard minigene assay.

Plasmid library construction

Plasmid was independently mutated by site-directed mutagenesis with a series of six mutagenesis primers (Generi-Biotech): SMN1_f1: AGC TAT TTT TTT TAA CTT CCT TTA TTT TCC TTA cag ggt ttc aGA CAA AAT CAA AAA GAA GGA AGG; SMN1_f2: CTT TAT TTT CCT TAC AGG GTT TCA gac aaa atc aAA AAG AAG GAA GGT GCT CA; SMN1_f3: AGG GTT TCA GAC AAA ATC Aaa aag aag gaA GGT GCT CAC ATT CCT TA; SMN1_f4: TCA GAC AAA ATC AAA AAG AAG Gaa ggt gct cac ATT CCT TAA ATT AAG GAG TAA GT; SMN1_f5: AAG AAG GAA GGT GCT CAC att cct taa aTT AAG GAG TAA GTC TGC; SMN1_f6: GAA GGT GCT CAC ATT CCT TAA Att aag gag taA GTC TGC CAG C and their complementary reverse counterparts, where lower-case letters correspond to degenerate positions with 12% of non-original nucleotides (4% of each nucleotide) and 88% of original wt nucleotides (depicted) [61]. After E. coli transformations, about 5 000 colonies were scraped out of the Petri dishes and pooled into one homogeneous glycerol stock and stored at −80°C. The plasmid library was isolated from these bacteria with QIAprep Spin Miniprep Kit (Qiagen), according to the manufacturer’s instructions. The first exon was replaced with a barcoded one as follows (Figure 1): using similar plasmid without an AscI restriction site as a template, a fragment of the first exon was re-amplified with primers BC_f_HindIII: GCA CCT CAA GCT TCC TG and BC_r_BglII: AAA AAA AGA TCT NNN NGN NNN GNN NNA CTT TCG CGG CGG G and subsequently inserted into the HindIII and BglII restriction sites of plasmids containing SMN1 exon 7. After transformation, about 5 000 colonies were collected, homogeneous glycerol stock prepared and stored at −80°C. To check the efficiency of barcoding, plasmid DNA was isolated and subsequently digested by NheI. According to the electrophoretic gel analysis, non-barcoded plasmids should not exceed 10% of the whole plasmid library. Additionally, 20 colonies were picked from the library and independently sequenced to confirm the presence of mutations inside SMN1 exon and correct barcoding. Individual plasmids containing the chosen mutations were mutated by appropriate mutagenesis primers and standard mutagenesis procedure as described above. Four mutations (18A>G, 19A>C, 19A>T, 21A>C) were missing in the plasmid library. The four missing mutations were added to the analysis afterwards to complete the whole set of single mutations. These later added minigenes were analysed separately by classical minigene assay using capillary electrophoresis to quantify minigene transcripts. None of these mutations affected the splicing pattern compared to the wt minigene.

Cell culture and transfection

Two different cell lines were used, HeLa and HEK-293. Both lines were obtained from European Collection of Authenticated Cell Cultures (ECACC). HeLa cells were maintained in an RPMI 1640 medium (Sigma-Aldrich) supplemented with 10% FCS (Sigma-Aldrich) and HEK-293 cells in DMEM (Sigma-Aldrich) supplemented with 10% FCS. Cells were seeded into a 12-well plate one day prior to the transfection. Both cell lines were transfected using 2.4 μl XtremeGene 9 (Roche) and 400–800 ng plasmid DNA. Cells were grown for a subsequent 24 h and processed to RNA analysis. It was done either separately with one unique plasmid, and subsequently analysed by capillary electrophoresis or with the whole plasmid library (overall six independent transfections were made in separated cultivation wells in two different days), and subsequently sequenced by NGS. With the exception of the NGS experiment, all transcript analyses were assessed at least in duplicates and all corresponding data are presented as average values.

RNA extraction and RT-PCR analysis

Total RNA from cells was isolated using the innuPREP Micro RNA Kit and innuPREP DNase I Digest Kit (Analytik Jena). Potential DNA contamination was excluded by PCR using PCR primers fpet: CAG GCT TTT GTC AAC AG and rpet: GCA CTG ATC CAC GAT G and pure RNA as a template. The RNA quantity was determined by spectrophotometry (NanoDrop 2000c, Thermofisher). cDNA were synthetized using random hexamers and Transcriptor First Strand cDNA Synthesis Kit (Roche).

PCR amplification and NGS library preparation

For the purposes of the standard minigene assay (one transfection – one plasmid), PCR with rpet and FAM-labeled fpet primer was performed using recombinant Taq DNA Polymerase (Thermofisher). Transcripts were detected using capillary electrophoresis (Applied Biosystems™ 3500 Series Genetic Analyzer, Thermofisher).

For NGS-based transcript isoform analysis, minigene transcripts were amplified by fpet and rpet primers and Q5 High-Fidelity DNA Polymerase (New England Biolabs) to ensure low error rates and purified using a QIAquick PCR Purification Kit (Qiagen). To get suitable products for genotyping, the barcoded plasmid library was digested with PacI and AsiSI. Longer plasmid fragments missing most of the first intron in the expression cassette were purified from the electrophoretic gel (Qiaquick Gel extraction kit, Qiagen) and re-ligated by T4 DNA ligase (NEB) in 100 μl of the reaction mixture with plasmid concentration 1 ng/μl to re-circularize the shortened plasmids. Re-ligated plasmids were electrophoretically separated and the open-circle plasmid form was purified from the gel. The exonic parts were subsequently amplified using fpet and rint primers: ATT CAC TTT CAT AAT GCT and Platinum Pfx DNA Polymerase (Thermofisher). NGS libraries from both plasmid and transcript PCR products were prepared using KAPA HTP/LTP Library Preparation Kit (Roche) and tagged with unique 6-mer barcode sequences. These libraries were pooled and diluted with 20% PhiX and sequenced on an Illumina MiSeq.

Endogenous SMN1 expression in HeLa cells was checked after PCR amplification using primers SMN1f6: GTC CAG ATT CTC TTG ATG ATG C and SMN1r8: GTG GTG TCA TTT AGT GCT GCT C.

Computational analysis

Raw sequence data quality was checked using FastQC. Unpaired reads from each sample were independently aligned to a corresponding reference sequence using the procedure in the Rsubread R/Bioconductor package [62]. Individual reads from plasmid sequencing samples were grouped according to the barcode sequence inside the first exon, using our own scripts and variants called using GATK HaplotypeCaller version 3.3–0. For the purposes of transcript analysis, the number of individual reads were counted for each plasmid and transcript isoform separately (Supplementary Table S1). Each plasmid was represented by 227 reads on average (median) and ranged from the tenth to several hundred reads (Supplementary Figure S4). ES frequency was calculated as a ratio of ES from all identified transcript isoform reads and the counterpart value corresponds to exon inclusion (EI). Identically, using the cryptic 5′ss was calculated as a ratio of corresponding reads from all reads carrying the SMN1 exon.

To test the similarity in the expression pattern in minigenes containing the same mutation but different barcode, the Mantel-Haenszel chi-squared test was used. This analysis confirmed that in most cases the splicing pattern in these minigenes is not significantly different (Supplementary Table S1).

To test the mutations’ impact on splicing, six independent transfections were carried out. The statistical significance of the difference in ES frequency or cryptic 5′ss use frequency was evaluated using the Mantel-Haenszel chi-squared test with continuity correction and an additional filter based on a 1% threshold for changes in frequency was applied.

The splice site strength was calculated using MaxEnt [35]. To assess the effect of SREs, the following algorithms were used: SpliceAid2 [33], EX-SKIP [34], HEXplorer [36] and ESRseq score matrix [6]. HEXplorer and ESRseq scores were calculated as the difference between the given score for wt and each mutation or double mutation. Both algorithms calculate the scores for the exact mutation position, but in the case of double mutants, the score is modified as well, if the second mutation shares the same hexamers. Therefore, the scores were calculated for both positions. Typically, HEXplorer and ESRseq scores were very similar for both positions and their average value was calculated.

Molecular dynamics (MD) simulations

The interaction between mRNA and U1 snRNA was studied using MD simulations. Therefore, we built RNA duplexes consisting of one strand from mRNA and a second one from U1snRNA by the NAB module of AMBER 14 [63]. The duplexes were built as a right-handed A-RNA using a parmbsc0 OL3 force field [64]. Mismatched base pairs were built in cis-Watson-Crick cis-WC/WC geometry as standard canonical pairs. For adenine/cytosine (A/C) pair and C/C pair we used protonated nucleotide versions, i.e. A+ and C+ from the AMBER package, to create cis-WC/WC A+/C and C+/C base pairs according to the RNA base pair classification [65]. Force field parameters for pseudouridine were taken from the modified residues library [66] and adjusted according to the force field used. Each RNA duplex was solvated by a periodic octahedral box of TIP3P water molecules extending 10 Å from the solute and neutralized with sodium counterions using AMBER’s Xleap module. The equilibration and production setup used in this study corresponds to the standard protocol applied by the ABC consortium for large-scale MD nucleic acid studies [67]. Equilibrated duplexes were simulated for 1 microsecond. The MD trajectories were processed using AMBER’s Ptraj module and visualized using the VMD program [68]. Stability of RNA duplex conformations was monitored using Root mean square deviation (Rmsd), which indicates the deviation in a superposition of MD structures and the initial A-RNA geometry. Base pair step parameters characterizing the mutual orientation of base pairs in the duplex, .i.e. three rotational parameters (roll, tilt and twist) and three translational parameters (shift, slide and rise) were analysed along the trajectories using program Curves [69].

Data Availability

Sequencing data are available at Sequence Read Archive (Bioproject ID: PRJNA507345).

Acknowledgments

We thank Dr Pavel Mazura and Dr Boris Tichý for sharing their experiences with mutagenesis procedure and NGS, respectively. We also thank Lucie Kopálková for her excellent technical help. This work was supported by the Ministry of Health of the Czech Republic (16-34414A), by the Grant Agency of the Czech Republic (GA16-11619S/2016) and by Ministry of Education, Youth and Sports of the Czech Republic (MUNI/A/1298/2018). NGS data were obtained in collaboration with CF Genomics CEITEC MU and supported by the NCMG research infrastructure (LM2015091 funded by MEYS CR). Computational resources were provided by the CESNET LM2015042 and the CERIT Scientific Cloud LM2015085, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures”.

Disclosure statement

The authors declare no conflict of interest.

Additional information

Funding

This work was supported by the Grantová Agentura České Republiky [GA16-11619S/2016]; Ministerstvo Zdravotnictví České Republiky [16-34414A]; Ministerstvo Školství, Mládeže a Tělovýchovy [LM2015085]; Ministerstvo Školství, Mládeže a Tělovýchovy [LM2015091]; Ministerstvo Školství, Mládeže a Tělovýchovy [LM2015042]; Ministerstvo Školství, Mládeže a Tělovýchovy [MUNI/A/1298/2018].