In this study, we analyzed 3,568 Hungarian HBOC patients by NGS method within the frame of routine genetic testing for BRCA1 and BRCA2 genes at the Department of Molecular Genetics of the National Institute of Oncology, Budapest, Hungary between 2015–2020. In Hungary, a national guideline was published in 2020 by the Board of Clinical Geneticists about the criteria for germline testing of patients with breast cancer (http://www.hbcs.hu/uploads/jogszabaly/3278/fajlok/2020_EuK_20_szam_EMMI_szakmai_iranyelv_2.pdf). In brief, breast/ovarian cancer under the age of 50, triple-negative breast cancer or ovarian cancer or male breast cancer at any age, breast cancer at any age with two or more first-degree relatives (1 ≤ 50) or at least one ovarian first-degree cancer relatives [13]. Prior to molecular genetic testing, clinical genetic counseling was performed in each case according to the Hungarian legal and ethical regulations, where personal and familial tumor history was registered. All participants gave written informed consent for the genetic testing. Genotyping was carried out for all coding exons and exon-intron boundaries of BRCA1 and BRCA2 genes with Multiplicom amplicon-based enrichment BRCA MASTR Dx or BRCA MASTR Plus Dx library preparation kit (Agilent Technologies, Santa Clara, CA) and sequenced on MiSeq Illumina platform (Illumina, San Diego, CA). Bioinformatics analysis was done with the MASTR Reporter software v.1.1 (Agilent Technologies). When a VUS was identified, the significance of this result was extensively explained and discussed with the patients during the post-test counseling. Then, if the patient agreed to participate in research studies (including in vitro characterization and family screening for studying segregation), a second sampling was performed and these samples were further used for functional characterization. Upon conclusive result, the patients were re-counseled in the light of the new result. The availability of genetic testing was offered for all at-risk relatives of the variant carriers’ families.

Splice alteration predictions for splice consensus regions (−3 to +8 at the 5′ splice site and −12 to +2 at the 3′ splice site) were taken from ADA_score and RF_score, arising from adaptive boosting [14] and random forest ensemble [15] learning methods integrated into the annotations of dbNSFP v4.0 [16]. Cutoff scores > 0.9 for ADA and > 0.7 for RF were considered. Possible splice effects of intronic variants outside of the consensus splice regions were queried by varSEAK, an online public access program (JSI Medical Systems, Ettenheim, Germany), based on mainly MaxEntScan calculations. Scores ≥ 4 were taken as a cutoff for plausible splice impact. RNABP (http://nsclbio.jbnu.ac.kr/tools/RNABP) [17] and LaBranchoR [18] predictors were applied for determining 3′ splice branchpoint positions. LaBranchoR defines the most probable position of the active adenine and RNABP predicts the odds for a nucleotide being a potential branchpoint site. Variants were selected for cDNA-level study if any of these predictions were suggestive of possible aberrant splicing and variant frequency was extremely low (< 0.001) in various populations. Variants, elected for transcript-level analysis are listed in Table 1. Phenotype and family characteristics of the probands carrying these variants are listed in Table 2 and S1 Fig.

RNA was isolated using the Tempus Spin RNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA) from peripheral blood taken in Tempus Blood RNA tubes (Thermo Fisher Scientific) according to the manufacturer’s recommendations. First-strand reverse transcription was carried out by SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific). Reverse transcription PCRs (RT-PCRs) amplifying the variant-containing exons along with at least two adjacent exons were designed individually (list of cDNA primers is given in S2 Table). Amplification products were visualized on 1% agarose gel next to Hyper Ladder 1 kb DNA sizing standard (Bioline, London, UK) and subsequently sequenced by conventional Sanger sequencing method on ABI3130 Genetic Analyzer using the BigDye v.1.1 Kit (Thermo Fisher Scientific). Sequencing was done for the whole RT-PCR product without separation of the respective bands to compare peak intensities of normal and aberrantly spliced products. Where it was necessary to remove the interfering predominant normal alternative splice product, fragments of different sizes were cut out and cleaned from the gel by Monarch Gel Extraction Kit (New England Biolabs, Ipswich, MA) and the purified product was sequenced as above.

Relative quantitation of the normal and aberrantly spliced isoforms was assessed by two analytical methods: quantitative real-time PCR (qPCR), as well as quantitative multiplex PCR of short fluorescent fragments (QMPSF), and the averaged results of the two methods were considered for the subsequent calculations. Selective amplification of the two types of transcripts was performed with specific primer pairs engineered to discriminate between the two splice forms. At least one primer of the pairs was designed so that it should span exon borders, to amplify only from the cDNA (exact primer sequences are listed in S2 Table). qPCR was run on QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific) in relative quantification mode using SYBRGreen chemistry (Xceed HRM 2× Mix, Institute of Applied Biotechnologies, Prague, Czech Republic). Since both types of transcripts were amplified from the same template cDNA, no calibrator sample was needed, the expressions of the normal and aberrantly spliced products were directly comparable. Each test was performed in technical triplicates and three independent measurements were done. Means of the measurements with standard deviations were calculated. For QMPSF, one of the primer pairs of the respective amplicons was labeled with FAM fluorescence. PCR was conducted using Qiagen Multiplex PCR kit (Qiagen, Hilden, Germany) for 24 cycles at 62°C annealing temperature. The resulted products were subjected to capillary electrophoresis on an ABI3130 Genetic Analyzer (Thermo Fisher Scientific) in Microsatellite analysis mode. Peaks were visualized using the Peak Scanner Software 2.0 provided for the instrument. Peak ratios were calculated based on the area under the curve (AUC). Three biological measurement replicates were done.

