Targeting BRF2 in Cancer with Repurposed Drugs Using In Sil- ico Methods

Overexpression of BRF2, the essential subunit of RNA polymerase III, is required for the development of a variety of cancers, including lung squamous cell carcinoma and breast cancer. BRF2 also acts as a central redox-sensing transcription factor (TF) and is involved in the regulation of oxidative stress (OS) pathway where its amplification enables cancer cells to evade OS-induced apoptosis. Here, we experimentally confirmed a link between BRF2 and DNA damage response (DDR) signaling. Focusing on its DNA binding capacity, we targeted BRF2-TATA binding Protein (TBP)-DNA complex interaction to functionally inhibit this transcription factor (TF). We found that BRF2 binding to TBP changes the conformation of this protein and therefore characterized these events in detail. Virtual screening allowed selection of the optimal drug based on binding energy, and intermolecular, internal, and torsional energy parameters. Molecular dynamics simulation confirmed the docking results, and all simulations were validated by root-mean-square deviation (RMSD). According to the simulation, the chosen drug was efficacious in targeting BRF2 which therefore requires further experimental validation to investigate the effect of this drug on BRF2 functions and its downstream pathways.


Introduction
RNA polymerase III (Pol III) is critically important for transcription initiation. TFIIBrelated factor 2 (BRF2) is a Pol III complex subunit required as a master regulator of oxidative stress, where its overexpression is linked to tumorigenesis [1]. Studies have shown high basal levels of BRF2 expression in breast and lung cancers [2,3], where it is a potential independent prognostic factor for the recurrence and metastasis of lung cancer. BRF2 was identified as a lineage-specific marker in lung cancer [4] and was later recognized as an oncogenic driver in breast cancer [5] and esophageal carcinoma [6]. The genetic activation of BRF2 has also been linked to lung squamous cell carcinoma and other cancers [7,8]. BRF2 knockdown can inhibit the migratory and invasive abilities of Non-Small Cell Lung Cancer cells (NSCLC) and induce loss of epithelial-mesenchymal transition [9,10]. Additionally, BRF2 has been shown to play a crucial role in breast cancers with heterogeneous HER2 gene amplification. Within such tumors, only 10 to 15% of the cells show amplification of a HER2 gene signature [3], whereas the remaining population is HER2-negative. By analyzing multiple TCGA cancer datasets, BRF2 has been shown to be significantly upregulated in a HER2-negative compartment of these tumors to compensate for the lack of HER2 amplification [3]. Multiple cancer models have also shown significant alterations in BRF2 copy number variations, specifically copy number gain and amplification [2]. Targeting specific Pol III components, such as BRF2 and other related molecular subunits, could provide some advantages due to its specific role in regulating certain aspects of transcription and/or the oxidative stress response. We and others have shown a contextdependent role of BRF2 as a cis-associated 'alternative driver oncogene' with a significant role in lung and ER-/HER2+ breast cancers, respectively [2,4].
BRF2 is a component of the Pol II I complex, which transcribes 5S rRNA, transfer RNA (tRNA), and other small non-coding RNAs (ncRNAs), such as small nuclear RNAs (snRNAs), and yet other non-coding RNAs collectively identified as type III genes, which are generally less than 400 base pairs long [11]. Pol III genes' encode essential house-keeping genes that regulate protein synthesis, hence cell proliferation and growth [12]. Since cancer cells require an increased protein biosynthesis for rapid proliferation, Pol III plays a central role in cancer progression to maintain a high rate of transcription, thus sustaining transcriptional machinery requirements. All Pol III subunits function together with three Transcription Factors (TFs): TFIIIA, TFIIIB and TFIIIC. While TFIIIA is a single protein that specifically recruits Type 1 promoter (5s rRNA), TFIIIB consists of three components: TBP (TATA binding protein), BDP1 and either of the TFIIB-related factors BRF1 (for tRNA) or BRF2 (only in type 3 or U6 promoter). The functions of TFIIIs are to recruit Pol III to the transcription start sites through a series of DNA-protein and protein-protein interactions. In order to target BRF2 effectively, the structure, interactions, functions, and regulatory role of this TF, among other subunits and components of the RNA polymerase III (Pol III) machinery, should be taken into consideration. BRF2 recruits RNA Pol III to type III gene external promoters, including the gene product selenocysteine (SeCys) tRNA, which is involved in the synthesis of selenoproteins, and acts to reduce reactive oxygen species (ROS) to maintain cellular redox homeostasis [1]. Absence or defective expression of selenoproteins induces apoptotic cell death [1]. Moreover, high baseline levels of BRF2 have been observed in cancer cell lines that are subject to prolonged oxidative stress [1]. Therefore, it is speculated that high levels of BRF2 expression in cancer cell lines supports sufficient expression of selenoproteins, in order to detoxify ROS and to preserve redox homeostasis, which is a hallmark of cancer cells.
In the current study, we have identified a novel function for BRF2 in regulating the DNA damage response, independent of its role in the oxidative stress response. Furthermore, we have targeted the BRF2-TBP-DNA complex with FDA approved/under clinical trial investigational drugs to inhibit its transcriptional activity and/or DNA binding to interrupt the functionality of this complex following by a detailed analysis of molecular events that disrupt BRF2-TBP complex. The ploidy/DNA content; 2 is the equivalent of 2n, the normal/diploid amount, whereas 4 would occur after whole-genome doubling. The larger the deviation from 2, the more aneuploid the cells are. Ploidy was determined using the ASCAT (allele-specific copy number analysis of tumors) algorithm, as per Van Loo et al., 2010 (PMID 20837533). (H) ploidy_groups: In this panel, the ploidy data shown in panel (g) are binned, with "diploid" referring to tumor samples with 2n DNA content, "aneup_lo" referring to tumor samples that are aneuploid but near diploid/2n (either lower or slightly higher than 2n DNA content), "aneup_high" referring to tumor samples that are aneuploid and with DNA content well above 2n, and "aneup_all" referring to all aneuploid tumor samples (sum of former two bins). (I) cal_scna_burden: The number of chromosome armlevel somatic copy number alterations. This is the total number of chromosome arms that are gained or lost, also known as chromosome arm aneuploidies (CAA). This was determined using SNP6 array copy number data, as per Shukla et al., 2020 (PMID 31974379). (J) wc_aneuploidy: The total number of whole chromosomes that are gained or lost. This was determined using SNP6 array copy number data, as per Thangavelu et al.

