Targeted Gene Panel Sequencing Unveiled New Pathogenic Mutations in Patients With Breast Cancer

Abstract

The increasing commercialization of new gene panels based on next-generation sequencing for clinical research has significantly improved our understanding of breast cancer genetics and has led to the discovery of new mutation variants. The study included 16 unselected Moroccan breast cancer patients tested with multi-gene panel (HEVA screen panel) using Illumina Miseq, followed by Sanger sequencing to validate the most relevant mutation. Mutational analysis revealed the presence of 13 mutations (11 single-nucleotide polymorphisms [SNPs] and 2 indels), and 6 of 11 identified SNPs were predicted as pathogenic. One of the 6 pathogenic mutations was c.7874G>C, a heterozygous SNP in HD-OB domain of BRCA2 gene, which led to the arginine to threonine change at codon 2625 of the protein. This work describes the first case of a patient with breast cancer harboring this pathogenic variant and analyzes its functional impact using molecular docking and molecular dynamics simulation. Further experimental investigations are needed to validate its pathogenicity and to verify its association with breast cancer.

Introduction

Breast cancer (BC) is the major cause of cancer deaths in women and the most diagnosed malignancy in the world. ¹ The incidence of BC in North Africa was higher in the continent with 29.3 per 100 000 (95% confidence interval [CI] = 20.0-38.7). ² The GLOBOCAN 2020 database of the International Agency for Research on Cancer (IARC) estimated that in Morocco BC represents 19.8 % of all cancer cases. Newly diagnosed BC cases were 38.9 % of female cancers in 2020; it was responsible for 10.5% of deaths in the country and it represents the second cause of death after lung cancer.

Breast cancer is believed to be a complicated disease that arises from a combination of environmental, microbiome, and genetics factors that increase its pathogenicity. ³ Indeed, at the genetics level, BC is associated with a number of genes with high- and moderate to low-penetrance susceptibility, implicated in DNA damage response. Breast cancer type 1 susceptibility protein (BRCA1) located on chromosome 17 and breast cancer type 2 susceptibility protein (BRCA2) located on chromosome 13 are the most important genes associated with BC, with a 60% and 55%, respectively, cumulative risks. ⁴ Mutations in other genes such as PTEN, CDH1, PALB2, STK11, and ATM were proven to confer an increased risk of BC. ⁵ Especially, PALB2 and PTEN are considered high-risk genes with a lifetime risk above 40%. On the other side, there are BC genes with low/moderate risk like NBN. ⁶ Furthermore, studies have identified a number of additional DNA repair genes that interact with BRCA1, BRCA2, and/or the BRCA pathways and confer about a 2-fold increase in BC risk, including CHEK2, ⁷ BRIP1 (BACH1), ⁸ ATM, ⁹ and PALB2. ¹⁰ Moreover, a mutation in a tumor suppressor gene, such as TP53, BRCA1, BRCA2, CHEK2, PTEN, or ATM, can be sporadic or germinal, which results in the promotion of breast carcinogenesis mechanism. Inherited mutations in certain genes such as TP53, BRCA1, and BRCA2 have been linked to an increased risk of BC. Similarly, sporadic mutations in TP53 have been detected in BC cells. These genes are crucial in maintaining genomic stability, and any mutations that affect their function can hamper the processing of DNA damage. This is likely the primary mechanism by which mutations in these 5 genes increase the risk of BC. ¹¹ Germline pathogenic variants in high-risk genes allowed to prevent the onset of cancer and to undertake risk-reduction strategies with greater efficiency, although germline pathogenic variants detected in low- or moderate-risk BC genes have also a significant clinical impact.^12,13 In addition, it has been shown that many variants may contribute to BC development, such as MSH6, RAD50, RAD51, APC, BARD1, MLH1, MSH2, MUTYH, NBN, and PMS2.^14-18 Furthermore, high-risk BC genes are those with a BC odds ratio >5.0 (BRCA1, BRCA2, PALB2, STK11, TP53, PTEN, and CDH1), according to a previous study. ⁶

With the advancement of next-generation sequencing (NGS) techniques, multi-gene panel testing is providing a cost-effective way of diagnosis, which allows the identification of genetic variants and novel mutations that have an impact on the cancer process.