Relative allelic expressions for allelic imbalance tests were measured on cDNA calculating the AUC ratios of exonic heterozygote positions in sequencing electropherograms [19]. The peak ratios defined on cDNA were normalized to the ratios of the same positions measured on gDNA.

Loss of the normal allele was tested in the tumor DNA of the probands, where it was available. DNA was extracted from the tumor using the Maxwell RSC DNA FFPE kit (Promega Corporation, Madison, WI). The PCR-amplicon of the variant-containing region was subjected to sequencing and allelic AUC ratio of the electropherogram peaks at the variant position was determined. LOH was calculated by normalizing the AUC ratios of the variant position of the tumor to that of the gDNA, and the tumor content was also taken into account by using the formula below:

We employed the VarSome software’s built-in pathogenicity calculator [20], corresponding with the statements of Goldgar et al. (2004) [21], for allocating variants into the 5-tier categories along with current ACMG guidelines [4]. As additional supporting evidence, at some variants we took into consideration co-segregation, LOH, family history, and proband phenotype characteristics to underpin their clinical relevance.

In the course of our routine BRCA1/2 diagnostic NGS-sequencing, we tested altogether 3,568 probands whom clinical presentation fulfilled the HBOC criteria for genetic testing (S3 Fig.) [9]. As a result of the comprehensive exon and exon-intron boundary sequencing of both genes, we detected 560 different variants, 130 of which were VUS according to the relevant ACMG criteria or not registered in BRCA1/2 locus-specific databases. The majority of them were extremely rare or absent in various population databases. We subjected these variants to diverse in silico splice prediction algorithms defining canonical splice disruptions or creation of novel splice sites (see “Materials and Methods”). Seven of the variants were suggestive for having spliceogenic effect by at least one of the in silico tools and for five variants out of them (with pan-population frequencies < 0.01 each), RNA samples were available for transcript-level analysis (Table 1). The study involved 10 nonrelated families altogether, carrying any of these five variants (Table 2, S1 Fig.).

Two probands of nonrelated breast cancer families (family 1 and family 2) carried a BRCA1 c.4484+4dupA variant, which was an insertion of an additional adenine nucleotide after the 4th basis of the BRCA1 intron 14, close to the canonical splice donor site. VarSeak prediction gave a high score for splice alteration (Table 1). Another proband in a different family was a c.4358-31A>C variant carrier (family 3). Branchpoint predictors anticipated that this latter variant affects the active adenine upstream the splice acceptor site of exon 14 (Table 1, S4 Fig.). RT-PCR amplification from the cDNA of the variant carriers with primers flanking exon 14 yielded a smaller-sized extra band in all three cases (Fig. 1A), which was sequenced and identified as an aberrant splice product with whole exon 14 skipping (Δ14) (Fig. 1B, S5 Fig.).

Semi-quantitative real-time PCR, as well as QMPSF analyses designed for specific amplification of the normal and aberrantly spliced RNA products showed that the Δ14 transcript is present in a lower quantity than the normal, full-length transcript (Fig. 1D and E). The average proportions were 0.32 (approx. ratio 1:3) for c.4484+4dupA carriers and 0.17 (approx. ratio 1:5.5) for c.4358-31A>C carrier calculated based on both detection methods. In addition, one c.4484+4dupA carrier proband carried numerous exonic heterozygote variants that allowed testing of allelic imbalance. Interestingly, AUC calculations of a heterozygote position in exon 16 (c.4837A>G) showed a 1:2 ratio of the two nucleotides in the electropherogram superposition (Fig. 1C). To resolve the discrepancy between the two ratios (1:3 and 1:2), we raised the possibility that incomplete aberrant splicing may interfere with the allele expression ratios measured. We tested this hypothesis using the c.4837A>G variant as a tagging polymorphism: we performed allele-specific RT-PCR with primers amplifying only the full-length, exon 14-containing transcript encompassing the c.4837A>G variant position. Concomitant Sanger sequencing of the resulted PCR product and measurement of the AUC of the electropherogram peaks in the tagging position yielded that the ratio was A:G=1:5.5 (Fig. 1F). Regarding the fact that Δ14-containing allele-specific PCR resulted exclusively in allele A, we could state that allele A is in “cis” position with 4484+4dupA variant, therefore the presence of allele A in the full-length transcript with 1:5.5 ratio markers that ~20% of the 4484+4dupA variant-containing allele also produced normally spliced RNA product. A remarkable portion of the remaining ~80% aberrantly spliced product may be partially degraded by nonsense-mediated decay (NMD), leaving only ~40% of the allelic expression as aberrant transcript (Fig. 1G). This indirectly suggests that NMD degrades approximately half of the BRCA1 Δ14 transcript. Unfortunately, the two other probands did not carry any exonic alterations in heterozygote form. Consequently, the calculation of the measure of incomplete aberrant splicing was feasible only indirectly in the case of the c.4358-31A>C carrier. Calculating with the same extent of NMD (~50% degradation of the BRCA1 Δ14 transcript), we yielded that nearly half of the variant-carrier allele might be normally spliced. Therefore, most probably, only half of the amount of aberrant splice product was detectable at c.4358-31A>C carrier relative to c.4484+4dupA carrier (Fig. 1G), which may be attributable to the larger incompleteness of this variant’s aberrant splicing effect.