TCGA Analyses for Cellular Signaling and Immunoregulation
The Spearman correlation analysis between the gene expression levels (mRNA levels from The Cancer Genome Atlas (TCGA) RNAseq datasets) and the variables listed below based on the dataset available in listed publication: (A) Pathway scores for each of these twelve pathways were determined by integration of multi-platform data from several types of 'omics' analyses, as described by Hoadley et al., 2014 (PMID 25109877

Virtual Screening Procedure
A drug library, containing 8.770 compounds in sdf format, has been extracted from the DrugBank database v. 5.1.5 [41] (https://www.drugbank.ca/). BRF2 receptor and drug structure files have been converted into pdbqt format using the prepare_receptor4.py and the prepare_ligand4.py tools of the AutoDockTools4 software [42,43]. Peptide drug structures and ligand containing Si atoms due to the absence of force-field parameters have been excluded from the analysis, The 8000 molecular docking simulations have been performed using an Autodock Vina 1.1.2 program [44]. Each drug was run through one docking simulation which included ten docking runs. Generated by the Lamarckian Genetic Algorithm, 10 binding poses representative of the cluster centroids of all conformations are selected by the the AutoDock Vina program. A box of size x = 66.98 Å, y = 46.90 Å, z = 44.83 Å has been placed over the full structure of BRF2 structure.

Molecular Dynamics Simulations for Protein-Ligand Complexes
Three-dimensional structure of BRF2 was extracted from the RCSB Protein Data Bank (PDB ID 4ROC), which is in complex with TBP and DNA molecules. Two sets of Molecular Dynamics Simulations were performed. In the first MD simulation, we investigate the stable conformation of ligand-free BRF2 structure before binding to TBP and DNA molecules. MD simulation and molecular mechanics (MM) minimization were performed using Gromacs 2019.6 package under an amber ff99SB forcefield. MD simulations were carried out with periodic boundary conditions. Van der Waals forces and the Particle-mesh Ewald method were both treated with a cut-off of 10 Å. Also, the frequency to update the neighbor list was 5. The protonation state of Gromacs package was used to calculate the total charge of protein. The protein was solvated with TIP3P water molecules with an 8 Å radius buffer zone around the protein in a truncated octahedron periodic box. The systems were neutralized by adding the corresponding number of counterions (Na + and CL -) using the genion module. MD simulation was accomplished in four steps. In the first step, the entire system was minimized using the steepest descent followed by conjugate gradients algorithms. In the second step or equilibration step, heavy atoms were restrained using a force constant of 1000 kJ/mol nm and the solvent and ions were allowed to evolve. This was done through minimization and molecular dynamics in the NVT ensemble for 200 ps and in the NPT ensemble for 200 ps. Equilibrium geometry was required to be at 298 K and 1atm. To achieve this, the the temperature of the system was increased and the velocities at each step were reassigned according to the Maxwell-Boltzmann distribution at that temperature and equilibrated for 200 ps. Temperature coupling was set to 0.1 ps and pressure coupling to 2 ps. The V-rescale and Parrinello-Rahman algorithms were used for the thermostat and barostat during the equilibration step, respectively. All bonds were constrained via the LINCS algorithm. In the final step or production phase, 100 ns MD simulation was performed under an NPT ensemble. Nosé-Hoover thermostat and Parrinello-Rahman barostat were used with removing position restraints to retain temperature and pressure stable in the production step. The temperature was at 310 K with a time step of 2 fs. Constraining the lengths of hydrogen-containing bonds was further improved by addition of the LINCS algorithm.
Then the output of the first MD simulation was used as input for virtual screening. One of the best-ranked complexes from virtual screening was minimized and used for second molecular dynamics simulations. To parametrize ligand molecule, GAFF forcefield which was assigned by the Antechamber program in AMBER tools was used. AM1-BCC charge model was used to calculate the atomic point. The coordinate file, in SYBYL mol2 format, was formerly loaded to Antechamber program in order to set these parameters. After generating AMBER topology and coordinate files, these files were converted to GROMACS topology and coordinate files using acpype conversion script. Other conditions are similar to first MD simulation.

Binding Free Energy Calculations
The molecular mechanics Poisson Boltzmann surface area (MM/PBSA) method is the widely used method for binding free energy calculations from the snapshots of MD trajectory. The binding free energies of the complexes between BRF2 and ligand were analyzed during equilibrium phase by taking 400 snapshots from 60 to 100 ns MD simulations, using g_mmpbsa tool of Gromacs [45].
The binding free energy of ligand-protein complex in solvent was expressed as: where Gcomplex is the total free energy of the protein-ligand complex, Gprotein and Gligand are total energy of separated protein and ligand in solvent, respectively. The free energy for each individual Gcomplex, Gprotein and Gligand were estimated by: where x is the protein, ligand, or complex. EMM is the average molecular mechanics potential energy in vacuum and Gsolvation is free energy of solvation. The molecular mechanics potential energy was calculated in vacuum as following: where Ebonded is bonded interaction including bond, angle, dihedral and improper interactions, and Enon-bonded is non-bonded interactions consisting of van der Waals (Evdw) and electrostatic (Eelec) interactions. ΔEbonded is always taken as zero.
The solvation free energy (Gsolvation) was estimated as the sum of electrostatic solvation free energy (Gpolar) and apolar solvation free energy (Gnon-polar): where Gpolar was computed using the Poisson-Boltzmann (PB) equation and Gnon-polar estimated from the solvent-accessible surface area (SASA) as equation following: where γ is a coefficient related to surface tension of the solvent and b is fitting parameter. The values of the contats are as follows: γ = 0.02267 Kj/Mol/Å 2 or 0.0054 Kcal/Mol/Å 2 b = 3.849 Kj/Mol or 0.916 Kcal/Mol

Umbrella Sampling (US) Simulation
The average structure from the last 10 ns of the simulation was selected for further US study. The initial system was prepared by using Gromacs2019.6 software The proteinligand complex was first made parallel to the y-axis. To pull the ligand along the y-axis by a 2.5 nm distance, a box was constructed with a y-axis length of 12 nm. The prepared box was solvated, neutralized, minimized and equilibrated at a temperature (NVT) and specific pressure (NPT) similar to the previous MD simulation step. The US simulation began with the center-of-mass-pulling method. The ligand was pulled from the protein pocket towards the solvent bulk over the course of 500 ps by using a 600 kJ/(mol*nm) force. The ligand pulled at the rate of 0.005 nm per ps. During this simulation, snapshots were saved at each ps, so in total 500 configurations were generated from the pulling simulations. Several configurations were used as starting configurations for each US simulation, where each was independently simulated by performing NPT equilibration for 100 ps. Next, 5 ns MDrun was performed for each individual configuration. Using the outputs from this US, the potential mean force (PMF) was calculated using the weighted histogram analysis method (WHAM), included in Gromacs as the gmx_wham utility. The force (kcal/mol) needed to pull the ligand from the binding pocket and the corresponding distance pulled are demonstrated on the PMF graphs (y and x axes, respectively). To attain equilibration, we used in total 27 configurations. The binding free energy (DG) was calculated for this ligand by taking the difference between the plateau region of the PMF curve and the energy minimum of the simulation.

Principal Component Analysis (PCA)
PCA was performed as described previously [45], to attain a mass-weighted covariance matrix of the protein atom displacement. This parameter is indicative of dominant and collective modes of the protein from the total dynamics of the MD trajectory. A set of eigenvectors and eigenvalues reflecting the molecules concerted motion were extracted by diagonalizing the covariance matrix. The g_covar which is the Gromacs in-built tool was used to yield the eigenvalues and eigenvectors by calculating and diagonalizing the covariance matrix. The g_anaeig tool was used to analyze and plot the eigenvectors. The gmx sham was used to calculated free energy landscape based on PC1 and PC2 eigenvectors.

Cell Culture
All cells were cultures based on ATTC instructions. Cultures were passaged every 4 days as per manufacturer's instructions (Stem Cell Technologies). All the cell lines were routinely tested for Mycoplasma infection by scientific services at QIMR Berghofer Medical Research Institute. or PEI (Polysciences, Inc., 23966-2, POL) (1:3 ratio) in Optimum media. At 5 h post-transfection, media was changed and the packaging cells incubated for 72 h. At 48 h and 72 h post-transfection, the media was filtered using a 0.45 µm filter onto target cells prior to spinfection (at 1000 X g for 1 h at 25° C) in the presence of hexadimethrine bromide (polybrene; Sigma Aldrich®, H9268-5G). Media was removed after spinfection and cells were allowed to recover for the next 48 h before selection with the corresponding antibiotic was carried out to select for transduced cells. Antibiotic selection was sustained until an untransduced control plate of cells had all died. Sequence:CCGGCCTCAGGAAGTTAGGGACTTTCTCGA-GAAAGTCCCTAACTTCCTGAGGTTTTT

Western Blot and Immunoblotting
Protein lysate was made by Urea/SDS buffer lysis buffer following by sonication and quantification of protein concentration using BCA assay kit (Thermo Scientific) according to manufacture protocol. For WB, 30 µg of protein was loaded in each well. The samples were incubated for 5 min at 95 ºC with 5X Laemmli loading buffer. Afterward, the protein lysates were run on the gel by 120 V. After a successful separation of the protein, the transfer of the protein on the membrane was performed. For the transfer, nitrocellulose membrane was used following by staining in Ponceau S to confirm the transferring.
For the immunodetection the membrane was incubated for 1 h in blocking buffer, following with overnight incubation (4 °C) with the primary antibody. After 3 times washing with PBS/Tween followed by 1 h incubation with peroxidase-conjugated secondary antibody and again 3 times washing, the immunodetection was undertaken. The detection of signals was based on the chemiluminescence reaction between peroxidase and luminol/hydrogen peroxidase mixture from Western Lightning ® Plus ECL following by development by ChemiDock Imaging System (BioRad). Table   Primary

Cisplatin Treatment
MCF10A cells (10 6 ) were treated, after overnight seeding in a 10 cm Petri dish, with 10 µM cisplatin. After different time points the cells were collected by using the cell scraper. The media was removed and PBS/vanadate buffer (2 ml) was added to the cells before the mechanical detachment began. Afterwards, the centrifugation (5 min, 4 °C, 1400 rpm) of the cell suspension, the supernatant was discarded and the pellet was resuspended in lysis buffer (Urea/SDS buffer).

Chromatin Extraction
The chromatin extraction was performed after the cells were treated with an irradiation dose of 6 Gy using an irradiator machine. The cell pellet was resuspended in 3 -5 times volumes of lysis buffer and incubated for 20 min on ice. Following centrifugation of 14000 rpm at 4 °C for 10 min, the supernatant, which contained the cytoplasmatic proteins, and the pellet, which contained the nuclei, were collected. The nuclei were washed once with lysis buffer and centrifuged (14000 rpm, 4 °C, 10 min). In order to receive nucleoplasmic proteins, the nuclei were mixed with 2 -5 volumes of low salt buffer + 1 % Triton X 100 and were incubated for 15 min on ice. Subsequently, the suspension was centrifuged (14000 rpm, 4 °C for 10 min), and the supernatant included the nucleoplasmatic proteins and the pellet the chromatin. After collecting the nucleoplasmic proteins, the pellet was resuspended with 2 -5 volumes of HCl 0.2 N and incubated for 20 min on ice. Following centrifugation (14000 rpm, 4 °C for 10 min), the supernatant, which comprised the acidsoluble proteins which is the chromatin fraction, was kept and neutralized with the same volume of 1 M Tris-HCl pH 8.

Statistical Analyses
All statistical analysis was performed using GraphPad Prism v 8.0, using a general linear statistical model as defined in each section. Error bars in all figures represent mean + standard error of the mean (SD). P values of < 0.05 were considered statistically significant. Statistical significance is designated with an asterisk (*), where *p < 0.05, **p < 0.01, ***p < 0.001 and ****p < 0.0001.

Bioinformatics Analysis Identifies BRF2 as a Promising Target in Cancer
BRF2, which localizes to the 8p11.23 chromosomal region, was initially identified as a lineage-specific marker in lung cancer [13] and later in HER2-positive breast tumors [3]. By performing copy number-altered network analysis, we identified BRF2 as a potential cis-associated driver gene in breast cancer [2]. Further examination of the literature revealed that aberrations in this region have been noted in breast tumors and these correlate with metastatic relapse within one year [14,15]. To further validate this, we performed pan-cancer analysis on BRF2 expression levels and copy numbers using TCGA datasets. Comparison of RNA Seq data of gene expression levels in normal and tumor tissues in 32 cancer types reveals that BRF2 is significantly overexpressed in many tumors, including malignancies of the breast, lung, liver and kidney (Fig 1A). In many cancers, such as esophageal squamous cell carcinoma, adrenocortical carcinoma, invasive breast carcinoma, NSCLC, bladder carcinoma and sarcoma, BRF2 is often amplified and/or gained through copy number alterations (Sup fig 1A). Analysis of TCGA data in lung cancer depicts amplification/copy number gain in lung cancer, squamous cell carcinoma, and adenocarcinoma by 32%, 41%, 23%, respectively (Sup fig 1B-D). Verification of BRF2 against an independent dataset of 597 PAM50-subtyped breast tumors showed that BRF2 was upregulated, gained, or highly amplified in 37% of cases of TCGA, with copy-number amplification strongly correlating to mRNA upregulation in a substantial subset within each of the PAM50 subtypes (Sup fig 2E)[16]. Furthermore, METABRIC dataset shows 14% amplifications in BC across 2509 patients (Sup fig 2F). Moreover, analysis of  2G). Correlation with survival data from these patients showed a significant association of BRF2 upregulation with poor survival, particularly in ER-positive tumors (Sup fig 2 H,I). Notably, we found that BRF2 amplification was mutually exclusive to loss of BRCA1 and BRCA2 in tumors both in TCGA and METABRIC data sets (Sup fig 2E,F). BRCA1 and BRCA2-associated tumors are defective in homology directed DNA double strand break (DSB) repair, BRF2 amplification is rarely seen in the context of DNA damage repair deficiency. Thus, it is possible that BRF2 amplification in context of BRCA-loss is synthetic lethal interaction that causes accumulation of DNA damage and consequent cell death. This analysis suggests that BRF2 is a potential therapeutic target in multiple cancers, including breast and lung cancer.  Next, we analyzed the correlation between BRF2 expression and 18 measures of genomic instability, a hallmark of cancer. These analyses revealed that there are no strong correlations between the expression of BRF2 in breast cancer dataset and selected features of genomic instability (Sup fig 2A-R).

Evaluation of the Role of BRF2 in the Regulation of DNA Damage Response (DDR) Pathway Utilizing Normal Mammary Epithelial Cells and Breast Cancer Lines
As BRF2 senses oxidative stress[1], and excess ROS causes severe damage to cellular macromolecules, especially proteins and DNA, we aimed to investigate its role in regulation of DDR in more detail. Upon exposure to oxidative agents, BRF2-depleted cancer cells undergo oxidative stress-induced cell death, but normal cells avoid this through upregulation of selenoproteins. ROS is known to activate DNA damage response signaling via induction of DNA damage. We therefore suspect that a relationship exists between DDR signaling and BRF2. Our data suggests that BRF2 is markedly increased in expression upon exposure of MCF10A cells, near normal mammary epithelial cell line derived from a patient with fibrosarcoma, to different dosages of γ-irradiation (IR), harvested cells 1 hr after radiation exposure ( Fig. 2A). γH2AX, the biomarker for DSBs and p53-S15 and pKap1 (S824) phosphorylation, a downstream substrate of ATM/ATR, indicated the level of DNA damage induction in the cells. We also detected a robust increase in BRF2 (within 1 hour after exposure to 0.5 Gy IR), which was more efficient than other markers of DDR signaling, including γH2AX and p-Kap1 ( Fig. 2A). Further, a time course experiment after exposing cells to 6 Gy IR revealed that BRF2 induction was very rapid and persisted during the experimental time course of 6 hrs in both MCF10A and MDA-MB-231 cells (Fig. 2  B,C). In another breast cancer line, SUM159PT, using a longer time course of up to 24h, we found similar rapid and persistent upregulation of BRF2 upon DNA damage (Fig 2D). This increase is not mediated due to an increase in mRNA levels of BRF2 (Fig 2E, data shown for MCF10A). Additionally, we used the chemotherapeutic drug cisplatin as another DNA damage-inducing treatment to exclude the possibility of induction of BRF2 by oxidative stress, which ionizing radiation is known to cause. BRF2 demonstrates a clear upregulation following cisplatin treatment (Fig. 2F), which implies a link between BRF2 and DNA damage induction. Notably, after cisplatin treatment BRF2 stabilization is independent of oxidative stress, as a slight decrease, rather than increase, in nuclear factor-E2related factor-2, NRF2 levels is evident under the same experimental conditions. NRF2 is Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 December 2020 doi:10.20944/preprints202012.0770.v1 a transcription factor that regulates the expression of antioxidant proteins and protects against oxidative damage [17]. This strengthened the possibility of the link between BRF2 upregulation and DNA damage induction, without the confounding effect of oxidative stress. Moreover, we found a robust stabilization of BRF2 that correlated with an increase in γH2AX in MCF10A cells treated with Tert-butylhydroquinone (tBHQ), the major metabolite of butylated hydroxyanisole, which induces an antioxidant response through the redox-sensitive NRF2 transcription factor (Fig. 2G). This experiment further links BRF2 to DDR through the induction of oxidative stress. Additionally, we found that knocking down of BRF2 using siRNA (MCF10A cell) causes an increase in baseline and IR-induced DNA damage, which was assessed by immunoblotting for γH2AX, however no difference was observed in IR-induced induction of pKAP1 (S824), a direct substrate of ATM signaling. Likewise, the γH2AX clearance kinetic after IR-exposure was also markedly delayed when MCF10A cells were stably BRF2-depleted using shRNAs ( Figure 2I,J), suggesting that DNA damage persists in BRF2-depleted cells, suggesting a potential role in regulation of DNA repair.

Targeting BRF2 to Interrupt DNA Binding
Given the upregulation of BRF2 in cancer and its link with DNA damage induction, we aimed to target BRF2 using available drugs. BRF2 interacts with TBP and this complex can bind to DNA[1] (Fig 3A). The structure of BRF2 consists of a Zn ribbon, N/C cyclin repeats, and an arch, which binds the C-terminal and TBP anchor domains to the rest of protein (Fig 3 B, C). Initially, we attempted to target different available pockets of this protein for inhibiting its function, particularly Cysteine 361 which is considered as an oxidative stress switch in BRF2[1]. Supplementary Table 1 lists the drugs from Drug Bank and NCI which showed the best energy binding affinity, and Supplementary Figure 4 shows the architecture of these at the C361 binding site according to the molecular docking. However, molecular dynamics simulation for all of these drugs in this pocket as well as the binding site to DNA failed to demonstrate a stable bond (data not shown).
Given the difficulty in targeting BRF2 at the C361 site, we instead aimed to inhibit the BRF2-TBP binding complex. For this, we first evaluated the binding architecture between BRF2 and TBP using molecular docking (Fig. 3D). Then, we used molecular dynamics (MD) simulation to investigate the stability of the bond. A 100-ns MD simulation was performed to investigate the dynamic stability of free BRF2 compared to its ligandbound form. The Root-Mean-Square Deviation (RMSD) can be applied to compare the initial and final structural conformations of a protein backbone. During the course of the simulation, the relative stability of the protein to its conformation is determined by the number of deviations, with fewer indicating greater stability. The RMSD value for the Cα Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 December 2020 doi:10.20944/preprints202012.0770.v1 backbone was calculated for 100 ns in order to evaluate the stability of both systems. The RMSD profile indicated that during the initial periods of simulations, the BRF2-free structure deviated considerably from the X-ray structure [1]. RMSD curves showed that the backbone trajectories of free-BRF2 structure were stable and reached equilibrium after the first 30 ns of the simulation, which ranged between 8 and 9 Å (Fig. 3E). Also, our 3Dstructural analysis showed that at the end of the MD simulation, domain 2 of the protein rotated clockwise (Fig. 3F, Video 1). This conformational change in BRF2 after binding TBP is an important feature that might play a key role in binding to DNA. Therefore, we aimed to investigate the change in conformation in more detail.

Principle Component Analysis (PCA) and Free-energy Landscape
Principle Component Analysis was carried out to understand the conformational changes of BRF2 during the MD simulation [18,19]. Figure 4A represents a plot of the eigenvalues obtained from the diagonalization of the covariance matrix of the Cα atomic fluctuations. The reduction in the eigenvalue amplitude indicates a shift from concerted motions to more constrained, localized fluctuations. This analysis suggested that the first two eigenvectors obtained from PCA can account for a higher percent of the total motions in all simulations. As seen in Figure 4A, the properties of motions for BRF2 described by the first two PCs are different.
To better understand the conformational changes of BRF2 protein at the free state, the conformational spaces of BRF2 were generated to gain significant information by using projections of MD trajectories on the first two PCs. Also, to deeply understand the motions, the displacements of eigenvector1 and 2 were calculated (Fig. 4B).
As shown in Figure 4B, the overall motion of the BRF2 protein has a different subspace of structures during 50 ns MD simulations, especially in PC1 mode. Based on Figure  4B, the protein visits three conformational clusters in PC1. RMSF analysis of PC1 showed that these clusters correspond to enhanced displacements in the regions of 1-30, 70-95 and 115-130 residues (Fig. 4C).
For studying the differences in motion behaviors between the PC1 and PC2 modes, two porcupine plots are displayed in Figure 4D which employ the first and second eigenvector using the VMD software. The direction of the arrow is an indicator of the collective motion direction, and the length of the arrow scales the strength of the movements. Based on Figure 6, the conformation of free BRF2 is different from TBP-DNA-BRF2 complex form of this protein.
In order to determine the low-energy basins (minima) and stability of the protein explored during the simulation, the FEL values of BRF2 were constructed using the projections of their own first (PC1) and second (PC2) eigenvectors (Fig. 4E). Dark blue spots indicate the energy minima and energetically favored protein conformations, and yellow spots reflect the unfavorable conformations. The shallow and narrow energy basin observed during the simulation revealed the low stability of BRF2. Based on this figure, BRF2 showed a cluster consisting of two connected energy minima basins close to each other.

Targeting BRF2-TBP binding using Drugbank Ligands
To interpose BRF2-TBP binding complex by an effective drug from the Drug Bank database, we performed a virtual screening using molecular docking with the rigid protein and flexible ligands. As mentioned above, the optimized output structure of BRF2 after MD was utilized as a target for the virtual screening step, and the top drugs were ranked based on the lowest binding-energy value. The list of drugs (Table 1) and architecture of binding in the top 10 drugs has been shown in Fig. 5. As can be seen in docking representations, the retinoic acid receptor gamma (CD564), which is under investigation, fits at the pocket of the site better than the large aromatic macrocyclic compound, Phthalocyanine. Alectinib is a second generation oral drug that selectively inhibits the activity of anaplastic lymphoma kinase (ALK) tyrosine kinase. Pranlukast is a cysteinyl leukotriene receptor-1 antagonist. It antagonizes or reduces bronchospasm caused, principally in asthmatics, by an allergic reaction to accidentally or inadvertently encountered allergens.
-9.6 Pranlukast Under investigation -9.5 DB07456 TMC647055 has been used in trials studying the treatment of Hepatitis C, Chronic Hepatitis C, Hepatitis C, Chronic, and Chronic Hepatitis C Virus. Therefore, we chose CD564 for simulation and performed MD for both BRF2 and drug (Fig. 6A, Sup Video 2). The RMSD value of the ligand was stable until 16000 ps followed by a sharp rise from 0.6 to 0.9 Å with fluctuations around this value until the end of the simulation. The analysis in Figure 6A indicates that the RMSD value of the protein reached equilibrium and oscillated around the 0.2 nm for approximately 7000 ps of simulation time. This evidence clearly indicates that the whole system was stable and equilibrated, and that CD564 is a suitable compound to inhibit the binding site of BRF2 and TBP. The final structure of the complex and ligand binding site is shown in Figure 6B-C and Video 2.

Thermodynamic Parameter Calculations
One of the most important analyses after MD is the calculation of binding free energy (ΔG). Using g_mmpbsa (Molecular Mechanic/Poisen-Boltzman Surface Area, MMPSA) software and molecular dynamic simulations, we calculated the relative binding free energy. Four hundred extracted snapshots were used to calculate the binding energy. The binding free energy and its corresponding components obtained from the MM/PBSA calculation of the ligand-BRF2 complexes are listed in Table 2. The results indicated that the ligand possessed a high negative binding free energy value of -412.059 kJ/mol. Further evaluation showed that electrostatic energy contributed most to the protein binding to the ligand, where the carboxyl group bound to Arg81 with electrostatic interactions, causing the total electrostatic energy to become negative ( Figure 5D). Van der Waals interactions also played a critical role in this binding as ∆Gsolv-polar was 197.664 KJ/Mol. The ∆GMMP-BSA was also indicative of strong, high affinity binding between protein and ligand. For calculating absolute ligand-protein complex binding energy, the ligand pulled 3 nm from the binding pocket. The binding free energy (DG) was calculated by taking the difference between the highest and lowest value of the PMF graphs. This unbinding process requires −10.29 kcal/mol energy to dissociate the ligand ( Figure XA). As shown in figure 6C , the ligand has carboxyl and aromatic ring sub-structures. Based on the figure  6D , the ligand has strong interactions with surrounding residues, especially TRP131, PHE180, GLN84 and ARG81. The pulling simulation revealed that an extended stay inside the binding pocket was the result of interactions with several amino acids. During the dissociation process, we have observed four energy minima ( Figure 6E 1a-1d). The conformation of the mentioned four energy minima were shown in figure XF. As shown in figure 6F, at the conformation of 1a the carboxyl group of ligand make hydrogen bond and electrostatic interactions with Tyr260 and Arg81, respectively. Also rings 3 and 4 established pi-pi stacking with Trp131 and Tyr176. At the second energy minimum (1b) the mentioned hydrogen and electrostatic interactions were broken, but the pi-pi stacking interactions still was remained. At the third energy minimum (1c), the carboxyl group of ligand make new electrostatic interaction with Lys181. The fourth energy minimum (1d), obtained after breaking off all mentioned contacts, and the ligand started to come out from the binding cavity. The molecule again pulled continuously, and PMF graph attained equilibration at around ~3 nm distance. At this position, the molecule was entirely unbound and became solvent-exposed. Also, figure 6G and Sup video 3 represent the unbinding pathway of the ligand from the binding site.

Discussion
It has been proposed that targeting transcription factors may offer a therapeutic option with lower toxicity and increased efficacy than current chemotherapy drugs, which opens opportunities for novel drug-development approaches [20]. Classical chemotherapy agents targeting DNA integrity and cell division have the potential to give rise to secondary malignancies. Current strategies target oncogenic molecules or signaling pathways that tumor cells are dependent on for growth and survival. This phenomenon is known as oncogene-addiction. However, there is emerging interest in targeting multi-component cellular machineries that are not directly involved in growth and proliferation to maintain the malignant properties that contribute to tumorigenicity (non-oncogenic addiction) [21,22]. BRF2 and the other class III genes transcribed by Pol III are known as 'housekeeping genes' and their expression is required in all cell types for cellular viability. Housekeeping genes have been suggested as anti-cancer targets with therapeutic potential given the likelihood of reduced associated toxicity [23]. TFs are targeted by approximately 10% of all FDA-approved anti-cancer drugs [24]. However, most TFs are not targetable because of their localization, or because their components and subunits do not possess enzymatic activity, which is usually easier to target with specific drugs. Therefore, we endeavoured to identify a mechanism by which BRF2 could be targeted. However, most of the accessible pockets and available domains in BRF2 were not targetable, and we could only target BRF2 and TBP binding sites to impair their DNA-binding capacity. The formation of a complex between various proteins has a key role in many biological processes, and consequently, preventing complex formation interrupts the associated processes. However, BRF2 plays a key role in protecting cancer cells against oxidative stress and our data additionally suggest that BRF2 plays a role in DDR regulation. Therefore, inhibition of this target can be more valuable than targeting other TFs or housekeeping genes.
Recently, there has been increased interest in targeting polymerases, such as the targeted destruction of the Pol I, to treat patients with advanced solid tumors, which ultimately expands therapeutic opportunities [25]. Notably, targeting components of Pol I using the selective inhibitor CX5461 to treat tumors with MYC amplification has shown selective killing of tumor cells through activation of non-canonical ATM signaling [26]. Additionally, different mode of action for CX5461 as TOPO isomerase II poison or G-quadruplex-mediated DNA-damage have been proposed [27,28]. Regardless of the difficulties in targeting RNA polymerases and their associated regulators, the recent discoveries of drugs that are able to target these proteins underline their potential in the future of cancer therapy. Oncogenes, such as MYC, or tumor-suppressors, such as RB and p53, are multiple regulators of Pol III and Pol II, which make them difficult to distinguish. Furthermore, these common regulators have a global effect on Pol III, which makes it more complex to assign a specific gene contribution to Pol III deregulation. Due to its pivotal role in cell growth, proliferation and tumorigenesis, Pol III activity is tightly regulated by a series of activators and repressors. For example, the Pol III repressor MAF1 is negatively regulated by a series of phosphorylation events on seven different serine residues [29]. Similarly, RB and the RB-related family members repress Pol III activity and heightened levels of PolIII transcripts are implicated in accelerated growth associated with RB dysfunctional tumors [30,31].
The Casein Kinase 2 (CK2) is an extremely conserved kinase that plays pivotal roles in many aspects of cell physiology; moreover, CK2 is often overexpressed in several human cancers [32]. Among its many substrates, CK2 also regulates several transcription factors and polymerases, including Pol III [33] and the TFIIIB sub-factor TBP [34]. Another kinase involved in Pol III phosphorylation is ERK. In response to mitogen activation, the mitogen-activated protein (MAP) kinase ERK regulates Pol III activity through direct binding and phosphorylation of BRF1 [35]. The phosphatase PTEN plays a fundamental role in tumor suppression by inhibiting the PI3K signalling, and it is often mutated or loss in human cancers [36]. It has been shown that PTEN represses the activity of Pol III by blocking the formation of functional TFIIIB complexes. Specifically, PTEN blocks the association between BRF1 and TBP: this inhibition occurs through the reduction of phosphorylation of BRF1 that induces the dissociation of TFIIIB [37]. We performed a comprehensive bioinformatic analysis, correlation studies and literature review on the effect of BRF2 on signalling pathways. However, our analysis reveals unlike POLIII transcript itself, BRF2 transcript is not linked to the above-mentioned signalling pathways. Notably, BRF2 mRNA expression in human tumors also did not correlate with any feature of genomic instability as an indicator of DNA damage (Suppl fig 3) which is supported by our in vitro data in human breast cancer lines (Fig 2) that nicely show that BRF2 protein rather than transcript level is regulated after exposure of cells to DNA damaging agents..Thus, it is likely that BRF2 might be post-transcriptionally regulated by above-mentioned signalling pathways.
We further characterized the relationship between BRF2 and DDR experimentally. BRF2 is a major oxidative stress regulator, and oxidative stress mechanisms can drive tumorigenesis [2]. A BRF2-regulated gene product, known as SeCys tRNA, is involved in the synthesis of selenoproteins. The main function of SeCys tRNA is to reduce reactive oxygen species (ROS) and maintain the cellular redox homeostasis through selenoprotein synthesis[1]. The absence or defective expression of selenoproteins induces apoptotic cell death. Studies have shown a high basal level of BRF2 expression in breast [2] and lung cancer [3], as well as in cancer cell lines such as HCC95 and which are subject to prolonged stress, such as oxidative stress [4]. Therefore, we speculated that the high levels of BRF2 expression in cancer cell lines assists in maintaining sufficient expression of selenoproteins in order to detoxify ROS and preserve redox homeostasis. DNA damage induction is linked to oxidative stress, which can in turn cause further DNA lesions. It was conspicuous that in the investigated cancer cell lines and near-normal breast cell line (MCF10A), the expression of BRF2 was upregulated after IR. The cisplatin experiment represents DNA damageinducing treatment to exclude BRF2 as possibly being activated by oxidative stress, which is triggered by radiation [38,39]. BRF2 demonstrates clear upregulation following cisplatin treatment, which implies a possible link between BRF2 and DNA damage. Additionally, the expression of NRF2 (a protein stimulated by oxidative stress [40]), was examined as well after cisplatin, and no significant changes were detected. This strengthened the possibility of the link between BRF2 and DNA damage, without the confounding effect of oxidative stress. To investigate the functional role of BRF2 in the DNA damage response, the downstream substrates, as well as the DNA damage sensor, were examined. We found that blocking ATR but not DNA-PK or ATM has direct role in BRF2 upregulation.
After characterizing the potential role of BRF2 in DDR and OS signaling, we initially aimed to target C361 in BRF2 which is a critical OS on/off switch. However, according to molecular simulations, we failed to introduce any effective drugs. Therefore, we concentrated on targeting the BRF2-TBP-DNA complex. The high protein RMSD value along the trajectory indicates that a significant change in the molecule, compared to the original structure, had occurred. However, the stability of the protein after 30 ns, persisting at a value of approximately 1 nm, illustrated its stability after the large conformational displacement. PCA is the most common method to reveal the critical motions of a protein which we used as confirmation for conformational changes. The PC analysis was performed to reveal the origin of the large conformational difference between the initial and final MD structures, which is supported by their different motion direction relative to each other. Analyses of residue displacements based on the projection of MD trajectories on the first and second component proved that larger displacements of some residues, located at the N-terminus of protein, make two conformations visiting different subspaces during MD simulations. This provided significant dynamics information involving conformational diversity of BRF2 protein in both the bound and free states. Free Energy Landscape analysis provided intimate details about the stable conformation and energy minima. BRF2 FEL analysis showed two energy minima that represent the two conformations of Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 December 2020 doi:10.20944/preprints202012.0770.v1 protein (bound and unbound). Finally, the stable conformation of this step was used as a receptor to perform virtual screening. A structure-based virtual screening method was used to select the ligand with a high binding energy value. Our analysis showed that the ligand has an acceptable binding energy value. The MD simulation method was used to study the stability of ligand at the binding site. RMSD plot confirmed that the protein and ligand are stable during the MD simulation. At the next step, the relative binding free energy of ligand to protein was calculated by the MMPBSA method. MMPBSA analysis showed that the van der Waals and electrostatic interactions and non-polar solvation energy negatively contribute to the total interaction energy while only polar solvation energy positively contributes to the total free binding energy. In terms of negative contribution, electrostatic interaction contributed more than Van der Waals interactions. Also, the interaction plot generated by ligplot showed that the carboxyl group at the ligand makes electrostatic interaction with Arg81 of protein. Collectively, using in silico approaches we have identified putative inhibitors of BRF2-TBP-DNA complex after virtual screening of more than 3500 FDA approved or under clinical trial investigational drugs. We confirmed the drug binding to BRF2 by using molecular dynamics simulation and further experimental validation is required to confirm in silico analysis.

Conclusions
In conclusion, we showed BRF2 as a potential therapeutic target in cancer by comprehensively analyzing publically available Pan-Cancer data and we also discovered regulation of BRF2 after DNA damage in breast cancer which its overexpression appears to be mutually exclusive to BRCA1 and BRCA2 loss in breast cancer. Notably, we using virtual screening on more than 3500 available drugs we found inhibitors for BRF2 -TBP-DNA complex binding which successfully validated by molecular docking, energy binding and MD simulation results. MD results revealed that BRF2 shows a conformational change after binding to TBP, which leads to changes in the available pockets. Using this structure, we performed virtual screening and identified CD564 as a promising candidate that can create a stable binding to the binding site of BRF2-TBP and consequently, interposes DNA binding with the functionality of BRF2. We validated these results by using MD simulation where RMSD charts indicated the stability of bound and thermodynamic parameters verified by the binding energy. The results require experimental validation to be pursued for therapeutic and clinical applications; however, this research identified a promising cancer target, implying new roles in DDR regulation and highlighting the potential for rational drug repurposing using in silico methods.
(from The Cancer genome Atlas (TCGA) RNAseq datasets) and 18 measures of genomic instability (GI). The 18 measures are labelled "a" through "r". For all panels, the y-axis shows the relative expression level of BRF2. Details on each of the GI variables on the x-axis of each panel is available in methods. In all scatter plots, linear regression analysis was performed and R and p values are given according to Spearman's rank correlation. Samples sizes, n, are also included. Panel (g) intentionally does not include linear regression analysis, because there is no scientific value for that for this data type. Instead, these data were binned and included in panel (h). Violin plots in panels (h), (m) and (r) show p values calculated according to Mann-Whitney U tests, Figure S3: A-D: The heatmaps show the correlations between gene expression levels and the variables listed on the y-axis in up to 32 cancer types, whose abbreviations are listed on the top x-axis. For more details see the methods. Each tile in the heatmaps shows the Spearman correlation between the gene expression levels (mRNA levels from The Cancer Genome Atlas (TCGA) RNAseq datasets)) and the variables listed below. The colour of each tile reflects the Spearman correlation coefficient (r), as indicated by the colour key on the right. Each tile shows the Spearman p value abbreviations in the following format: blank tiles: p > 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001, Figure S4: Immune cell infiltration: The levels of indicated tumour-infiltrated immune cell types were estimated using the tumor immune estimation resource, for more details see the methods and for statistics methodology see Sup 3, Figure S5: A-L: Docking shows architecture of binding of indicated drugs to C361 amino acid of BRF2 (black arrowed). Figure S6: A-E. Molecular docking and interaction between amino acids of BRF2 with ligands as indicated in figure (compounds listed 10-15 in ranking), Table S1: Binding energies of top-ranking docked compounds from Drug Bank against C361 amino acid of BRF2), Table S2: Table 2. Binding energies of top-ranking docked compounds from NCI against C361 amino acid of BRF2.