In this respect, the aim of our study was to assess further potential mutations that may be related to BC by sequencing and analyzing the main 22 genes related to BC (ATM, APC, CHEK2, BARD1, BRCA1, BRCA2, MSH2, BRIP1, CDH1, PALB2, EPCAM, MLH1, MSH6, MUTYH, NBN, PMS2, PTEN, RAD50, RAD51C, RAD51D, STK11 and TP53) grouped in the HEVA screen kit. This analysis was conducted on 16 Moroccan patients with BC.

Materials and Methods

Patients

In this study, fresh breast tumor tissue was collected from 16 patients who had been diagnosed with primary BC and had undergone surgery for mastectomy at the National Institute of Oncology, Ibn Sina University Hospital of Rabat, Morocco, between June and December 2020. This study was approved by the Ethical Committee of the Faculty of Medicine and Pharmacy in Rabat, Morocco, with reference number 103/17.

DNA extraction

Genomic DNA was extracted from fresh breast tissue using QIAamp DNA Mini Kit (Qiagen GmbH, Germany). For validation of the observed mutation, DNA from blood was also isolated using the QIAamp DNA Mini Kit Blood kit (Qiagen) according to the specification of the manufacturer. DNA was quantified by Qubit fluorometer 2.0 (Thermo Fisher Scientific, Waltham, MA, USA) and its quality was assessed using 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and stored at −20°C until analyzed.

Amplicon sequencing gene panel approach

We used a commercial kit HEVA screen (4 bases SA) to analyze 22 genes associated with BC (ATM, APC, BARD1, BRCA1, BRCA2, BRIP1, CDH1, CHEK2, EPCAM, MLH1, MSH2, MSH6, MUTYH, NBN, PALB2, PMS2, PTEN, RAD50, RAD51C, RAD51D, STK11, TP53). A molecular protocol based on NGS has been used. We used 20 ng of DNA that was amplified in 3 independent multiplex polymerase chain reactions (PCRs). Then products were pooled for each sample, 1 µL of Reagent A was added, and the samples were incubated at 50°C for 10 minutes, 55°C for 10 minutes, and 60°C for 20 minutes. After that, a ligation mix was added, and samples continued in incubation at 22°C for 30 minutes and 72°C for 10 minutes. Library purification was performed using Agencourt Ampure beads (Beckman Coulter) and then eluted into a volume of 40 µL. Purified libraries underwent the insertion of adapters and indexes in the Index PCR profile and then was again purified with Agencourt Ampure beads (Beckman Coulter). The quantity and quality of prepared libraries were assessed by Qubit 2.0 fluorometer (Thermo Fisher Scientific). Each library was diluted to a concentration of 4 nM according to the manufacturer’s instructions and samples were pooled to 1 final library at equimolar ratios. A further denaturation was necessary using 0.2 N NaOH, diluted with HT1 buffer to 12 pM concentration, and 20% of PhiX Control v3 was added (Illumina). Library products were analyzed by massive parallel sequencing using MiSeq Reagent Kit v3 (600-cycle) (Illumina) on the Illumina platform MiSeq in 2 × 251 cycle profile.

Variant calling and genetic annotation

The study processed the raw paired-end FASTQ files using Trimmomatic v0.33 to filter out poor-quality reads and aligned them to the reference genome (GRCh37/hg19) using Burrows-Wheeler Aligner (BWA) v0.7.10. The subsequent steps involved removing PCR duplicates, overrepresented sequences, and low-quality reads. The identification of variants from aligned reads was carried out using the GATK software package v1.6, and coverage analysis was done in parallel with Samtools ¹⁹ . The study excluded variants such as single-nucleotide variants (SNVs) and insertions/deletions (InDels) located outside the exon region (except for splicing sites) and synonymous SNVs. Finally, the ANNOVAR software ²⁰ was used to annotate the variants, and the Integrative Genomics Viewer ²¹ (IGV; Broad Institute, Cambridge, MA, USA) was used to visualize variants against the reference genome.

Gene data analysis

