Identification of Effective Anticancer G-Quadruplex-Targeting Chemotypes through the Exploration of a High Diversity Library of Natural Compounds

In the quest for selective G-quadruplex (G4)-targeting chemotypes, natural compounds have been thus far poorly explored, though representing appealing candidates due to the high structural diversity of their scaffolds. In this regard, a unique high diversity in-house library composed of ca. one thousand individual natural products was investigated. The combination of molecular docking-based virtual screening and the G4-CPG experimental screening assay proved to be useful to quickly and effectively identify—out of many natural compounds—five hit binders of telomeric and oncogenic G4s, i.e., Bulbocapnine, Chelidonine, Ibogaine, Rotenone and Vomicine. Biophysical studies unambiguously demonstrated the selective interaction of these compounds with G4s compared to duplex DNA. The rationale behind the G4 selective recognition was suggested by molecular dynamics simulations. Indeed, the selected ligands proved to specifically interact with G4 structures due to peculiar interaction patterns, while they were unable to firmly bind to a DNA duplex. From biological assays, Chelidonine and Rotenone emerged as the most active compounds of the series against cancer cells, also showing good selectivity over normal cells. Notably, the anticancer activity correlated well with the ability of the two compounds to target telomeric G4s.


Introduction
Advanced technologies have allowed an in-depth understanding of the structural and biological features of G-quadruplex (G4) structures and evidence pointed to the biological relevance of G4 nucleic acids, particularly as targets in anticancer strategies [1][2][3][4]. To date, several synthetic compounds have been identified as selective ligands able to bind and stabilize G4 structures, with some of them showing effective anticancer activity in vivo and, therefore, being evaluated in advanced clinical trials [5][6][7][8][9]. Nevertheless, none of the most promising G4 ligands have been approved as a drug thus far. For the real progression of anticancer therapies based on G4-targeting ligands the investigation of large libraries of small molecules endowed with high chemical diversity is, therefore, strongly needed to identify novel chemotypes. In this context, natural compounds have been poorly studied as G4 ligands compared to synthetic compounds [5,[10][11][12], though representing, in principle, appealing candidates due to the remarkable structural diversity of their scaffolds.
Plants are a rich source of structurally diverse secondary metabolites, which can be exploited in the development of new drug candidates [13,14]. Due to their high biodiversity, medicinal plants provide a huge number of natural compounds [15][16][17]. However, less than 1% of this biodiversity has been exploited in drug discovery due to several factors including the lack of a proper multidisciplinary view. The advent of powerful and efficient methods, such as the integrated combination of combinatorial chemistry and High-Throughput Screening (HTS), as well as user-friendly informatics tools, such as computer-aided drug design, which met the demand of major pharmaceutical companies to accelerate the research process, promoted the revolution of natural products screenings in drug discovery [18,19]. In addition, the development of new techniques for the isolation and characterization of novel compounds significantly improved the efficiency of the processes, in which the major challenges can currently be identified in the generation of high-quality libraries of diverse natural products that might allow the fast identification of lead compounds of pharmacological relevance [20][21][22].
A unique high diversity library composed of ca. one thousand individual natural products, isolated mainly from indigenous plants collected in biodiversity-rich countries, especially in tropical and subtropical areas, and enlarged with their semi-synthetic and synthetic derivatives, is available from the Organic Chemistry Laboratory of the Department of Chemistry and Technology of Drugs of Sapienza University of Rome, Italy [23]. During the years, the exploitation of this in-house collection of natural products offered a unique chance to identify unexpected new scaffolds for the development of therapeutically relevant molecules. Furthermore, the successful application of computer-aided methods in screening this unique and diverse in-house library provided some lead compounds that have been developed and, in some cases, patented as anticancer and antimicrobial agents [24][25][26][27][28]. Here, a docking-based virtual screening has been carried out, using both telomeric and oncogenic G4 models as targets [29,30], to evaluate the ability of the in-house natural products to target G4 grooves and identify novel G4-targeting chemotypes. Groove and loop binders are expected to be more selective than compounds that stack on top of the guanine quartets, although structural details of the highly flexible G4 loops are generally not precise enough to allow a reliable docking-based virtual screening in these regions compared to G-tetrads [31]. The virtual screening process identified 28 potential selective G4 ligands (Table S1) which, to the best of our knowledge, have not been previously investigated in their interaction with G4 structures, with the only exception of Aloin and Chelidonine for which only preliminary studies are reported in the literature [32,33]. Then, the ligands prioritized in silico have been experimentally screened by exploiting the G4-CPG (G-quadruplex on Controlled Pore Glass) assay, an affinity chromatography-based method to efficiently and quickly identify G4 selective ligands [11,[34][35][36][37][38][39]. The compounds showing the highest affinity and selectivity for the selected G4 models have been studied for their interaction with the G4-forming sequences of choice in solution, in parallel with a duplex structure as control, by using circular dichroism (CD) and fluorescence spectroscopies. Additionally, molecular dynamics (MD) simulations and biological studies have been carried out for the best G4 ligands in order to get a deeper insight into their binding behavior towards G4 structures as well as their antiproliferative activity on cancer and normal cells.

Docking-Based Virtual Screening of the In-House Library of Natural Compounds
The unique high diversity in-house library investigated here consists of fully characterized natural products and their derivatives belonging to different classes of organic compounds, including variably substituted flavonoids, benzophenones, xanthones, anthraquinones, alkaloids, steroids, terpenoids etc. It was then enlarged with natural compounds from commercially available sources and semi-synthetic and synthetic compounds [23]. A docking-based virtual screening was carried out with the AutoDock program, using the solution structure of both a telomeric and an oncogenic G4 target as rigid receptors [29,30], and the solution structure of a DNA duplex as a negative control for the identification of potentially G4 selective virtual hits. As previously described [36], the binding site was centred on the G4s groove in the search for ligands that might be more selective than compounds that stack on top of the G4 tetrads. In detail, the rectangular binding site was centred in the groove formed by G4-G6, T8 and G22-T24 in tel 26 and the groove formed by G2, A3, G5, G6 and G17-G19 in c-myc covering the groove and the loops. Although we are aware that this selection might not consider potential binders to other sites such as G-tetrad stackers, we believe that groove binders might be more profitable than unspecific stackers for further development. In contrast, the entire surface of the unspecific DNA duplex was scanned in the docking-based virtual screening. After docking, compounds were ranked based on their theoretical affinity within the groove of the G4 targets. The ligand-binding mode was visually inspected, and the compounds were further filtered by chemical diversity through a custom Python script for compounds clustering based on a combination of fingerprints and substructure search [40]. This latter step was aimed at reducing chemical redundancy within the test set and exploring the largest portion as possible of the natural products space represented in the in-house library. Specifically, the compound with the best score of each chemical cluster was selected for further processing. Finally, a comparison of the ligands scores in binding to the target G4s and the unspecific duplex sequence led to the final selection of 28 candidate hits (Table S1). A considerable part of these compounds belongs to the chemical class of alkaloids (1-13), one of the largest and most intriguing families of natural compounds. Indeed, alkaloids are characterized by vast structural diversity with no uniform classification. Several candidate hits belong to the chemical class of polyphenols (anthranoids, flavonoids and benzophenones). Within the anthranoids subclass, which can be chemically described as dihydroxy-anthraquinones, -dianthrones and -anthrones, three ferruginines (14, 15 and 16), two anthrones (17 and 18) and one vismione (19) were selected. Among the flavonoids subclass, two rotenoids (20 and 21), containing a cis-fused tetrahydrochromeno [3,4-b]chromene nucleus, were identified. Among the polyphenols included in the library, a benzophenone (22) was identified. Three compounds were steroids (23, 24 and 25), a subclass of terpenes featured by a characteristic molecular structure composed of 17 carbon atoms arranged in four rings. A Diels-Alder type adduct (26) and a naturally occurring dibenzofuran derivative (27), together with a cyanogenic glycoside (28), complete the test set.

Experimental Screening of the Library of 28 Natural Compounds by the G4-CPG Assay
The 28 natural compounds in silico selected as G4 ligands were experimentally evaluated for their ability to interact with G4 structures by the G4-CPG assay, an affinity chromatography-based method for the screening of putatively selective G4 ligands [11,[34][35][36][37].
As in molecular docking, two cancer-related G4-forming DNA sequences originated from human telomeres (tel 26 ) or c-myc oncogene promoter (c-myc) were used as the targets [29,41]. In parallel, the interaction with a model unimolecular duplex-forming DNA sequence (ds 27 ) was also examined to determine if the analyzed compounds could discriminate G4 vs. duplex DNA structures [37].
The 28 compounds were dissolved in DMSO to prepare stock solutions. Then, they were evaluated for their solubility in the concentration and in the washing/releasing solutions used in the G4-CPG assay (see Materials and Methods for details). All tested compounds proved to be fully soluble and stable in the assay experimental conditions, except for Amygdalin, Clusiacitran B, Digitonin, Diosgenin, Solanidine and Usnic acid.
Successively, for each of the remaining 22 natural compounds, tests were carried out to evaluate, firstly the possible unspecific binding on the nude CPG support, and then their ability to bind G4-and duplex-forming DNA. The results of the G4-CPG assay are summarized in Table 1. Generally, no significant unspecific binding was observed, thus, not precluding further analysis on the oligonucleotide-functionalized supports. However, some unspecific binding on the nude CPG support was detected for Ferruanthrone, Ferruginin A, Ferruginin B, Hydrastine, Kuwanon G, Narceine and Rotenone. This was taken into due consideration in the following selection of the best G4 ligands. Table 1. Summary of the binding data obtained for the 28 natural compounds by the G4-CPG assay and selectivity indexes calculated as the ratio between the percentages of ligand bound to the indicated G4-and duplex-functionalized supports. Bound ligand was calculated as a difference from the unbound ligand, recovered with 50 mM KCl, 10% DMSO, 10% CH 3 CH 2 OH washing solution, and is expressed as % of the amount initially loaded on each support. The errors associated with the % are within ±2%.  Overall, 20-OH-Ecdysone, Aloin, Aspidospermine, Ferruginin B, γ,γ'-OH-Ferruginin A, Rotenolone, Veratrine, Vindoline, Vismione B and Yohimbine proved to weakly interact with G4-and duplex-functionalized supports. On the other hand, Emetine, Ferruanthrone, Ferruginin A, Hydrastine, Jervine and Narceine showed a good affinity for both G4-and duplex-functionalized supports, with no significant difference in terms of the percentage of bound ligand to the different secondary structures of the DNA investigated.
Conversely, Bulbocapnine, Chelidonine, Ibogaine, Kuwanon G, Rotenone and Vomicine showed a good affinity towards G4-functionalized supports and low-to-null inter-actions with the duplex-functionalized support, as also evidenced by the related indexes of G4/duplex selectivity (Table 1). However, Kuwanon G was not considered for further investigations due to its affinity towards nude CPG similar to tel 26 -functionalized CPG.
In summary, Bulbocapnine, Chelidonine, Ibogaine, Rotenone and Vomicine ( Figure 1) were selected as the best G4 ligands in terms of affinity and selectivity over the DNA duplex, according to the results of the G4-CPG assay and were further analyzed in solution by biophysical techniques in order to achieve additional information on their interaction with G4 structures.

Circular Dichroism Studies
Considering the results of the G4-CPG assay, the ability of Bulbocapnine, Chelidonine, Ibogaine, Rotenone and Vomicine to interact with telomeric and oncogenic G4 models as well as a control duplex, was investigated in solution by CD experiments. All DNA oligonucleotides were prepared by overnight annealing tel 26 , c-myc and ds 27 solutions at 2 µM DNA concentration, in 5 mM KCl, 5 mM phosphate buffer, 5% DMSO (pH 7). In full agreement with the literature data, in the above conditions, we found that: (i) tel 26 folded into a hybrid G4, featured by a double hump-band, with maxima centered at 290 and 265 nm [42], (ii) c-myc adopted a parallel G4 topology, with a maximum centered at 262 nm and a minimum at 242 nm [42] and (iii) ds 27 showed a positive band at 280 nm along with an intense minimum at 251 nm, characteristic of a B-DNA duplex structure (see Figures S1-S3, black lines) [43].
Then, the three secondary structure-forming oligonucleotides were titrated with increasing amounts of the five selected compounds (up to 10 molar equivalents), and the corresponding CD spectra were recorded after each addition ( Figures S1-S3).
Bulbocapnine, Chelidonine, Ibogaine, Rotenone and Vomicine present one or more chiral centers. Aiming at evaluating the contribution of these ligands to the CD spectra obtained from titration experiments with the oligonucleotides, control CD spectra of these five compounds were recorded by adding increasing amounts of each ligand to the buffer alone, thus reproducing the above titration experiments but in the absence of DNA ( Figure S4).
Thus, the contribution of each ligand was subtracted from the CD spectra obtained upon titrations of the tel 26 and c-myc G4s or ds 27 duplex, obtaining a more accurate picture of the spectral changes of the DNA oligonucleotides induced by each ligand (Figures S5-S7).
After these subtractions, it emerged that for the tel 26 /Ibogaine, tel 26 /Rotenone and tel 26 /Vomicine systems only slight spectral changes of tel 26 G4 were detected ( Figure S5C-E). On the other hand, relevant spectral changes were observed for tel 26 G4 upon addition of Bulbocapnine or Chelidonine ( Figure S5A,B).
For all the c-myc/ligand systems, a slight dose-dependent decrease of the CD intensity of the 262 nm band was observed ( Figure S6). In addition, a slight dose-dependent increase of the CD intensity of the 242 nm band was found for the c-myc/Bulbocapnine and c-myc/Rotenone systems ( Figure S6A,D).
As far as titrations of ds 27 duplex are concerned, no relevant spectral changes were detected by inspection of the CD spectra after ligand contribution subtraction ( Figure S7). Only slight spectral changes were found for the ds 27 /Bulbocapnine and ds 27 /Rotenone systems, with a dose-dependent decrease in CD signal intensity of the 280 nm band ( Figure S7A,D).
To semi-quantitatively evaluate the spectral changes detected by titration of the tel 26 and c-myc G4s or ds 27 duplex, the differences (∆CD) between the CD intensity of DNA/ligand 1:10 ratio systems (upon ligand contribution subtraction) and the CD intensity of DNA alone were calculated considering the CD values at 290, 262 and 251 nm for tel 26 , c-myc and ds 27 systems, respectively. The obtained ∆CD values-to be intended as easy-to-handle parameters indicative of a trend, and not as quantitative data for the description of an observed effect-are reported in Figure 2. In detail, the most significant effects on tel 26 and c-myc G4 structures were induced by Bulbocapnine and Chelidonine, while the most relevant spectral changes on the ds 27 duplex were detected for Bulbocapnine and Rotenone. Significantly, the highest ∆CD values for Bulbocapnine, Chelidonine, Ibogaine and Vomicine were found for the two investigated G4s compared to the control duplex structure, confirming their good G4 vs. duplex selectivity evidenced by the G4-CPG assay. On the other hand, Rotenone seemed to more markedly affect the duplex than the G4 structures. This finding is in line with its lower G4 vs. duplex selectivity compared to the other four compounds, as observed by the G4-CPG assay.
CD-melting experiments were also performed on all the DNA/ligand mixtures in the 5 mM KCl, 5 mM phosphate buffer, 5% DMSO (pH 7) to evaluate if stabilizing or destabilizing effects on the G4 and duplex structures were obtained upon incubation with each of the five ligands. CD-melting curves of tel 26 and c-myc G4s or ds 27 duplex in the absence or presence of each ligand (DNA/ligand 1:10 ratio) were recorded by following the CD changes at the wavelength of maximum intensity (290, 262 and 251 nm for tel 26 , c-myc and ds 27 , respectively) ( Figure S8A-C). Melting temperatures for all the investigated systems are summarized in Table 2. Table 2. Melting temperatures (T m ) of tel 26 and c-myc G4s or ds 27 duplex in the absence or presence of the ligands (10 molar equivalents) as measured by CD-melting experiments. Moderate to low stabilizing effects on tel 26 G4 were found for Chelidonine (∆T m = +4 • C), Ibogaine (∆T m = +2 • C), while no effect was detected for Bulbocapnine, Rotenone and Vomicine ( Figure S8A and Table 2).
As concerns the c-myc systems, stabilizing effects were observed for all the investigated ligands. However, the T m values could not be accurately determined for these systems due to the absence of good sigmoidal behavior of the related melting curves in the 10-95 • C range ( Figure S8B). In order to overcome this drawback, melting experiments for the c-myc systems were also performed in a buffer containing a ten-fold lower potassium ions concentration, i.e., 0.5 mM KCl, 0.5 mM phosphate buffer, 5% DMSO (pH 7) ( Figure S8D). Even at this very low K + ion concentration, c-myc proved to fold in a G4, featured by parallel topology and a T m of 45 • C ( Figure S9). Under these conditions, melting temperatures were obtained for the c-myc/ligand systems showing a strong stabilizing ability of these compounds on c-myc G4 with ∆T m = +13, +9, +5, +5 and +4 for Chelidonine, Bulbocapnine, Ibogaine, Vomicine and Rotenone, respectively ( Figure S8D and Table 2).Notably, no or even destabilizing effects on ds 27 duplex were observed for all the five compounds ( Figure S8C and Table 2).
Overall, in full agreement with the G4-CPG assay and CD titrations, the results of the CD-melting experiments demonstrated that all the five investigated compounds showed an excellent ability to interact with G4s, also discriminating well the G4 vs. duplex structures.

Fluorescence Spectroscopy Studies
In order to get a deeper insight into the affinity of the selected natural compounds to the selected G4 and duplex DNA, fluorescence experiments were performed. First, the fluorescence behavior of the selected compounds was evaluated. In detail, fluorescence spectra of 2 µM solutions of each ligand were recorded in 5 mM KCl, 5 mM phosphate buffer, 5% DMSO (pH 7). Bulbocapnine, Chelidonine and Ibogaine showed strong emission bands at 460, 330 and 356 nm, respectively (Figures S10-S12, left panels, black lines). Conversely, Rotenone and Vomicine did not show appreciable fluorescence intensity, thus hampering further investigations on DNA binding by fluorescence spectroscopy.
Then, fluorescence titrations were carried out for Bulbocapnine, Chelidonine and Ibogaine (Figures S10-S12, left panels) at a fixed concentration of ligand (i.e., 2 µM), by adding increasing amounts of tel 26 and c-myc G4s or ds 27 duplex, previously annealed in 5 mM KCl, 5 mM phosphate buffer, 5% DMSO (pH 7). Upon each addition, the corresponding fluorescence spectrum was recorded after the stabilization of the signal. Successively, for each system, the fraction of bound ligand was calculated from the obtained fluorescence intensity values and plotted as a function of the DNA concentration (Figures S10-S12, right panels). These data were then fitted with an independent and equivalent-sites model [44] to calculate the binding constants (K b ) ( Table 3). Table 3. Binding constants (K b ) obtained by fitting of fluorescence data for DNA/Bulbocapnine and DNA/Ibogaine systems by using an independent and equivalent-sites model [44]. 26 c-myc ds 27 Bulbocapnine Titrations of Bulbocapnine with increasing amounts of tel 26 and c-myc G4s or ds 27 duplex showed an overall fluorescence quenching (Figure S10A-C). K b values of 1.0, 1.1 and 1.2 × 10 6 M −1 were obtained for Bulbocapnine binding to tel 26 , c-myc and ds 27 systems, respectively (Figure S10D-F and Table 3).
As far as Chelidonine is concerned, an overall quenching effect was observed on increasing tel 26 G4 concentration ( Figure S11A). On the other hand, the c-myc/Chelidonine system ( Figure S11B) exhibited an initial fluorescence enhancement (from 0 to 1 µM DNA concentration) followed by fluorescence quenching (from 1 to 10 µM DNA concentration). Notably, the fluorescence for the ds 27 /Chelidonine system did not significantly change by increasing DNA concentration ( Figure S11C). Unfortunately, by plotting fluorescence intensity at the wavelength of Chelidonine intensity maximum vs. added DNA increasing concentration (Figure S11D-F), the experimental data did not follow a well-defined behavior, thus not allowing adopting fitting protocols to obtain binding constants.
As in the case of Bulbocapnine, a significant fluorescence quenching was observed upon titration of Ibogaine with all the investigated oligonucleotides (Figure S12A-C). K b values of 4.1, 3.9 and 3.7 × 10 5 M −1 were obtained for Ibogaine interaction with tel 26 , c-myc and ds 27, respectively (Figure S12D-F and Table 3).
Overall, while higher binding constants were found for Bulbocapnine with tel 26 and c-myc G4s as well as ds 27 duplex compared to Ibogaine, the latter one showed slightly higher G4 vs. duplex selectivity than Bulbocapnine. Notably, the presence or absence of significant fluorescence intensity variations for Chelidonine upon titration with the G4s or the duplex, respectively, further evidenced a good G4 vs. duplex selectivity for this compound, in agreement with G4-CPG assay and CD results.

Molecular Dynamics Simulations
The tandem application of in silico and experimental screenings proved very useful in identifying five hit binders of tel 26 and c-myc G4 structures out of many natural compounds. However, the molecular docking-based virtual screening relied on the use of static receptor structures, solved in non-physiological conditions, which can fail to address loops flexibility of the G4 targets. Thus, to investigate the coherence and persistency of binding modes identified by molecular docking, MD simulations were run on all the examined G4/ligand complexes. In addition, MD simulations were run on duplex/ligand complexes to explore structural determinants for the observed selectivity for G4s. For each complex, MD trajectories were produced for 500 ns without positional restraints, and the ligands' theoretical affinity was estimated by the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) approach [45]. Finally, MD frames were clustered, and the most representative structure was used for further discussion and structural speculations.
In binding to tel 26 G4, compounds showed a peculiar behavior (Figure 3). Chelidonine and Ibogaine preserved the binding site identified by molecular docking within the G4 groove and established H-bond interactions with the phosphate oxygen of G residues from G-quartets ( Figure 3B,C). Bulbocapnine moved towards the 3 -end of the DNA strand and established an H-bond with the OH oxygen of 3 -T, although a strong stabilization seems to be determined by stacking interactions with a T from the edgewise (or lateral) loop ( Figure 3A). Similarly, Rotenone and Vomicine detached from the groove to find their preferred binding location in a region on top of the G-quartets characterized by sequence specificity. Specifically, thanks to its L-shape, Rotenone was H-bonded to the OH oxygen of 3 -T as well as stacked on the same nucleobase ( Figure 3D), while Vomicine bound within the 5 -end of the sequence by H-bonding the phosphate oxygen of T at position 2 and being stacked by 5 -T and an A from the propeller loop ( Figure 3E).
In contrast to tel 26 systems, MD simulations on the complexes between selected ligands and c-myc provided conformational results that were more coherent with docking outcomes. Indeed, Bulbocapnine, Chelidonine and Ibogaine stably bound within the G4 grooves ( Figure 4A-C). Bulbocapnine and Chelidonine established H-bond interactions with the phosphate backbone. In contrast, Rotenone detached from the groove and moved to the propeller loop where it was sandwiched between the A and T of the loop, highlighting the propensity of this compound to establish aromatic interactions with single-stranded oligonucleotide tracts ( Figure 4D). Vomicine moved slightly towards the 3 -end of the c-myc sequence where it bound to the terminal part of the groove being stacked onto a G from the G-tetrad as well as H-bonded to the phosphate backbone ( Figure 4E).
MD simulations were also run on the docking complexes between selected ligands and the DNA duplex model, to provide a rational structural explanation of the G4 vs. duplex selectivity observed by the G4-CPG assay, as well as by CD and fluorescence studies. Notably, only Bulbocapnine was found to be firmly bound to the duplex groove in correspondence of an AT-rich sequence, establishing an H-bond with the phosphate backbone ( Figure 5A). All other compounds moved to the terminal ends of the duplex, where they bound non-specifically to the terminal base pair ( Figure 5B-E). Since, in living cells, the DNA duplex does not have frequent chain breaks, this behavior might be interpreted as a weak affinity of the ligands for this nucleic acid structure. Only Rotenone was able to perform a rather specific interaction with the terminal end of the duplex, consisting of a π-stacking with 3 -G of a strand and an H-bond with G at position 2 of the opposite strand ( Figure 5D).
Theoretical affinity of ligands to nucleic acid sequences tel 26 and c-myc was estimated by the MM-GBSA approach along with the most populated cluster of MD frames. Theoretical affinity to the DNA duplex was not reported, as the compounds were bound in a pose that is not consistent with their possible behavior in a physiological context. Results were reported in Table 4 and remarkably highlight an affinity scenario that is highly comparable to that observed by CD and fluorescence spectroscopy. Overall, theoretical affinity data suggest a stronger binding to c-myc compared to tel 26 G4, as highlighted by CD and thermal melting experiments, as well as a tighter affinity for Bulbocapnine and Chelidonine followed by Ibogaine and Vomicine, whereas Rotenone is the weakest binder of the series in agreement with experimental results.
Taken together, results of MD simulations suggest that selected ligands bind specifically to tel 26 and c-myc G4 structures with a peculiar interaction pattern, being unable to bind the groove of a DNA duplex model, except for Bulbocapnine. The agreement between computational and experimental trends corroborates the predicted binding mode and sheds further light on the mechanism of action of these natural modulators of tel 26 and c-myc G4s.

Evaluation of Biological Activity of the Identified G4 Ligands
Starting from the results of the biophysical and in silico studies, we proceeded to explore the biological activity of the identified G4 ligands (Bulbocapnine, Chelidonine, Ibogaine, Rotenone and Vomicine). To this aim, the antitumor potential of the candidate molecules was firstly assessed. In detail, BJ-EHLT cells, a line of human transformed fibroblasts, were treated with growing doses (from 0.1 to 10 µM for 72 h) of each of the five compounds and cell viability was evaluated by crystal violet assay. Notably, while Bulbocapnine, Ibogaine and Vomicine were almost ineffective, Chelidonine and Rotenone were found to produce a dose-dependent effect on cell viability ( Figure 6A), with an estimated IC50 of 0.64 µM and 0.15 µM, respectively. These results prompted a specific focus on the biological characterization of the latter two compounds. Thus, the selectivity of these biologically active compounds against malignant cells was then tested. To address this point, BJ-EHLT cells and their normal counterpart, BJ-hTERT, were exposed to treatment with Chelidonine and Rotenone for 72 h at the indicated doses, and the effect on cell viability was evaluated ( Figure 6B). Interestingly, the growth curves, besides confirming the efficacy of the two treatments on transformed fibroblasts, showed very poor activity of both the ligands on normal cells, providing evidence of the selectivity of these molecules against tumor cells. Moreover, considering that anti-tumor activity of G4 ligands mainly depends on their capability to induce DNA damage [46], BJ-EHLT and BJ-hTERT were treated for 24 h with 1 µM of either Chelidonine or Rotenone and the amount of phosphorylated histone H2AX (γH2AX), a typical hallmark of DNA double-strand breaks [47], was estimated by immunofluorescence (IF) microscopy. Quantitative analysis of fluorescence intensity of γH2AX signal, evaluated on at least 100 nuclei ( Figure 6C), showed that both the ligands are more effective in transformed than in normal cells ( Figure 6D). Notably, for both treatments, the mean of the γH2AX signal in BJ-EHLT was three times higher than in BJ-hTERT, with an increase in the signal that in BJ-hTERT was approximately 30% compared to untreated cells (CTR) but reached 75% in BJ-EHLT. These results suggested that the effect of the G4 ligands on cell viability could be due to the capability of these compounds to induce selective DNA damage in transformed cells.
Additionally, to assess the ability of the two biologically active compounds to target G4 structures in cells, a fluorescence in situ hybridization (FISH) assay was performed. The analysis of telomeric damage, evaluated by quantification of the co-localization spots of γH2AX with a fluorescent telomeric probe (Telomere Induced Foci, TIF), confirmed, on one hand, the capability of Chelidonine and Rotenone to elicit DNA damage in a dosedependent manner ( Figure 7A,D) and, on the other hand, demonstrated that a large part of this damage was telomere-located ( Figure 7B,D). Interestingly, both compounds exhibited a capability of inducing TIFs similar to pentacyclic acridine derivative RHPS4, a well-known G4 ligand used here as a positive control [48]. In particular, Chelidonine appeared as the best of the two ligands, reaching the highest percentage of TIF-positive cells with a mean number of five TIFs per nucleus ( Figure 7C,D).
Finally, the viability assays were extended to a human breast cancer cell line (MDA-MB-231; Figure 8A,B) over-expressing TRF2 (pTRF2; Figure 8C), a telomeric protein playing a key role in telomere protection [49,50]. Here, TRF2 over-expression was used as a tool aimed at definitively proving that the anti-tumor activity of Chelidonine and Rotenone was related to telomere targeting. Interestingly, under TRF2 over-expression, cells were preserved from the cytotoxic effect of the two compounds, as demonstrated by IC50 values that increased from 0.451 ± 0.066 to 0.934 ± 0.074 µM for Chelidonine and from 0.112 ± 0.006 to 0.306 ± 0.048 µM for Rotenone.
Altogether, our biological data led to the identification of two natural compounds that exhibited a potent DNA damage-mediated cytotoxic activity, selectively targeting telomeric G4s. Notably, between the two biologically active compounds, the highest selectivity of action on cancer over normal cells along with the highest specificity in telomere targeting and damaging was found for Chelidonine, in full agreement with its higher selectivity for G4 over duplex DNA than Rotenone, as determined by the G4-CPG assay and the CD-melting data.

Conclusions
Aiming at searching selective G4 ligands as putative candidate drugs for anticancer targeted therapies, a unique high-diversity library of natural compounds has been investigated. A molecular docking-based virtual screening identified 28 putative G4 ligands that were then evaluated by the G4-CPG experimental screening assay, whereby five molecules were confirmed as effective G4 ligands, i.e., Bulbocapnine, Chelidonine, Ibogaine, Rotenone and Vomicine. Then, CD and fluorescence spectroscopy indicated that the five investigated compounds can interact with G4s, also selectively stabilizing the G4 vs. duplex structures.
A detailed inspection of biophysical data revealed that the highest stabilizing effects and affinity on G4 over duplex structures were detected for Chelidonine, while the lowest ones were observed for Rotenone, in line with its lower G4 vs. duplex selectivity compared to the other four compounds as observed by the G4-CPG assay. For Chelidonine a stable binding to grooves of G4 structures was suggested by MD simulations along with its inability to firmly bind to duplex DNA, while Rotenone was proven to mainly target G4 flanking or loop residues, as well as to perform a rather specific interaction with the terminal end of the duplex DNA.
Moreover, Chelidonine and Rotenone, whose anti-tumor potential has been evaluated also in recent reports [51,52], were found to produce a potent anticancer activity mediated by their capability to bind and stabilize telomeric G4 structures. In particular, demonstrating that both compounds are effective at low µM doses and show selectivity for tumor cells, our results pave the way for the design of novel and even more effective synthetic analogs of these natural compounds that could find their application field in anticancer therapies. In this regard, NMR studies on the interaction of the selected natural compounds with G4 models are currently underway in our labs aiming at obtaining high-resolution structures of their complexes with G4s and conclusively establishing their ability to target the G4 grooves/loops.

Chemistry
All the tested compounds (1-28) are known structures belonging to the in-house library of natural products available from the Organic Chemistry Laboratory of the Department of Chemistry and Technology of Drugs of Sapienza University of Rome, Italy. The chemical identity of compounds was assessed by re-running NMR experiments and proved to be in agreement with the literature data reported below for each compound. The purity of all compounds, checked by reversed-phase High-Performance Liquid Chromatography (HPLC), was always higher than 95%.

Molecular Docking
The NMR structures of tel 26 and c-myc G4s and a 12-mer DNA duplex were retrieved from the Protein Data Bank [68] under the accession code 2JPZ, 1XAV and 1NAJ, respectively [29,30,69]. The first NMR model was extracted from PDBs and used as a rigid receptor in molecular docking simulations carried out with the AutoDock4.2 program [70]. According to previous studies, the ligand-binding site was centered in the G4s groove [36]. Specifically, in the first NMR model of the tel 26 structure, the binding site was centered in the groove formed by G4-G6, T8 and G22-T24 having 50 × 50 × 50 points dimension with a point-spacing of 0.375 Å. In the first NMR model of the c-myc structure, the ligandbinding site was centred on the loop formed by G2, A3, G5, G6, and G17-G19 having 60 × 60 × 60 points dimension with a point-spacing of 0.375 Å. In the first NMR model of the 12-mer DNA duplex, the binding site was centered in the mass center of A6-T19 and T7-A18 base pairs with a dimension of 60 × 100 × 60 points and a point-spacing of 0.375 Å to cover the entire surface of the unspecific duplex. Default AutoDock4.2 settings were used, and ten docking poses of each compound were stored. Ligands chemical diversity was evaluated by a Python-based clustering algorithm based on a combination of fingerprints and maximum common substructure search [40,71].

G4-CPG Assay
Nude CPG, tel 26 -, c-myc-and ds 27 -functionalized CPG supports were synthesized as previously reported [11,[34][35][36][37]. The stock solutions for each ligand were prepared by dissolving a known amount of the sample in pure DMSO, thus, obtaining 4 mM solutions (with the only exception of Bulbocapnine used as a 3.3 mM solution). A measured volume was taken from the stock solution to obtain a 600 µM ligand solution in a 50 mM KCl, 10% DMSO, 10% CH 3 CH 2 OH aq. solution. The detailed general procedure adopted for the assays is described as follows: weighed amounts of nude CPG or G4-and duplexfunctionalized CPG supports were left in contact with 300 µL of the 600 µM ligand solution in a polypropylene column equipped with a polytetrafluoroethylene frit, a stopcock and a cap. After incubation on a vibrating shaker for 4 min, each support was washed with defined volumes of the washing solution (50 mM KCl, 10% DMSO, 10% CH 3 CH 2 OH aq. solution) or the releasing solution (2.5 M CaCl 2 , 15% DMSO aq. solution or pure DMSO) and all the eluted fractions were separately analyzed by UV measurements. After treatment with the releasing solution, inducing G4s and duplex denaturation, the G4-and duplex-functionalized CPG supports were suspended in the washing solution and then subjected to annealing, by taking them at 75 • C for 5 min and then slowly cooling to room temperature.
The UV measurements were performed on a JASCO V-550 UV-vis spectrophotometer. A quartz cuvette with a path length of 1 cm was used. The UV quantification of the ligands was determined by measuring the absorbance relative to the λ max characteristic of each ligand and referring it to the corresponding calibration curves. The errors associated with the % of bound ligand were within ±2%.

Circular Dichroism
CD spectra were recorded in a quartz cuvette with a path length of 1 cm, on a Jasco J-715 spectropolarimeter equipped with a Peltier-type temperature control system (model PTC-348WI). The spectra were recorded at 20 • C in the range from 240-600 nm, with 2 s response, 200 nm/min scanning speed and 2.0 nm bandwidth, and were corrected by the subtraction of the background scan with buffer. All the spectra were averaged over three scans. The oligonucleotides d[(TTAGGG) 4 TT] (tel 26 ), d(TGGGGAGGGTGGGGAGGGTGGGGAA GGTGGGGA) (c-myc) and d(CGCGAATTCGCGTTTCGCGAATTCGCG) (ds 27 ) were purchased from Biomers as HPLC-purified compounds with a purity of >99%. The oligonucleotides were dissolved in a 5 mM KCl, 5 mM phosphate buffer, 5% DMSO (pH 7) or 0.5 mM KCl, 0.5 mM phosphate buffer, 5% DMSO (pH 7), thus, obtaining 2 µM solutions, which were then annealed by heating at 95 • C for 5 min, followed by slow cooling to room temperature. CD titrations were obtained by adding increasing amounts of the ligands (up to 10 molar equivalents, corresponding to a 20 µM solution in ligand) to tel 26 and c-myc G4s or ds 27 duplex. For the CD-melting experiments, the ellipticity was recorded at 290, 262 and 251 nm for tel 26 , c-myc and ds 27 systems, respectively, with a temperature scan rate of 0.5 • C/min, in the range from 10-95 • C.

Fluorescence Spectroscopy
Fluorescence spectra were recorded at 20 • C on HORIBA (Bensheim, Germany) JobinYvon Inc. FluoroMax ® -4 spectrofluorometer equipped with F-3004 Sample Heater/ Cooler Peltier Thermocouple Drive, by using a quartz cuvette with a 1 cm path length. For the fluorescence titration experiments with Bulbocapnine, Chelidonine and Ibogaine, excitation wavelengths of 307, 289 and 293 nm were used, respectively. The spectra were registered in the range from 315-600, 295-550 and 300-550 nm for Bulbocapnine, Chelidonine and Ibogaine, respectively. Titrations were carried out at a fixed concentration (2.0 µM) of ligand. Increasing amounts of tel 26 and c-myc G4s or ds 27 duplex (up to 10 µM concentration) were added from 120 µM annealed stock solutions of each DNA sample dissolved in a 5 mM KCl, 5 mM phosphate buffer, 5% DMSO (pH 7). After each addition, the system was allowed equilibrating 10 min before recording the spectra.
The fraction of bound molecules was calculated from the fluorescence intensity at 460, 330 and 356 nm for Bulbocapnine, Chelidonine and Ibogaine, and reported in a graph as a function of the DNA concentration.
The fraction of the bound ligand was determined using the equation: where Y, Y 0 and Y b are the values of fluorescence emission intensity at the maximum at each titrant concentration, at the initial and final state of the titration, respectively. The points were fitted with an independent and equivalent-sites model using the Origin 8.0 program. The equation of the independent and equivalent-sites model is as follows: where α is the mole fraction of ligand in the bound form, [L] 0 is the total ligand concentration, [DNA] is the added DNA concentration, n is the number of the equivalent and independent sites on the DNA structure and K b is the binding constant [44].

Molecular Dynamics Simulations
Docking complexes were solvated in a rectangular box of TIP3P-type water molecules [72] buffering 10 Å from the macromolecule surface. The total charge was neutralized by K + ions. The OL15 force field was used for DNA [73,74], while the GAFF2 force field was used for small molecules [75]. The ligands' partial charges were computed at the am1-bcc level. The Amber18 program was used to run MD simulations [76]. In agreement with previous works [77][78][79], the following protocol was implemented herein to generate robust MD trajectories: (i) the solvent was energy minimized for 500 steps with the steepest descent algorithm (SD) followed by 2500 steps with the conjugate gradient algorithm (CG); (ii) the solvated solute was energy minimized for 1000 steps with the SD followed by 5000 steps with the CG; (iii) heating from 0 to 300 K was achieved by the Langevin thermostat at constant volume for 1 ns; iv) system density was equilibrated for 1 ns using the Berendsen barostat at constant pressure [80]; (v) the system was preliminarily relaxed for 50 ns at constant pressure; vi) final production of MD trajectories was carried out for 500 ns on each system at constant pressure without positional restraints. Time step was 2 fs in all MD simulations, which were run using the GPU version of pmemd. Analysis of MD trajectories was carried out with the CPPTRAJ program [81], while the ligands' theoretical affinity was computed with the MMPBSA.py program [45].

Cells and Culture Conditions
Human fibroblasts (BJ-hTERT) were obtained by infecting primary BJ cells with a retrovirus carrying hTERT (Addgene plasmid #1773) resulting in a telomerase immortalized cell line, while BJ-EHLT were derived from the transformation of the same BJ cells with an hTERT and SV40 early region resulting in p53 and pRB silencing [82]. The human breast cancer cell line (MDA-MB-231) was purchased from ATCC. BJ-hTERT, BJ-EHLT and MDA-MB-231 were grown in Dulbecco's Modified Eagle Medium (D-MEM, Invitrogen Carlsbad, CA, USA) supplemented with 10% Fetal Bovine Serum (FBS), 2 mM L-glutamin and antibiotics at 37 • C in a 5% CO 2 -95% air atmosphere. Stable TRF2-overexpressing cells (pBabe-puromycTRF2) and the control counterpart (pBabe-puro-Empty) [83] were obtained by infecting the cells with amphotropic retroviruses generated into Phoenix packaging cells transfected with retroviral vectors, using the JetPEI reagent (Polyplus, New York, NY, USA), according to the manufacturer's instructions.

Viability Assay (Crystal Violet)
BJ-EHLT and BJ-hTERT cells were seeded in a 24-well plate at a density of 5 × 10 4 and 10 × 10 4 for well, respectively. After 24 h, cells were treated with the compounds at different doses for 72 h. Then, cells were washed twice in 1× phosphate-buffered saline (1× PBS) and fixed with 4% formaldehyde for 15 min at room temperature (RT). After washing, 500 µL of crystal violet staining solution (Sigma-Aldrich, St. Louis, MO, USA) was added to each well and incubated for 30 min at RT. Finally, the plates were rinsed twice with water, air-dried, and cell pellets were dissolved in 250 µL of acetic acid, 10% aqueous solution. 100 µL of each sample were transferred to a 96-well plate and the optical density was measured at 570 nm (OD570) with an ELISA reader (Thermo Scientific, Waltham, MA, USA). The average absorbance in each condition was used to calculate the viability expressed as a percentage of treated vs. untreated conditions. The half-maximal viability inhibitory concentration (IC50) was calculated by CalcuSyn Version 2.1.

Immunofluorescence (IF) and Fluorescence in situ Hybridization (FISH) Assays
For γH2AX fluorescent signal analysis, cells were fixed in 4% formaldehyde in 1× PBS for 15 min at RT and then permeabilized by treatment with 0.5% Triton X-100, 0.1% Na-Citrate (1× PBS) for 5 min at RT. Cells were blocked for 1 h in 3% BSA, 0. For immunofluorescence combined with DNA FISH, cells were fixed in 4% formaldehyde (1× PBS) followed by permeabilization with 0.1% Triton X-100 in (1× PBS) for 7 min at (RT). Cells were blocked for 1 h in 3% BSA (1× PBS) and incubated overnight with an anti-phospho-Histone H2AX antibody. After incubation with a secondary antibody, cells were fixed in 4% formaldehyde (1× PBS) and subjected to standard telomere DNA FISH. For quantitative analysis of γH2AX positivity, 300 cells on triplicate slides were analyzed. For TIFs quantification, 30 γH2AX-positive cells were scored at 63× magnification. Nuclei showing at least three colocalizations of TelC-Cy3 telomeric probe (Panagene, Daejeon, Korea) and γH2AX were considered as TIF-positive. Fluorescence signals were recorded by using a ZEISS LSM 880 confocal laser-scanning microscope by Carl Zeiss Ltd. (Oberkochen, Germany) Images were elaborated by ZEN Black 2.3 SP1 and γH2AX signal intensity was quantified using ImageJ 1.53e.

Clonogenic Assay
MDA-MB-231 cells were seeded in a 6-well plate at the clonogenic density of 500 cells/well and after 24 h cells were treated with DMSO (negative control) or the indicated doses of Chelidonine and Rotenone. After 24 h, fresh medium was replaced in each well and cells were allowed to grow for 10 days to form colonies. Then, cells were stained with 2% methylene blue in 50% ethanol and the number of colonies was counted. Surviving fractions were calculated as the ratio of the treated vs. negative control. IC50 values were calculated by CalcuSyn Version 2.1.

Statistical Analysis
The experiments were repeated three times and the obtained results are presented as mean ± standard deviation (S.D.). Statistical analysis was assessed by using the Student's t-test for unpaired data, and the values of p < 0.05 were considered statistically significant. Data analysis was performed with GraphPad Prism 6 (San Diego, CA, USA).

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/pharmaceutics13101611/s1, Table S1: Chemical structures, features and natural sources of the 28 natural compounds here investigated, Figure S1: CD spectra of 2 µM solutions of tel 26 27 duplex as a function of the DNA concentration. The black squares represent the experimental data; the red lines represent the best fit obtained using an independent and equivalent-sites model, Figure S11: Left: Fluorescence emission spectra obtained by adding increasing amounts of (A) tel 26 G4, (B) c-myc G4 and (C) ds 27 duplex to 2 µM solutions of Chelidonine. Arrows indicate the variation of fluorescence intensity on increasing DNA concentration. Right: Fluorescence intensity at 330 nm vs. concentration of (D) tel 26 G4, (E) c-myc G4 and (F) ds 27 duplex, Figure S12: Left: Fluorescence emission spectra obtained by adding increasing amounts of (A) tel 26