The pedigrees of the three families are depicted in S1 Fig. All showed characteristic personal and familial HBOC features (Table 2). Family members in family 1 were available for genetic testing allowing genotype-phenotype co-segregation analysis. The proband’s mother, who was nonaffected turned out to be a non-carrier. The paternal grandmother, who had breast cancer at the age of 54, carried the variant. The proband’s father turned out to be also a variant carrier, but without clinical symptoms.

The BRCA1 c.5407-10G>A variant was detected in two independent probands in our tested cohort (family 4 and family 5) and both harbored characteristic personal and familial HBOC features (Table 2, S1 Fig). The variant changes a G nucleotide to an A in intron 22, ten nucleotides ups-tream of the exon 23, which was predicted in silico to disturb canonical splicing (Table 1). RNA sample was available from the proband of family 4. RT-PCR amplification yielded a fragment, which was indistinguishable from the wild type in length. Nonetheless, sequencing analysis of the fragment revealed aberrant splicing with retention of eight nucleotides of intron 22 upstream of the BRCA1 exon 23 (Fig. 2A). The nucleotide change introduced a novel AG acceptor dinucleotide within the AG exclusion zone [22], which acted as a novel strong acceptor site. Since the mutant transcript generated stop codon only in the last exon (exon 24), the resulting aberrant transcript was not subject to NMD. This was verified on the cDNA by a heterozygote exonic position, which actually did not show allelic imbalance (Fig. 2B). Therefore, the relative abundance of the normal and alternative transcripts reflects reliably the original ratio of the two splicing events. This was also confirmed by QMPSF technique (Fig. 2C). The completeness of the aberrant splicing was also studied applying a tagging variant c.4837A>G in exon 16, which was present in heterozygote form in one of the carriers. The tagging variant was co-amplified in a specific PCR reaction, which was designed for selective amplification of the normal, wild-type transcript. Sequencing electropherogram of the tagging variant position yielded only the G allele, no traces of the A allele (which was “in cis” with c.5407-10G>A) was detectable (Fig. 2D). This result ascertained that all the transcripts generated from the c.5407-10G>A variant-carrier allele were aberrant, so the aberrant splicing induced by this variant was complete. Furthermore, LOH test of the breast tumor tissues of index cases was available both in family 4 and family 5. Loss of the normal allele was demonstrated in both cases with R=0.26 and R=0.3 scores, respectively (Fig. 2E).

Based on in silico predictions, we selected two different BRCA2 variants for cDNA analysis, for which it was anticipated that the canonical splice sites might be affected. Of these, BRCA2 c.8487G>T, positioned in the last nucleotide of the BRCA2 exon 19, occurred in four unrelated probands in our cohort 4/3,568 (0.11%). Blood RNA sample was available from only one proband (family 6). RT-PCR–amplification of the region flanking the variant carrier exon yielded two products: one corresponded to the full-length transcript, while the other, shorter fragment proved to be an aberrant transcript with whole exon 19 skipping (Fig. 3A). The peak intensities of the sequencing electropherogram at the superposition of the normal and aberrant transcript sequences were equal (Fig. 3B); therefore, we suspected that the aberrant splicing was complete. Indeed, a tagging exonic variant position (c.7242A>G), which was also present in the proband in heterozygote state, we detected only the allele G, when amplifying the normal transcript selectively (Fig. 3C). No traces of allele A was present, implying that no full-length transcript was generated from the c.8487G>T variant carrier allele. Tumor sample DNA was available from three probands (family 6, 7, 9). LOH was demonstrated in all three samples with a mean Z score=0.25 (standard deviation, 0.03) (Fig. 3D).

BRCA2 c.793G>A affected the last nucleotide of BRCA2 exon 9 and was prognosticated as potentially spliceogenic variant affecting canonical splice donor site by varSEAK program (score 4). Opposed to this prediction, we observed neither aberrant transcript nor allelic imbalance when tested cDNA of the variant carrier (Fig. 3E and F).