The identified variants were graphically presented using OncoPrint, and Mutation Mapper was implemented in cBioportal (http://cbioportal.org). These tools provided information about genomic alterations and gave a graphical representation of the mutations found in each gene queried, survival analysis, patient-centric queries, and network visualization and analysis.

Mutation validation

Sanger sequencing was performed for tissue and blood DNA samples to confirm the presence of the R2625T mutation and to determine whether the mutation is germline or somatic. The targeted region was amplified using the following primers: BRCA2_F: GCCACCATGCTCAGCAATGAA and BRCA2_R: GTCACTGACAACTGGCTTGTG; those primers were designed using NCBI/PrimerBLAST. The PCR reaction volume was 25 µL, composed of MyTaq Mix (Bioline, BIO-25041) 12.5 µL, ddH₂O 7.5 µL, forward and reverse primers 2 µL, and template DNA 3 µL. The PCR reaction program was as follows: pre-denaturation at 95°C for 3 minutes, denaturation at 94°C for 30 seconds, annealing at 60°C for 30 seconds, extension at 72°C for 30 seconds for 30 cycles, final extension at 72°C for 10 minutes, and storage at 4°C. The purified PCR products were sequenced using the ABI3500 Genetic Analyzer (Applied Biosystems) and the BigDye X-terminator v3.1 Cycle Sequencing Kit (catalog number 4376486; Applied Biosystems) according to the supplier’s protocol. Integrative Genomics Viewer software ²² and chromatogram visualization were used to analyze and visualize the R2625T mutation.

Functional protein prediction

The functional impact of mutations was evaluated using 4 different bioinformatics algorithms: SIFT (Sorting Intolerant From Tolerant), ²³ PolyPhen-2 (PP-2), MutationTaster, ²⁴ and LTRpred. ²⁵ A SNP was considered potentially pathogenic if it was classified as deleterious by all 4 algorithms.

3D structure modeling of helical and OB1 domain of BRCA2

The 3-dimensional (3D) model of helical domain (HD) and OB1 domain (amino-acid sequence: residues 2481-2832) of BRCA2 was generated using AlphaFold software. The Mutant model was built by making point mutations R2625T on the generated AlphaFold model with the Rotamers tool of UCSF Chimera. ²⁶ Then, the mutant model was relaxed by 1000 steps of SD followed by 1000 steps of CG minimizations, keeping all atoms far by more than 5Å from the mutated residue fixed.

Molecular docking between HD-OB1 and DSS1

To gain insight into the impact of mutation on the BRCA2 HD-OB1 and deleted in split hand/split foot protein 1 (DSS1) complex, a docking study was done using HDOCK web server. ²⁷ This server supports protein-protein docking and automatically predicts their interaction through a hybrid algorithm of template-based and template-free docking. For the BRCA2 HD-OB1, the wild-type (WT) and mutated-type (MT) models generated by Alphafold were used. For DSS1, we used the Homo sapiens model extracted from 1MIU which was in complex with mouse BRCA2. The resulting complexes were analyzed by the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) method using the HawkDock server to predict the global binding free energies of the protein-protein systems and to highlight the key residues in the binding interface with the per residue MM/GBSA free energy decomposition function. ²⁸ To validate our docking protocol, 2 known pathogenic mutations located at BRCA2 HD-OB1 were used as reference. These mutations were W2626C and I2627F; the last one is known to affect the interaction between BRCA2 and DSS1 complex by reducing significantly the affinity, while W2626C has been reported to have no effect on the affinity of the complex. ²⁹

Molecular dynamics simulation

The molecular dynamics (MD) simulation was done for the WT and 3 MT (R2526T, W2626C and I2627F) helical HD-OB1 models in complex with DSS1. Molecular dynamics was carried out for 50 ns with a recording interval of 100 ps using Desmond module. ²⁸ In system builder, the OPLS3e force field was selected and TIP3P was used as a solvent model with a 10Å orthorhombic box and then the system charge was neutralized by adding 0.15 M of sodium (Na+) and chloride ions (Cl^–). The generated system was subjected to energy minimization and equilibrated via NPT ensemble at a constant temperature of 310 K and 1.01325 bar pressure. All other Desmond parameters were kept at their default values. Once MD simulation was done, simulation trajectories were analyzed using the simulation interaction diagram included in Desmond to obtain root mean square deviation (RMSD) and root mean square fluctuation (RMSF) data.

Results

