RAD51D Aberrant Splicing in Breast Cancer: Identification of Splicing Regulatory Elements and Minigene-Based Evaluation of 53 DNA Variants

Simple Summary In the BRIDGES project, the breast/ovarian cancer gene RAD51D has been sequenced in >113,000 women. In the present study, we focused on the impact that 11 pre-selected RAD51D variants at the intron/exon boundaries had on the splicing process (intron removal). For this purpose, we developed a splicing reporter minigene, containing RAD51D-exons 2–9 wherein any variant could be introduced and functionally assayed for splicing alterations. All variants impaired splicing, 10 of which caused complete splicing aberrations. Moreover, we developed a minigene-based strategy to search for non-canonical, spliceogenic variants that disrupted splicing enhancers/silencers in the non-constitutive exon 3. Twenty-six BRIDGES and 16 artificial exon 3 variants were also tested. Thirty variants impaired splicing by producing variable amounts of the FL transcript. In total, up to 9 variants were classified as Likely Pathogenic, and therefore were clinically actionable. Carriers may benefit from tailored prevention protocols and therapies. Abstract RAD51D loss-of-function variants increase lifetime risk of breast and ovarian cancer. Splicing disruption is a frequent pathogenic mechanism associated with variants in susceptibility genes. Herein, we have assessed the splicing and clinical impact of splice-site and exonic splicing enhancer (ESE) variants identified through the study of ~113,000 women of the BRIDGES cohort. A RAD51D minigene with exons 2–9 was constructed in splicing vector pSAD. Eleven BRIDGES splice-site variants (selected by MaxEntScan) were introduced into the minigene by site-directed mutagenesis and tested in MCF-7 cells. The 11 variants disrupted splicing, collectively generating 25 different aberrant transcripts. All variants but one produced negligible levels (<3.4%) of the full-length (FL) transcript. In addition, ESE elements of the alternative exon 3 were mapped by testing four overlapping exonic microdeletions (≥30-bp), revealing an ESE-rich interval (c.202_235del) with critical sequences for exon 3 recognition that might have been affected by germline variants. Next, 26 BRIDGES variants and 16 artificial exon 3 single-nucleotide substitutions were also assayed. Thirty variants impaired splicing with variable amounts (0–65.1%) of the FL transcript, although only c.202G>A demonstrated a complete aberrant splicing pattern without the FL transcript. On the other hand, c.214T>C increased efficiency of exon 3 recognition, so only the FL transcript was detected (100%). In conclusion, 41 RAD51D spliceogenic variants (28 of which were from the BRIDGES cohort) were identified by minigene assays. We show that minigene-based mapping of ESEs is a powerful approach for identifying ESE hotspots and ESE-disrupting variants. Finally, we have classified nine variants as likely pathogenic according to ACMG/AMP-based guidelines, highlighting the complex relationship between splicing alterations and variant interpretation.