In this study, we evaluated the mutational spectrum of 22 genes among 16 Moroccan BC patients using the HEVA screen panel. Among the 16 patients, 2 were filtered out due to their poor DNA sequencing quality. According to hormonal receptors status (progesterone receptor [PR], estrogen receptor [ER]) and human epidermal growth factor receptor 2 (HER2) of the remaining 14 patients, 2 were ER+/PR+/HER2+, 2 were ER+/PR+/HER2 equivocal, 6 were ER+/PR+/HER2–, 1 was ER+/HER+, and 3 were triple-negative breast cancer (TNBC) (Table 1).

Table 1.

Characteristics of patients with breast cancer stratified by age and hormone receptor status.

Patient	Age	Hormonal receptor status		HER2	Ki67 (%)
		ER	PR
1	79	+	+	–	18
2	72	+	+	Equivocal	–10
3	55	+	+	–	–1
4	53	–	–	–	NA
5	47	+	+	+	NA
6	59	–	–	–	NA
7	65	+	+	–	NA
8	37	+	+	+	NA
9	33	+	–	–	NA
10	65	+	+	–	30
11	52	–	–	–	NA
12	42	+	–	+	70
13	49	+	+	–	NA
14	41	+	+	Equivocal	NA

Abbreviations: ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; PR, progesterone receptor.

Mutational screening

A total of 11 SNPs and 2 insertions were identified from the targeted sequencing regions of 12 of 14 patients. According to the Oncoprint analysis, the SNP alterations comprised missense and splice site mutations. TP53 emerged as the gene with the highest frequency of mutations, with 4 SNPs, followed by APC with 3 SNPs, while BRCA2, CHECK2, RAD50, and RAD51C harbored 1 SNP each. For NBN and ATM genes, only 1 insertion was detected (Figure 1A).

Figure 1.

Mutational landscape of the 14 patients with BC. (A) OncoPrint showing the distribution and frequency of genomic alterations. The oncoprint was obtained using OncoPrinter (cBioPortal) for Cancer Genomic. Percentage in blue corresponds to variant allele frequency, and the letter P indicates patients’ IDs. (B) Lollipop plots showing the distribution of identified mutations in BC patients on a linear representation of the proteins and their domains. The plots were obtained by the informatics tool Mutation Mapper-cBioPortal for Cancer Genomics. BC indicates breast cancer.

Among the 13 identified variants, only 6 were known and found in different human germinal/somatic variant databases. This includes all the 4 TP53 mutations (c.528C>A, c.785G>T, c.853G>A and c.991C>T), identified, respectively, with cosmic and/or rs numbers (COSM43734, COSM11198, [rs112431538, COSM10722], and COSM11354), the c.8044G>A APC mutation (COSM4559410), and the c.7874G>C BRCA2 SNP (COSM7403374) (Table 2).

Table 2.

Nonsynonymous SNPs predicted with SIFT, PolyPhen-2, MutationTaster, and LTRPred programs.

			Mutation			Prediction
Patient ID	Gene	Gene change	Proteine change	Variant allele frequency(VAF) %	Known mutation	SIFT	PolyPhen-2	Mutation taster	LTRPred	Total deleterious
P1	CHEK2	c.580A>C	p.S194R	43		T	B	B	N	0 of 4
P2	APC	c.2965G>A	p.D989N	19.63		D	D	D	D	4 of 4
P7	APC	c.7748A>G	p.K2583R	24		T	D	D	D	3 of 4
P7	BRCA2	c.7874G>C	p.R2625T	39.64	COSM7403374	D	D	D	D	4 of 4
P8	RAD50	c.2776A>G	p.K926E	14.33		D	D	B	D	3 of 4
P9	RAD51 C	c.137T>A	p.L46H	8.51		D	D	D	D	4 of 4
P12	APC	c.8044G>A	p.V2682M	28.12	COSM4559410	D	D	D	D	4 of 4
P12	TP53	c.785G>T	p.G262 V	25.71	COSM11198	D	D	D	D	4 of 4
P13	TP53	c.853G>A	p.E285K	35.48	rs112431538 COSM10722	D	D	D	D	4 of 4

Abbreviation: SIFT, Sorting Intolerant From Tolerant.

To evaluate the impact of the 11 SNPs in terms of pathogenicity, we used, in a preliminary step, 4 functional prediction webservers (Sift, Polyphen2, MutationTaster, and LTRpred). The obtained results show that 6 mutations have been predicted as pathogenic by the 4 servers, including 2 mutations of APC (c.2965G>A and c.8044G>A), 2 of TP53 (c785G>T and c.853G>A), and 1 mutation of BRCA2 (c.7874G>C), while the c.7748A>G of APC and c.2776A>G of RAD50 were predicted as pathogenic by only 3 servers. The A580C was predicted as tolerant (Table 2). The stopgain SNPs at TP53 (c.528C>A, c.991C>T) and RAD51C (c.196A>T) as well as the stopgain insertion at ATM (7543-7544 insTCACCGACTGCCCATAGAG) are considered pathogenic as they directly affect the length of the protein. The NBN frameshift insertion (1675dupG) was also considered deleterious as it changes the reading frame. Subsequently, Mutation Mapper was used to map the identified mutations on linear proteins and their domains to determine whether these mutations affect the functional domain of proteins (Figure 1B). The results show that 6 mutations are located in functional domains, namely those of TP53 (G262V, E285K, C176*, Q331*) which are situated in the DNA-binding domain (DBD) and the P53 tetramerization domain, while the V2682M mutation of APC gene is located at the EB-1-binding domain and finally the R2625T mutation of BRCA2 is located at the helical domain (HD-OB1). This last mutation has been cited in a previous publication as a somatic bladder tumor mutation. We studied this mutation to see whether it was a somatic mutation of the breast tumor and also to see its impact on the protein. First, we validate the R2625T mutation using Sanger Sequencing, and then its impact on protein stability was assessed by molecular docking and MD.

The validation of the presence of R2625T using Sanger sequencing

To validate the presence of the R2625T mutation and to determine whether this mutation is somatic or germinal, we performed Sanger sequencing for tissue and blood DNA. The results showed that the mutation was present in both tissue and blood with 2 peaks at 7874 position of BRCA2 gene (Figure 2B and 2C). The peaks correspond to G/C alleles; this was in agreement with the Illumina result where the depth of the mutated allele was 113/285 with a percentage of 40% (Figure 2A). The presence of 2 alleles in the Sanger and Illumina results indicates that this mutation is germinal and heterozygous.

Figure 2.

Integrative Genomics Viewer (IGV) and chromatogram visualization of c.7874G>C/p.R2526T mutation of BRCA2. (A) IGV screenshot of bam file analysis of a heterozygous substitution (c.7874G>C) of the BRCA2 gene. Among the total count of 767 reads in this position, there are 452 reads with G (reference allele) and 313 reads with C (alternative allele). This distribution showed that the 7874G>C mutation is a heterozygous. The confirmation of BRCA2 mutation c.7874G>C with Sanger sequencing in tumoral tissue sample (B) and in blood sample (C) revealed that the c.7874G>C is a heterozygous germinal mutation. BRCA2 indicates breast cancer type 2 susceptibility protein.

The impact of R2625T mutation on the interaction between HD-OB1 and DSS1

The germinal R2625T mutation was predicted to be deleterious by all 4 functional prediction webservers. This mutation is located in the HD-OB1, a critical domain of BRCA2 through which it interacts with the DSS1 proteins. This led us to assume that the R2625T mutation may affect this interaction. In this regard, after the construction of the WT and MT 3D models of HD-OB1 in complex with DSS1 by molecular docking, we calculated their interaction free energy using the MM/GBSA method. To validate the used modeling/docking protocol and the obtained results, 2 known pathogenic mutations (W2626I, I2627F) of BRCA2 HD-OB1were used as reference. The calculated free energy values obtained were –300.67, –299.79, –299.31, and –280.14 for WT, R2625T, W2626C, and I2627F, respectively (Figure 3A). These results revealed that the R2625T mutation does not affect the interaction between HD-OB1 and DSS1 as there is no significant difference in the calculated free energy value between the MT R2625T (–299.79 kcal/mol) and the WT model (–300.67 kcal/mol). The calculated free energy value of I2627F (–280.14 kcal/mol) showed a significant diminution compared to WT, while for the W2626C, the free energy value was approximately equal to that of WT. It was reported in a previous experimental study that the pathogenic W2626C mutation does not affect the binding affinity between BRCA2 and DSS1, while the I2627F mutation significantly decreases this affinity. ²⁹ This was in agreement with our results and validates our docking/dynamic protocols. Furthermore, visual inspection of the interaction surface (Figure 3C) and calculation of the per residue MM/GBSA free energy decomposition showed that the R2625, W2626 and I2627 residue are not located in the main interaction surface with DSS1 and do not contribute to this interaction (Figure 3B). However, the interaction could be affected indirectly in the presence of 2627 mutation.