RAD51D
[MIM#602954] is one of the five RAD51 paralogs (RAD51B, RAD51C, RAD51D, XRCC2, XRCC3) that play an important role in the repair of DNA double-strand breaks via homologous recombination [1,2]. RAD51D loss-of-function variants confer risk of breast and/or ovarian cancer [3,4]. In a cohort of 6690 families, RAD51D pathogenic variants were associated with relative risks of 7.6 and 1.83 for tubo-ovarian cancer and breast cancer, respectively [3]. Recently, two large-scale studies have estimated an overall breast cancer relative risk of 1.8 and 1.72, respectively [5,6]. The association with triple-negative (TN) breast cancer (RR = 6.1) appears to be particularly strong [5].
During RNA splicing, introns from eukaryotic genes are removed from pre-mRNA, and consecutive exons are precisely joined together to form mature mRNA [7]. Exon recognition is a crucial step controlled by a large number of trans-acting factors, including ribonucleoproteins and splicing factors as well as cis-acting sequences identified by splicing machinery. The array of splicing signals is extensive and includes 5 and 3 splice sites, the polypyrimidine tract, the branch point and exonic and intronic splicing enhancers (ESE/ISE) and silencers (ESS/ISS) that stimulate or repress exon inclusion in mature mRNA [8]. It has been proposed that~60% of disease-causing variants disrupt pre-mRNA processing [9], generating aberrant transcripts that can affect protein function and correlate with an increased risk of a given genetic disease. Consequently, it is imperative to assess the biological and clinical significance of spliceogenic variants in order to improve genetic diagnosis and counseling as well as to trigger new therapeutic interventions [10]. Unfortunately, the lack of accurate in silico predictors makes functional assays vital to doing so. Next, splicing functional assays, either from patient or minigene RNAs, provide critical data to evaluate the pathogenicity of a particular genetic variant [11][12][13][14][15].
Our study was conducted using the framework of the BRIDGES project (Breast Cancer After Diagnostic Gene Sequencing; https://bridges-research.eu/, accessed on 1 April 2021), which has sequenced a panel of 34 known or suspected breast cancer susceptibility genes in 60,466 cases and 53,461 controls [5]. We bioinformatically analyzed 47 RAD51D variants from the intron/exon boundaries, 11 of which were selected and functionally tested in a RAD51D-splicing reporter minigene including exons 2-9 (mgR51D_ex2-9). In addition, we functionally mapped ESE-rich regions of the non-constitutive exon 3 by overlapping microdeletions and analyzed 42 candidate BRIDGES and artificial ESE variants. Finally, we made a tentative clinical classification of BRIDGES variants according to ACMG/AMP (American College of Medical Genetics and Genomics and the Association for Molecular Pathology)-based guidelines.

Materials and Methods
Ethical approval for this study was obtained from the Ethics Committee of the Spanish National Research Council-CSIC (28 May 2018).

Functional Assays
MCF-7 cell growth, transfection and inhibition of the nonsense-mediated decay w performed as previously described [16]. The triple-negative breast cancer (TNBC) cell MDA-MB-231 was cultured in Dulbecco's Modified Eagle Medium (DMEM) supp mented with 10% fetal bovine serum, 2 mM glutamine, 1% non-essential amino acids a To estimate the relative proportions of each transcript, semi-quantitative fluorescent RT-PCRs were performed using a FAM-labeled primer under standard conditions, except that 26 cycles were herein applied [14,27,28]. Fluorescent products were run with LIZ-1200 Size Standard by Macrogen (Seoul, Korea) and analyzed using Peak Scanner software V1.0 (Life Technologies). The mean peak areas of three independent experiments of each variant were used to calculate the relative proportions of each transcript as well as the standard deviations.

ACMG/AMP-Like Classification of 37 RAD51D Variants Detected in BRIDGES Samples
No RAD51D expert panel specifications of the American College of Medical Genetics and Genomics or the Association for Molecular Pathology (ACMG/AMP) variant curation guidelines were currently available (www.clinicalgenome.org, last accessed on 29 December 2020). Therefore, we performed a tentative classification based on: (i) generic ACMG/AMP guidelines [29], (ii) specific recommendations of the ClinGen Sequence Variant Interpretation Working Group for interpreting the loss-of-function PVS1 [30] and functional PS3/BS3 evidence codes [31], (iii) some non-gene-specific approaches proposed by the ClinGen CDH1 variant curation expert panel [32] and (iv) ad hoc rules based on expert judgment. The latter were essential to assign a specific functional splicing code (PS3/BS3) to complex minigene readouts (Supplementary Methods). In addition to PS3/BS3, only the rarity code (PM2) made a significant contribution to the classification process.

Bioinformatics Analysis and Functional Assays
Forty-seven RAD51D variants detected in BRIDGES subjects were collected from the intron-exon boundaries of exons 2 to 9 and analyzed using MaxEntScan (MES) or NNSplice (Table S2). Eleven potential splice-site disrupting variants were selected to be assayed in MCF-7 cells (Table 1).

ESE Mapping
Effective exon recognition in alternative splicing requires the binding of trans-acting factors to exon-splicing enhancer (ESE) sequences [34]. Exon 3 is a non-constitutive exon, excluded from several RAD51D naturally occurring, alternatively spliced isoforms [33] and consequently, it is likely regulated by splicing regulatory elements (SRE) [35]. We proceeded to functionally map SREs in exon 3 by introducing four overlapping 30-35-bp deletions (c.147_176del, c.172_206del, c.202_235del and c.231_260del) along this exon (119 bp), excluding the first two and last three nucleotides to preserve splice site-conserved positions [24,25]. All four deletions had an impact on splicing, and different levels of exon 3 skipping were observed (Figure 2A,B; Table S4; Figure S5), pointing out that these sequences might have contained regulatory motifs. Deletion c.202_235del affected splicing the most, observing no trace of the FL transcript (Figure 2A,B).
We then decided to functionally assay 25 BRIDGES variants reported at the most spliceogenic segment, c.202_235del, and surrounding sequences. In this study, variants with a FL transcript proportion ranging from 65.8 to 73.1% (<10% reduction of the FL transcript in the wt minigene) were not considered spliceogenic (see Section 4). In total, 16 out of 25 variants (64%) effectively disrupted splicing (Table 2; Figures 2A,B and S5), wherein exon 3 skipping (not observed in the wt minigene) was detected and ∆(E3_5) was increased (>17.5%) in all the spliceogenic variants. In addition, c.202G>A also expressed a new isoform, ∆(E3p36), through the use of an internal cryptic acceptor site. Remarkably, c.202G>A was the only variant that did not produce FL transcripts, whereas c.214T>C produced 100% full-length transcripts.
To refine the SRE map and identify potential hotspots for spliceogenic variants, five additional overlapping 10/11-bp deletions of the c.202_235 interval were tested. All deletions impaired splicing, being particularly significant in c.210_219del, wherein the FL transcript was absent (Table S4, Figure S5). Next, we tested all the possible changes of the highest spliceogenic interval, c.210_219. We selected the six specific nucleotides, i.e., c.212-217, within this segment that were not shared with the overlapping microdeletions c.202_211del and c.218_227del (see Figure 2A). BRIDGES variants c.213C>T, c.214T>C, c.216C>T and c.217G>A of this segment had been already tested in the preliminary SRE assays (Table 2). Next, the 14 remaining possible nucleotide substitutions in segment c.212c.217 were introduced into the minigene and evaluated. They represented the so-called "artificial" (ad hoc) variants that had not been identified in BRIDGES subjects. Eleven artificial variants within this interval impaired splicing, reducing the FL transcript by ≥10%, inducing the aberrant transcript ∆(E3) and increasing the relative proportion of the alternative isoform ∆(E3_5). Putting together the artificial and BRIDGES variants of c.212-217, the final proportion of spliceogenic variants was 66.7% (12/18).
In order to assess the accuracy of four ESE/ESS prediction tools (Human Splicing Finder (HSF), HEXplorer, Hot-Skip and ∆tESRseq), we pooled together and analyzed the 25 BRIDGES variants and the 14 artificial ones of the c.212-217 segment ( Table 2; Tables  S5 and S6) [36,37]. HEXplorer showed the highest sensitivity (74.1%, 20/27) and accuracy (61.5%) but also showed the lowest specificity (33.3%). The sensitivity of the other three predictors was very low (<37%). On the other hand, Hot-Skip presented the highest specificity (83.3%), but also 19 false negative variants. Remarkably, the spliceogenic variant c.199A>G would not have been predicted by any tool. Moreover, Hot-Skip analysis of the full exon 3 sequence suggested three changes with the highest probabilities of splicing disruption: c.163C>T (detected in BRIDGES subjects), c.163C>G and c.178C>T. We tested them in our RAD51D minigene, where they had weak-to-moderate effects, generating different quantities of the FL isoform (38.7-62.0%; Table 2).
Summarizing, we evaluated a total of 42 candidate ESE variants (26 identified in BRIDGES subjects), 30 of which impaired RAD51D exon 3 splicing.

ACMG/AMP-Like Classification of 37 RAD51D Variants Identified in the BRIDGES Cohort
We developed ACMG/AMP-based specifications for the interpretation of RAD51D spliceogenic variants (see Supplementary Methods, Figures S6-S8). With this approach, based on minigene read-outs (PS3/BS3 codes) and gnomADv2.1 global population frequencies (PM2/BS1/BA1 codes), we classified 37 BRIDGES variants ( Table 3). As expected, most variants located at canonical splice sites (9 of 11) ended up as likely pathogenic. By contrast, many of the remaining variants ended up as variants of uncertain significance (VUS), despite demonstrating an effect on splicing.

Discussion
The progress in DNA sequencing technology has resulted in the development of multi-gene panels for clinical genetic testing of diverse genetic diseases, including hereditary breast cancer. Although next-generation sequencing is a powerful tool to identify potentially deleterious variants, most variants are classified as VUS because there are no existing datasets correlating them with their impact on gene function and cancer risk [39]. The goal of the BRIDGES initiative was to establish reliable breast cancer risk estimates for 34 putative breast cancer susceptibility genes commonly included in multigene panels by sequencing more than 60,000 breast cancer cases and 53,000 controls [5]. As a partner in this project, we have carried out the most comprehensive functional evaluation of candidate splice variants of RAD51D performed so far.
Despite the fact that assays using patient-derived RNA are considered the most suitable strategy to evaluate splicing outcomes, minigene assays have proved to be a useful, simple and robust approach to assess potential spliceogenic variants [12,40]. This method presents several significant advantages such as (1) analysis of single allele events without the meddling of the WT allele; (2) accurate quantification of isoforms by inhibiting the nonsense-mediated decay; (3) high versatility, allowing the study of different variants with a single minigene; (4) high reproducibility of physiological and pathological splicing patterns; (5) analysis in a cell type relevant for the disease.
In relation to the latter, we have performed all analyses in breast cancer cell line MCF-7, and we have performed a sub-analysis in the TN breast cancer cell line MDA-MB-231, replicating findings ( Figure S3). Current data supports a shared heritability for ER-negative (or TN) breast cancer and ovarian cancer susceptibility [3,[41][42][43]. Interestingly, epidemiological and molecular evidence indicates that high-grade serous ovarian cancer arises not from ovarian epithelial cells, but from cells in the fimbriae of the fallopian tubes [44]. At any rate, the RAD51D alternative splicing profiles in fimbria, ovary and breast samples are similar [33], lending further support to the relevance of our analyses. As indicated, we have performed all experiments in malignant breast epithelial cell lines, despite the fact that the carcinogenic process may target the splicing machinery [45]. Yet, in a previous study, non-malignant cells (MCF-10A) demonstrated clear disadvantages such as very slow growth and high cell lethality after transfection [46]. Further, in the same study, we observed similar minigene splicing outcomes in malignant and non-malignant cell lines.
Many studies, including previous work by our group, support the ability of minigene assays to accurately reproduce splicing alterations as observed in RNA from carriers [16,24,25,[47][48][49]. Unfortunately, as far as we know, none of the RAD51D variants here assessed have been tested previously in RNA from carriers (the only exception being a sub-optimal study performed in a c.738+1G>A carrier, see footnote to Table 1). Yet, our RAD51D minigene (Figure 1) mimicked several physiological alternative splicing events [33], supporting the accuracy of the approach. Of note, splicing assays performed in RNA from carriers do not necessarily provide an accurate description of the actual spliceogenic effects caused by a genetic variant, as conflicting data in carriers of the same variant do exist (e.g., BRCA2 c.7976+5G>T or PALB2 c.3113+5G>C) [17,[50][51][52]. In brief, caution should be taken when interpreting the biological and clinical implications of any splicing assay, irrespective of the experimental approach. Bearing this in mind, minigenes are certainly valuable tools for initial risk assessment, particularly if RNAs from carriers are not available.
In summary, we tested 53 RAD51D variants (37 from BRIDGES subjects and 16 artificial variants) in the RAD51D minigene: 11 splice-site and 42 putative, SRE-disrupting variants. Forty-one of them (77.4%; 11 splice-site and 30 ESE/ESS variants) impaired normal splicing patterns, supporting the high frequency of this deleterious mechanism.

Splice-Site Variants
Historically, a splice variant was considered critical when it affected positions ±1, 2, but other nucleotides are crucial for exon recognition, as well. Here, we analyzed 11 changes of the intron/exon boundaries, comprising seven ± 1, 2 variants, as well as other changes at intron positions −3, −4, −8 and +6 and the antepenultimate exon nucleotide. Our results showed that 10 variants dramatically disrupted splicing because they showed no trace, or residual amounts, of the FL transcript. Five variants disrupted the 3'SS, and remarkably, three of them induced the use of de novo or cryptic acceptor sites. Six variants disrupted or weakened the 5 SS wherein exon skipping was the main outcome, although two of them also provoked the use of alternative donors. Only variant c.263+6T>C showed a significant proportion of the FL transcript (49%). Finally, it is worthy of mention that MES correctly anticipated disruptions of the canonical splice sites.

SRE-Spliceogenic Variants
Alternative splicing is a major mechanism for increasing gene expression complexity, producing multiple mRNA and protein isoforms, and is regulated by a huge array of factors [53]. RAD51D exon 3 is a non-constitutive exon [33], so it is conceivable that it is controlled by ESE/ESS elements [54]. Given the low accuracy of ESE/ESS predictors, we mapped exon 3 by overlapping exonic microdeletions. This approach revealed the existence of ESE motifs throughout this exon because all the microdeletions influenced exon 3 recognition, being particularly significant at interval c.202_235. This result confirmed that this region was characterized by a high density of cis-regulatory motifs that controlled exon 3 inclusion. Fine mapping using internal 10/12-bp deletions showed that interval c.210_219 was crucial in exon 3 recognition because its deletion strongly triggered exon skipping and did not show any trace of FL transcript. Consequently, potential spliceogenic ESE variants might have appeared in these ESE-rich segments. Altogether, we tested 42 potential ESE/ESS variants mapping to these ESE-rich intervals. Of note, in a region with high levels of naturally occurring alternative splicing and most variants having limited impacts, the definition of a spliceogenic variant is controversial. In this study, we defined an arbitrary 10% reduction cutoff for the contribution of the FL transcript to the overall expression (i.e., variants with 65.8% to 73.1% of the FL transcript were not considered spliceogenic). Based on this cut-off criteria, 30 variants (71.4%; Table 2) impaired splicing, thus supporting our SRE-screening strategy. All spliceogenic SRE-variants but one (c.202G>A) produced non-negligible levels of FL transcripts (65.1% to 27.4%), with variants c.216C>G, c.215A>C, c.200_218del and c.187G>C producing 36.4%, 28.7%, 28.7% and 27.4%, respectively. c.202G>A was the only tested SRE variant that did not generate FL transcripts. Curiously, c.214T>C presented the opposite effect, strengthening exon 3 inclusion to produce only FL transcripts.
The partial impact of these variants might have been due to the disruption of independent active ESEs, which, in physiological conditions, might synergistically cooperate for efficient exon recognition, as previously reported [27]. In effect, some deletions such as c.202_235del and c.200_219del (Table S4), which remove the crucial c.212_217 hotspot, display strong effects (no trace of the FL transcript) because they would simultaneously remove several active ESEs. In fact, any variant of the interval c.202_235 should be considered as potentially spliceogenic and functionally assayed. For instance, the recognition of atypical GC donors, which are usually associated with alternative splicing, requires the binding of several SR proteins such as SC35, Tra2β, 9G8, SF2/ASF and SRp40 [24,48,55]. Alternatively, there have been reported bi-functional SREs with enhancer and silencer properties such as the composite exonic splicing regulatory elements (CERES) [56,57], wherein the splicing impact is extremely difficult to predict. Exceptionally, one variant can decrease and increase in parallel the enhancer and silencer strengths, respectively, resulting in complete splicing anomalies. This is the case in the synonymous variant c.891C>T of the SPINK5 gene (MIM #605010), which induces exon 11 skipping. This change reduces the interaction with enhancing factor Tra2β and concurrently increases the binding of repressor hnRNPA1 [58]. Hence, RAD51D c.202G>A, which did not generate any FL transcripts, in contrast to other close variants (i.e., c.198G>T or c.208G>A), might behave like this SPINK5 variant. Likewise, the spliceogenic variants that we report here might differently alter the ESE/ESS balance such that they would induce different impacts. Nevertheless, the participation of any activator/repressor splicing factors should be demonstrated by RNA binding assays.
With regard to the ESE/ESS prediction approaches, all of them showed advantages and disadvantages (Tables 2, S5 and S6). HEXplorer showed the highest sensitivity (74.1%) and accuracy (61.5%) but the lowest specificity (33.3%); in fact, HEXplorer would have selected most variants of this critical segment (28 variants). However, the estimations of HSF, HEXplorer, Hot-Skip and ∆tESRseq would have failed in detecting 19, 7, 19 and 17 spliceogenic variants, respectively. For example, c.184T>A, which presented a moderately strong impact (62.8% of abnormal transcripts), would only have been detected by the less sensitive program (HSF, 29.6%). Tubeuf et al. [59] evaluated four ESE/ESS predictors, including ∆tESRseq and HEXplorer, in a large dataset of >1300 variants, wherein they found greater accuracy of both tools, HEXplorer specificity and ∆tESRseq sensitivity than in our study. This may be due to the small size of our sample or specific exon/gene parameters that may affect the performance of SRE predictors, as indicated by these authors. Nevertheless, the low accuracy of ESE/ESS prediction tools lays bare the need to improve the detection algorithms. The cumulative knowledge of the splicing regulation and functional tests could provide a basis for increasing the precision of SRE algorithms and the prediction of variant impacts. On the other hand, ESE functional mapping by exonic microdeletions appears as an optimal and reasonable strategy for discovering SRE hotspots [25,27]. Consequently, we propose to test all the reported variants that are located in these specific, highly spliceogenic intervals. Furthermore, the identification of spliceogenic ESE variants should be focused on exons that are targeted by prevalent alternative events.

ACMG/AMP-Based Classification of Spliceogenic Variants
Because RAD51D expert panel specifications of the ACMG/AMP guidelines were not yet available (clinicalgenome.org/, last accessed on 29 December 2020), we developed our own tentative specifications (Supplementary Methods). Upon doing that, we faced several challenges, including: (i) combining predictive and functional splicing evidence, (ii) assigning specific PS3/BS3 code strengths to complex readouts (i.e., two or more mRNA transcripts with different predicted impacts on the coding sequence), and (iii) integrating a high level of naturally occurring alternative splicing into the classification system, which is of particular relevance when evaluating the pathogenicity of VUS [60]. As far as we know, no specific ClinGen-SVI recommendations addressing these issues have been yet released [31].
We tackled them by incorporating into the classification system some ad hoc rules based on structural/functional predictions and expert judgment (e.g., when assigning a PS3 code strength to complex splicing readouts, we ignored transcripts representing <10% of the overall expression). Therefore, we do not intend to provide here (see Table 3) a definitive clinical classification of genetic variants, but rather to highlight some limitations of the current ACMG/AMP framework to interpret complex spliceogenic outcomes (as in those observed for several RAD51D variants, see Tables 1 and 2).
However, we think that, overall, our approach produces meaningful classifications, with spliceogenic variants ending up as variants of uncertain significance if expressing either a significant proportion of full-length transcripts (e.g., c.263+6T>C, c.185C>T) or in-frame transcripts of unknown impact on protein function (e.g., c.476_480+1dup, c.481-8C>A).
RAD51D c.202G>A is remarkable in that, contrary to most exon 3 variants, full-length transcripts are not produced. Yet, we have ended up classifying this variant as of uncertain significance (rather than likely pathogenic) based on: (i) the minigene readout showing a significant proportion of in-frame ∆(E3p36) transcripts, and (ii) the variant being detected in control samples (Table 3). Interestingly, RAD51D c.214T>C produced only full-length transcripts (therefore, suppressing the level of alternative splicing observed in wt alleles). The clinical relevance, if any, of suppressing naturally occurring alternative splicing, is far from obvious. At any rate, canonical transcripts code for a rare missense change (p.Tyr72His) that warrants classification as of uncertain significance. Indeed, for this and other RAD51D variants predicted to introduce missense changes and/or in-frame alterations, additional analyses are required to assess consequences for gene function.
Current ACMG/AMP guidelines [29] have not been developed to identify "intermediate risk variants". Yet, we think it is worth considering this possibility for certain variants expressing variable proportions of (likely) functional and (likely) non-functional mRNAs (see Table 3 and Supplementary Methods). It will be extremely challenging to evaluate risk for individual RAD51D variants that are exceedingly rare in the population, but it should be possible, in principle, to design burden analyses for variants displaying similar splicing effects.

Conclusions
We demonstrated through minigene-based analysis that splicing disruptions were induced by 28 out of 37 RAD51D variants (75.7%) originally identified in breast cancer cases and/or healthy controls of the BRIDGES cohort. Functional mapping by exonic microdeletions was a reliable strategy to identify SRE variants that impaired splicing. Indeed, using this approach, we identified RAD51D exon 3 sequences in which any genetic variant (regardless of its predicted coding impact) was a candidate spliceogenic variant. Finally, we have shown that current ACMG/AMP guidelines raise several issues when applied to variants associated with complex spliceogenic alterations.