Figure 3.

Analysis of the interaction between DSS1 and HD-OB1 in wild and mutated types (R2625T, W2626C, I2627F). (A) Binding free energy (MM/GBSA) of DSS1 in complex with HD-OB1. (B) Per residue free energy decomposition of the binding surface from HD-OB1-DSS1 complex (wild type and mutated type). Only the top 25 residues with high MM/GBSA score were presented. The positive and negative values indicate the unfavorable and favorable contribution for the binding, respectively. (C) Localization of R2625, W2626, and I2627 on the interaction surface of DSS1 and HD-OB1. The main top 25 residues implicated in the HD-OB1-DSS1 interaction are colored in yellow. HD indicates helical domain; MM/GBSA, molecular mechanics/generalized Born surface area.

Indeed, a detailed observation of the top 25 residues of the protein-protein binding interface with the highest MM/GBSA scores (Figure 3B) showed a significant reduction in the MM/GBSA value at residue Q2655 in the mutated I2627F model; this value changed from a negative value (–2.96) in the WT to a positive value (+0.5) in the mutated I2627F model. This indicates that although residue I2627 does not have a direct impact on the interaction surface between BRCA2 and DSS1, it indirectly influences this interaction through residue Q2655.

Impact of R2625T mutation on the dynamics of HD-OB1-DSS1 complex

The results of molecular docking-MM/GBSA were further validated under dynamic conditions through MD simulation. In this respect, we performed a comparative study of MD simulation between WT and 3 MT (R2625T, W2626C and I2627F) of HD-OB1 of BRCA2 in complex with DSS1. The MD results were analyzed using (1) protein RMSD that measures the conformational changes of given complex over time and describes whether the simulation is in equilibrium (Figure 4A) and (2) protein root mean square fluctuation (P-RMSF) that characterizes local changes along the protein chain (Figure 4B).^30,31

Figure 4.

RMSD (A) and RMSF (B) analyses of C-alpha atoms for wild and mutated models (R2625T, W2626C, I2627F) of HD_OB1 in complex with DSS1 during 50 ns of MD simulations. The blue arrow indicates the position of residues 2625. HD indicates helical domain; MD, molecular dynamics; RMSD, root mean square deviation; RMSF, root mean square fluctuation; WT, wild type.

A gradual increase in the RMSD fluctuations for all systems was observed during the initial 2 ns of the simulations. Subsequently, the systems started to stabilize around 4Å after 5 ns of the simulation, with WT, R2625T, and W2626C reaching thermodynamic equilibrium after 30 ns and remaining stable until the end of the simulations with average RMSD values of 3.96 ± 0.38Å, 3.81 ± 0.41Å, and 3.75 ± 0.47Å, respectively. This indicates that these mutations do not affect the stability of the HD-OB1-DSS1 complex as the RMSD results showed no significant difference between the R2625T, W2626C, and the WY, whereas, compared to WT, the RMSD of the I2627F system continued to increase slightly until 35 ns where it reached equilibrium and showed a higher RMSD value with an average of 4.96 ± 0.8Å, indicating that this mutation could affect the HD-OB1-DSS1 complex stability (Figure 4A). On the other hand, the RMSF values do not showed any significant differences between the 4 studied systems, indicating that all these mutations do not affect neither locally nor globally the protein flexibility (Figure 4B).

Taken together, the result of molecular docking-MM/GBSA and MD showed a significant concordance, indicating that the R2625T mutation has no potential impact neither on the interaction between HD-OB1 and DSS1 nor on the stability HD-OB1-DSS1 complex. However, we believe that further in vitro analyses are needed to validate these results.

Discussion

This study describes somatic/germinal mutations detected in Moroccan BC patients using a panel of 22 genes. We identified 11 SNPs and 2 indels in TP53, APC, BRCA2, CHECK2, RAD50, RAD51C, NBN, and ATM genes. Among these SNPs, 6 were predicted as deleterious by 4 software predictors. Among these 6 deleterious mutations, 4 were found in COSMIC database (TP53 c.853G>A(E285K) and c.853G>A(E285K) mutations, APC c.8044G>A(V2682M) mutation, and c.7874G>C(R2625T) BRCA2 mutation). The TP53 c.853G>A mutation is a germinal deleterious mutation that has been reported in previous studies; it is located within the DBD of TP53 and it is a temperature-sensitive mutation inducing the loss of its function at 37°C. ³² E285K was shown to be associated with shorter survival in advanced non–small cell lung cancer (NSCLC) treated with tyrosine kinase inhibitors (TKIs) ³³ ; furthermore, Blanden et al demonstrated in their study that E285K was the most unstable variant of TP53 and resistant to the zinc metallochaperone ZMC1. ³⁴ In this study, we investigated one of the deleterious SNPs, c.7874G>C, that causes the change of arginine (R) to threonine (T) at position 2625 of the BRCA2 protein. This mutation has only been identified in a single study in a patient with bladder cancer and has been described/confirmed as a somatic mutation. ³⁵ In contrast, in our study, we found that this mutation is germinal and heterozygous using Sanger sequencing for both tissue and blood DNA.

The R2625T mutation is located at the helical-binding domain of the DBD of BRCA2 protein, specifically at OB1 motif. It was demonstrated that mutations in this domain disrupt the interaction between BRCA2 and DSS1, a small acidic protein required for the stability of BRCA2 which interrupts homologous DNA repair via BRCA2 in cells with DBD mutations, thus causing cancer evolution.^29,36 These findings led us to assume that the mutation R2625T might decrease DSS1-binding capacity with BRCA2 and inhibit its function. Molecular docking between DSS1 and WT/R2625T BRCA2 indicates that this mutation had no impact on HD-OB1-DSS1 interaction. Moreover, MD also showed that this mutation had no impact on HD-OB1-DSS1 complex stability. Indeed, the RMSD results showed no significant difference between the R2625T and the WT, with average RMSD values of 3.96 ± 0.38 and 3.75 ± 0.47, respectively.

Molecular docking and MD results were referenced to 2 known pathogenic mutations W2626C and I2627F. Experimental data from these mutations showed that W2626C does not affect or slightly affects the interaction between DSS1 and BRCA2, while I2627F highly affects this interaction. ²⁹ Using our molecular docking and MD protocol, we found comparable results to those obtained from experimental data, which validates our docking and dynamics protocol and provides a baseline for comparison with R2625T results.

Although the results for R2625T showed that this mutation has no impact on the interaction between DSS1 and BRCA2, we believe that further experimental analysis is are needed, particularly by studying the impact of this mutation on the dislocation of BRCA2-DSS1 complex from nucleus to cytosol. Indeed, this phenomenon of dislocation alteration is observed especially in the presence of some mutations, that even retain residual DSS1-binding capacity have similar affect to that, ablate it more severely in impairing the nuclear localization of full-length BRCA2 protein. ²⁹ For example, the W2626C mutation, which is adjacent to R2625T, retains significant DSS1 binding and also forms intracellular oligomers; however, its nuclear localization is impaired, which causes a lot of homologous DNA recombination. ²⁹

Collectively, this study provides an insight on the germinal/somatic mutational profile of Moroccan BC patients, which will help in the optimization of specific gene panel testing for Moroccan population. It is necessary to increase the number of the population size to have a global view. Furthermore, the in silico techniques used in this study could be used as proof-of-principle in the preliminary step for studying the impact of HD-OB1 mutations on HD-OB1-DSS1 complex stability.

Conclusion

This preliminary study aims to assess further potential germline /somatic mutations that may be related to BC using targeted gene panel sequencing. Collectively, the results of this study revealed the presence of 11 SNPs in various genes, such as TP53, APC, BRCA2, CHECK2, RAD50, RAD51C, NBN, and ATM. It is noteworthy that 6 of these SNPs were predicted to be deleterious. R2625T mutation was thoroughly investigated in this study as it is located in the BRCA2 HD-OB1 domain, a critical domain through which BRCA2 interacts with a stabilizing protein, the DSS1. Docking and MD demonstrated that this mutation had no impact on the affinity and stability of the HD-OB1-DSS1 complex. An in-depth experimental study is necessary to better understand the pathogenic mechanism of this mutation.