An In-Silico Pipeline for Rapid Screening of DNA Aptamers against Mycotoxins: The Case-Study of Fumonisin B1, Aflatoxin B1 and Ochratoxin A

Aptamers are single-stranded oligonucleotides selected by SELEX (Systematic Evolution of Ligands by EXponential Enrichment) able to discriminate target molecules with high affinity and specificity, even in the case of very closely related structures. Aptamers have been produced for several targets including small molecules like mycotoxins; however, the high affinity for their respective target molecules is a critical requirement. In the last decade, the screening through computational methods of aptamers for their affinity against specific targets has greatly increased and is becoming a commonly used procedure due to its convenience and low costs. This paper describes an in-silico approach for rapid screening of ten ssDNA aptamer sequences against fumonisin B1 (FB1, n = 3), aflatoxin B1 (AFB1, n = 2) and ochratoxin A (OTA, n = 5). Theoretical results were compared with those obtained by testing the same aptamers by fluorescent microscale thermophoresis and by magnetic beads assay for their binding affinity (KD) revealing a good agreement.


Introduction
Aptamers are nucleic acid biopolymers and represent a new generation of receptors that are selected in vitro by the Systematic Evolution of Ligands by EXponential Enrichment (SELEX) [1,2]. Aptamers are composed of single-stranded DNA (ssDNA) or RNA usually in the region of 35-100 nucleotides in length and have the ability to recognize different classes of targets such as mycotoxins, amino acids, antibiotics, pesticides, proteins and even whole cells with high specificity and affinity [3][4][5][6]. Compared to antibody generation, SELEX allows greater control over binding conditions and allows selection under non-physiological conditions. In addition to the potential for low production costs, benefits of aptamers with respect to antibodies include short time of development process, detection of non-immunogenic small molecules, stability during storage and improved tolerance to solvents or extreme pH or ionic strength.
The past few years have seen a dramatic increase in aptamer development and interest for diagnostic and therapeutic applications [7]. Recently, there has been an increase in demand for

Reagents and Apparatus
The ssDNA sequences of aptamers FB1_39 [25], FB1_39t3 [26], FB1_10 [27], AF_AB3 [20] and AF_APT1 [19], and all aptamers modified with 5 -6-carboxyfluorescein were purchased from Bio-Fab Research Srl (Rome, Italy) as lyophilised form and with HPLC purity. Then, each aptamer was solubilized in the respective buffer solution. The ssDNA aptamer sequences for FB1 and AFB1 and the relevant composition of buffer solutions are reported in Table 1. All chemicals for buffer solutions and all other purposes were purchased from Sigma-Aldrich (Milan, Italy). The Dynabeads™ M-270 Amine and the DynaMag™-2 Magnet were obtained from Invitrogen (Life Technologies, Monza, Italy). FB1 and AFB1 standards were purchased from Vinci-Biochem Srl (Florence, Italy).
Mycotoxins stock solutions were prepared by dissolving AFB1 powder in toluene:acetonitrile (90:10, v/v) at concentration of 0.5 mg/mL, FB1 powder in acetonitrile:water (1;1, v/v) at concentration of 1 mg/mL and OTA powder in toluene:acetic acid (99:1; v/v) at concentration of 0.5 mg/mL. Mycotoxins working solutions were prepared by drying appropriate aliquots of stock standard solutions and redissolving them in the respective buffer solutions at concentrations ranging from 0.1 nM to 200 µM. Information relevant to OTA aptamers and relative standard preparation are reported elsewhere [31,39].

Hardware and Software
For the aim of the present research, the following free computational software were chosen: SimRNA [40,41], RxDock [42], BioPython [43] and Open Babel [44]. SimRNA and RxDock where run on a 24 Intel ® Xeon ® E5-2640 core machine. All the computations were bounded by central processing unit (CPU) speed. The relevant timings can be found in the Supplemental Information (S.I.).

Prediction of the 3-D DNA-Aptamer Structures
Our protocol for the determination of the three-dimensional structure of ssDNA aptamers consisted of six consecutive steps as shown in the procedural scheme of Figure 1. First, the DNA sequence of residues was mutated to RNA by changing thymine to uracil nucleotides. Second, the RNA structure was annealed by means of SimRNA. This software can build an initial regular, essentially cyclical, arrangement of RNA coil from the sequence, even though, this is very far from realistic, therefore a long annealing process lets this evolve towards the stability region of the potential landscape. Third, a large Monte Carlo dynamic was performed at normal temperature (T = 1 in SimRNA units). The simulation consisted of 50 M steps, with frames being saved every 5000 steps. Fourth, the resulting 10,000 three-dimensional structures were classified into a small number of clusters. On the coarse-grained representation of SimRNA, a clustering parameter of 7 Å brings to 2-50 clusters per molecule. For each cluster, only the molecule with lowest energy was considered for further treatment. Fifth, the coarse-grained representation was translated back to an all-atom representation. We used the program SimRNA_trafl2pdbs, which was found bundled with SimRNA. The translation can bring to some deviation of the P-O5' bond between residues from the normal distance, sometimes as large as 0.5 Å, so that some software does not correctly recognize the chain as a single piece. Aside from problems with visualization software and interfaces to other computational programs, these deviations, that could be reduced with the help of molecular mechanics minimization, should not produce major effects on the length scale of docking potentials. Finally (sixth), the resulting RNA structures were mutated to DNA structures. A Python script, provided in the S.I., making use of the BioPython library [43], was written for the purpose. The program builds a standard Protein Data Bank (PDB) description of the residues, e.g., compatible with Gromacs expectations for the Assisted Model Building and Energy Refinement (AMBER) force fields. H12 GGGAGGACGAAGCGGAACCGGGTGTGGGTGCC TTGATCCAGGGAGTCTCAGAAGACACGCCCGACA Same as H8 [16] Polymers 2020, 12, x FOR PEER REVIEW 5 of 18 Figure 1. Scheme of the computational pipeline bringing from aptamer sequence to affinity estimate.

The Docking Approach
For each DNA molecule, the most stable configuration along with all similar configurations within 3 kcal/mol, were further processed and docked with the relevant toxin ligand molecule. To this purpose we used RxDock [42,45]. This software requires that the receptor and ligand be in specific formats, respectively, Sybil mol2 and MDL mol format. The conversion was achieved by means of Open Babel [44], which also adjusted the hydrogen atoms. RxDock defines the scanning room for docking as a sphere. However, defining a single sphere encompassing the entire receptor, tends to produce very large grids for ssDNA aptamers, which usually are not very compact. Indeed, the largest of our samples had for example an extension of almost 140 Å. To address this issue, a program was written that classifies the residue positions of a DNA configuration into a small number of clusters. For each of these a small sphere was obtained, typically of 15-25 Å, to which we added a 5 Å padding region.
RxDock allows customization of the docking through parameter files detailing the scoring function and each of the intermediate phases of the docking search, even providing some stock files for predefined purposes. For the present work, the complete docking protocol of RxDock was used, with the default parameters as detailed in dock.prm, coming with the RxDock package [42].
Docking trials were also carried out to test the selectivity of investigated aptamers towards mycotoxins different from their cognate target.

Docking Calculations and Relation to Experimental Data
Docking calculations are partially driven by random number generation. Replicating the calculation several times provides information about the effectiveness of space exploration and, hence, the quality of the resulting poses. This is basically reflected in a distribution of docking scores.
Our strategy for deciding the number of repetitions consisted in obtaining at least four replicas yielding docking scores within 3 kcal/mol of the lowest value. This means that in the worst and unrealistic case of a uniform docking score distribution, the probability to find a pose better than those already obtained by more than 3 kcal/mol was less than 1/8. However, for a more realistic bellshaped distribution, such probability is actually much lower.
In the case of OTA and AFB1 docking, the distribution was extremely tight and most repetitions yielded identical optimal results, therefore, the number of repetitions was fixed to 5. For FB1 the docking score distribution was extremely wide, therefore, to meet the mentioned criterion about 100 repetitions were needed.
The best docking score was then added to the aptamer configurational energy, yielding what we for simplicity call a ΔG 0 estimation (though actually more related to enthalpy), that was correlated with the experimental binding ΔG 0 . This is of course straightforwardly obtained from dissociation constant (KD) values through the relation ΔG = − ⁄ .

The Docking Approach
For each DNA molecule, the most stable configuration along with all similar configurations within 3 kcal/mol, were further processed and docked with the relevant toxin ligand molecule. To this purpose we used RxDock [42,45]. This software requires that the receptor and ligand be in specific formats, respectively, Sybil mol2 and MDL mol format. The conversion was achieved by means of Open Babel [44], which also adjusted the hydrogen atoms. RxDock defines the scanning room for docking as a sphere. However, defining a single sphere encompassing the entire receptor, tends to produce very large grids for ssDNA aptamers, which usually are not very compact. Indeed, the largest of our samples had for example an extension of almost 140 Å. To address this issue, a program was written that classifies the residue positions of a DNA configuration into a small number of clusters. For each of these a small sphere was obtained, typically of 15-25 Å, to which we added a 5 Å padding region.
RxDock allows customization of the docking through parameter files detailing the scoring function and each of the intermediate phases of the docking search, even providing some stock files for predefined purposes. For the present work, the complete docking protocol of RxDock was used, with the default parameters as detailed in dock.prm, coming with the RxDock package [42].
Docking trials were also carried out to test the selectivity of investigated aptamers towards mycotoxins different from their cognate target.

Docking Calculations and Relation to Experimental Data
Docking calculations are partially driven by random number generation. Replicating the calculation several times provides information about the effectiveness of space exploration and, hence, the quality of the resulting poses. This is basically reflected in a distribution of docking scores.
Our strategy for deciding the number of repetitions consisted in obtaining at least four replicas yielding docking scores within 3 kcal/mol of the lowest value. This means that in the worst and unrealistic case of a uniform docking score distribution, the probability to find a pose better than those already obtained by more than 3 kcal/mol was less than 1/8. However, for a more realistic bell-shaped distribution, such probability is actually much lower.
In the case of OTA and AFB1 docking, the distribution was extremely tight and most repetitions yielded identical optimal results, therefore, the number of repetitions was fixed to 5. For FB1 the docking score distribution was extremely wide, therefore, to meet the mentioned criterion about 100 repetitions were needed.
The best docking score was then added to the aptamer configurational energy, yielding what we for simplicity call a ∆G 0 estimation (though actually more related to enthalpy), that was correlated with the experimental binding ∆G 0 . This is of course straightforwardly obtained from dissociation constant (K D ) values through the relation ∆G 0 = −RTln(K D /M).

Fluorescent Microscale Thermophoresis (MST)
Binding assays by fluorescent Microscale Thermophoresis (MST) for FB1 and AFB1 aptamers were carried out by the company 2Bind, GmbH (Regensburg, Germany), according to the procedure described by McKeague et al. [39]. Briefly, 5 µL of each serial dilution of standard mycotoxin (from 1.464 nM to 48 µM) in the appropriate binding buffer were mixed with 5 µL of 10 nM of fluorescently labelled aptamer. The final reaction mixture, containing a constant concentration of fluorescent aptamer (5 nM) and a variable mycotoxin standard solution (ranging from 0.732 nM to 24 µM) was filled in capillaries. Samples were analyzed on a Monolith NT at 25 • C, with 35% LED power and 80% Laser power. The binding between aptamer and mycotoxin was monitored by measuring the changes in the thermophoretic behavior through a temperature gradient. For each aptamer two independent experiments were carried out at each tested concentration. The dissociation constant K D value (highest affinity, lowest K D ) was obtained by fitting the binding curve with the fraction of fluorescent aptamer that formed the complex, calculated for each independent experiment (n = 2) from the law of mass action [46]. Details on the MST assay for OTA aptamers are reported in a paper by McKeague et al. [39].

Magnetic Beads (MBs) Assay
The functionalization of magnetic beads with FB1 was performed according to the procedure described by Frost et al. [26]. Briefly, an aliquot (450 µL) of amine-activated Dynabeads ® (10 8 beads, 30 mg/mL) was washed three times with PBS buffer (50 µL) and separated on the magnet each time. Then, the beads were re-suspended with 5% glutaraldehyde (Sigma Aldrich) in PBS buffer and incubated (for 2 h at room temperature) under slow stirring. After three washes with PBS buffer, the beads were re-suspended into 450 µL of FB1 standard (20 µM) in PBS buffer solution, and incubated (at room temperature) under slow stirring overnight. At the end of the process, unbound NH 2 groups, still present on the bead surface, were capped by incubation with a solution of Sulfo-NHS acetate (Thermo Scientific) in buffer (for 2 h at room temperature) under slow stirring. The beads were washed other three times with PBS buffer and stored at 4 • C till use. All the washing solutions from the coupling were collected and analyzed by HPLC for unbound FB1 quantification according to De Girolamo et al. [47]. The analyses indicated that more than 99% of FB1 was conjugated to the beads.
The functionalization of magnetic beads with AFB1 was performed according to the manufacturer's protocol (Invitrogen by Life Technologies) and the coupling was achieved by Schiff base (imine) formation and reductive amination. Briefly, an aliquot (300 µL) of amine-activated Dynabeads ® (10 8 beads) were washed three times with 0.1 M sodium borate buffer, pH 9.4 (300 µL) and separated on the magnet each time. The washed beads were mixed with 300 µL of AFB1 standard solution (180 µM) in 0.1 M sodium borate buffer; then, an aliquot (6 µL) of 5 M sodium cyanoborohydride in 1 M NaOH was added to the reaction mixture, vortexed and incubated (for 2 h at room temperature) under slow stirring. After removal of supernatant, an aliquot (600 µL) of 0.1 M ethanolamine (pH 7.4) was added to the reaction mixture and incubated (for 15 min at room temperature) by gentle mixing. The coated beads were then washed three times with 0.5% BSA, 0.01% Tween ® -20, and 0.02% sodium azide, and beads were stored in buffer (100 mM NaCl, 20 mM Tris-HCl (pH 7.6), 2 mM MgCl 2 , 5 mM KCl, 1 mM CaCl 2 at 4 • C) at 4 • C until use. All the washing solutions collected after the functionalization reactions were collected and analyzed by fluorescence spectroscopy (LS-55, Perkin Elmer, Waltham, MA, USA) for AFB1 quantification using a standard 10 mm path length quartz cuvette. Spectra were acquired in the 350-600 nm spectral range with λ ex set at 494 nm and λ em set at 520 nm. It was calculated that more than 99% of the mycotoxin was conjugated to the beads.
The affinity of the selected aptamers for their target mycotoxin was determined by incubating increasing concentrations of 5 'FAM labelled aptamer (100 µL) prepared in the respective buffer with a fixed amount of mycotoxin-coupled beads (25 µL). In particular, concentrations were between 0.5-100 nM for aptamer FB1_39, between 5-2000 nM for aptamer FB139_t3, between 20-2000 nM for aptamer FB1_10 and between 0.5-250 nM for AF_Apt1, between 0.5-4000 nM for aptamer AF_AB3. After a 1-h binding reaction (at room temperature) the unbound aptamers were removed by several washing steps with binding buffer. The amount of DNA aptamer eluted from the beads was determined by fluorescence measurements with λ ex set at 494 nm and λ em set at 520 nm and calculations using a calibration plot. For each aptamer two independent experiments were carried out. The value of K D was calculated by applying a nonlinear regression equation using the function of "receptor for a non-specific binding site" (SigmaPlot v. 12.3). Details on the MBs assay for OTA aptamers are reported in the paper by McKeague et al. [31].

Results and Discussion
In the present study three families of aptamers for three target mycotoxins, namely, fumonisin B1 (FB1), aflatoxin B1 (AFB1) and ochratoxin A (OTA), were screened by using an in-silico pipeline approach for their binding affinity towards target mycotoxins. The reliability of the proposed approach was verified by comparing theoretical results with those obtained by using two experimental approaches, i.e., Fluorescent Microscale Thermophoresis (MST) and Magnetic Beads (MBs) assays.

In-Silico Binding Affinity of Mycotoxins to Aptamers
The strategy presented herein for gauging the binding affinity of an aptamer-toxin couple consists of obtaining a set of possible conformations for the aptamer, within a reasonable range of conformational energies. Then, the toxin is docked to the aptamer to obtain a docking score that summed to the conformational energy provided an estimate of the binding free energy.
As much as the first task is concerned, i.e., obtaining the three-dimensional aptamer configurations, our six-step procedure may look somewhat convoluted. The problem is that, to our knowledge, there is no software that straightforwardly addresses this problem for ssDNA. Contrarily, there is a good choice of software specialized in making previsions of structure for RNA fragments. Therefore, we were somewhat obliged to take an indirect route, that is, studying a corresponding RNA fragment and then rigidly mutating this to a DNA fragment by leaving fixed all the atoms in common. Zhang et al. [37] opted for a similar strategy, but used Mfold and RNA Composer for the realization of RNA tertiary structure. Our choice, SimRNA, though certainly more computationally expensive, presents a distinct advantage.
The main duty of SimRNA is to perform a Monte Carlo dynamic of a simplified, coarse grained representation of an RNA fragment. Such coarse-grained representation consists of 5-6 freedom degrees per residue, dealing with nucleotides as essentially rigid entities. The energy landscape is provided by a model potential, the most important terms of which were written so that simulations reproduce geometrical relationships observed in nature. These were extracted from a curated database of experimental well-resolved structures. The main effect of this approach is that SimRNA provides an evolving view of the molecule rather than a single geometry. The molecule movie is then characterized with clustering techniques, by means of which characteristic configurations are extracted. While all this results in extra work, it may produce more reliable results, considering the fact that the subsequent docking is conducted holding the aptamer rigid. For example, for aptamer FB1_39, we considered four different configurations. According to SimRNA, these configurations are pretty much equally plausible, can be found with similar frequency during the simulation and, possibly, are all similarly present in a real solution.
As for most docking software, the RxDock scoring function is a sum of interaction energies between receptor and ligand. These are expressed in kcal/mol and are meant to mimic the enthalpy of the binding process: the lower the value, the better the binding. It must be noted that, notwithstanding the physics-based recipe behind it, the primary purpose of the scoring function for RxDock is to provide a ranking of the docking poses and not to provide an energetic estimate of the binding process. In the same way, the main role of the SimRNA potential for SimRNA is to assess a relative probability of a configuration compared to another, the authors of SimRNA are explicit about this. Therefore, while our estimate of the binding affinity as E dock + E SimRNA is somewhat natural and obligatory, its correspondence to ∆G 0 cannot be taken for granted but needs to be verified by means of a sort of experimental verification. Table 2 reports a set of predictions for each of the selected aptamers towards target mycotoxins, the relevant docking score values and the DNA configurational energy. In general, the pool of DNA configurations was restrained to those within 3 kcal/mol of the most stable one, but often only one configuration survives this selection criterion. However, in the case of FB1_10, we decided to include more configurations. For OTA and FB1, the docking energy reported pertains to the best of five different docking runs. This repetition was deemed necessary because, depending on the exhaustiveness of the configuration space sampling, docking calculations are not necessarily deterministic; indeed, docking calculations somewhat rely on random number generation. However, in the case of OTA and FB1, the resulting docking calculation was highly reproducible, with small dispersion and rare cases of badly guessed poses. Figure 2 shows a histogram of the docking score values for more than one hundred calculations of the FB1_10-FB1 system in a docking sphere of 15.1 Å, with identical input. The dispersion was very high, so that only five poses were within the lowest 5 kcal/mol. Therefore, one can expect so many samples to be necessary to obtain reliable results for systems like fumonisins. This is probably due to the large number of free dihedrals and therefore configurational space of such molecules. Incidentally, each run took 2-3 times more CPU time for these flexible molecules.   [48], which is partially supported in our modelling ( Figure 3). Note that RNA prefers to adopt an Aform duplex which would account for the differences in our model. In contrast, in the case of OTA aptamers, while the modelling of aptamer 1.12.2 predicts a stem-loop with a short single-stranded overhang ( Figure 5), circular dichroism experiments show the characteristics of an anti-parallel Gquadruplex structure [49]. This discrepancy could be due to the difficulty in modelling G-quadruplex structures by SimRNA and should be kept under consideration when comparing binding affinity results.  A large variety of binding possibilities were obtained, with hydrogen bonding, hydrophobic repulsion and other kinds of electrostatic interactions being involved with differing number and intensity and no road to a simple or uniform interpretation of the interaction energy. Fumonisins, however, evidenced a higher predisposition for bonding to the aptamer backbone. Notably, the majority of these aptamer structures compare well with what has been predicted for the sequences. For example, circular dichroism experiments on FB1_39 show indications of a B form duplex structure [48], which is partially supported in our modelling ( Figure 3). Note that RNA prefers to adopt an A-form duplex which would account for the differences in our model. In contrast, in the case of OTA aptamers, while the modelling of aptamer 1.12.2 predicts a stem-loop with a short single-stranded overhang ( Figure 5), circular dichroism experiments show the characteristics of an anti-parallel G-quadruplex structure [49]. This discrepancy could be due to the difficulty in modelling G-quadruplex structures by SimRNA and should be kept under consideration when comparing binding affinity results.
Aptamer selectivity is also an important characteristic that needs to be understood before a sequence can be determined to be fit for purpose. The queried aptamer systems should all be expected to show limited cross-reactivity to mycotoxins other than their cognate target. The docking of AFB1 with all the OTA-binding aptamers revealed docking energies that were consistently less stable for AFB1, as expected (Table S2). In particular, A08 and A08min showed a destabilization of 4.28 and 6.95 kcal/mol, respectively, when the incorrect target was docked. Considering that a difference of 1 kcal/mol equates to a factor 5.4 in the binding constant, the selectivity of the OTA-binding aptamers is confirmed. In contrast, when OTA was docked onto the AFB1-binding sequences, similar or improved docking scores were observed (Table S1). This finding is in agreement with previous studies reporting a comparable or slightly lower selectivity towards other target mycotoxins for aptamers specific for AFB1 [18,37]. Further computational investigations will be carried out to perform other cross-reactivity studies to predict the affinity of the investigated aptamers in situations which may arise in multi-mycotoxin-contaminated samples. Further investigations will focus also on verifying if the AF_AB3 and AF_APT1 aptamers use a different binding site for OTA compared to that of AFB1.      Aptamer selectivity is also an important characteristic that needs to be understood before a sequence can be determined to be fit for purpose. The queried aptamer systems should all be expected to show limited cross-reactivity to mycotoxins other than their cognate target. The docking of AFB1 with all the OTA-binding aptamers revealed docking energies that were consistently less stable for AFB1, as expected (Table S2). In particular, A08 and A08min showed a destabilization of 4.28 and 6.95 kcal/mol, respectively, when the incorrect target was docked. Considering that a difference of 1 kcal/mol equates to a factor 5.4 in the binding constant, the selectivity of the OTA-binding aptamers is confirmed. In contrast, when OTA was docked onto the AFB1-binding sequences, similar or

Experimental Binding Affinity of Aptamers to Mycotoxins
Fluorescent Microscale Thermophoresis (MST) is a recently described method that leverages the physical phenomenon of molecular movement within temperature gradients. Each molecule (or aptamer) has a very specific "thermophoresis" that is dependent on the size, charge and hydration shell; therefore, upon binding to a target molecule, at least one of these parameters will be altered and can be measured. Because the movement of the interacting partners can be monitored with fluorescence, in-solution binding information can be obtained [50]. MST assay was chosen as an experimental approach because it is a powerful tool to characterize aptamer-target interactions independent in a highly sensitive cost-and time-efficient performance and with a free choice of buffers. Furthermore, the low sample consumption and the flexible assay setup are great advantages of this in-solution method [50]. Binding assays results by MST indicated a good binding of aptamers FB1_39 and FB1_10 towards FB1 and of aptamer AF_AB3 towards AFB1, with K D values lower than 200 nM, even though a certain variability of results was observed (up to 85%). This variability was probably related to the poor solubility of the mycotoxins in the buffer solutions at the tested concentrations. The aptamer AF_APT1 did not show any binding towards AFB1 (Table 3). In the case of OTA aptamers, all of them showed a good binding affinity towards the target mycotoxin in the nanomolar range, with A08min showing the highest affinity as previously reported [39] (Table 3).
Magnetic beads (MBs) assay is a fast and easy-to-use tool for isolation of analytes from complex matrices thanks to the colloidal stability of magnetic nanoparticles, homogenous size distribution, high and uniform magnetite content and presence of surface functional groups [30]. Results of the MBs assay confirmed the good binding affinity of aptamer FB1_39 towards FB1, while the other two aptamers, i.e., FB1_39t3 and FB1_10, did not show any binding towards FB1. Conversely to results obtained from MST assay, a good binding affinity towards AFB1 was obtained for AF_APT1, while aptamer AF_AB3 did not show any binding towards AFB1. All the aptamers towards OTA showed K D values comparable to those obtained by MST [31,39] (Table 3).
The different results observed herein between the aptamers belonging to the same family could be explained by the typology of the assays. Indeed, a previous study carried out by McKeague et al. [44] demonstrated that a given aptamer may have limitations or opportunities that make it better suited for certain molecular recognition schemes or applications. In the present study, in the MST assay both the target mycotoxin and the aptamer were free in solution, while the MBs assay was based on the interaction between free aptamer in-solution and bounded mycotoxin on the beads' surface. The fact that the binding studies carried out by MBs did not show any binding for three out five aptamers as compared to MST assay indicated that the labelling at 5 end of the aptamer may have abolished the binding ability of the aptamer (Table 3). Furthermore, when the target mycotoxin is bounded on a surface, the assays may suffer from nonspecific bindings between the aptamers and the surface, thus reducing their functionality in the investigated application. These results were in agreement with those reported by McKeague et al. [31,44] that performed a comparison of aptamer binding assays for a family of seven different aptamers towards the OTA using eight different assays, including MBs. As important differences beyond binding affinity between the three families of aptamers were observed, they concluded that aptamer binding varied when either the target or aptamer is used in solution vs. immobilized, and the sensitivity of each technique affected the apparent affinity [31,39]. A similar variability in terms of K D values calculated by using three different experimental methods was also reported by Mousivand et al. [38] for a family of six aptamers against AFB1. Altogether these results indicate that to screen aptamers it is important to implement high-throughput methods in order to screen more candidate ones. In general, methods allowing aptamer-target binding occurring free in solution, as the MST assay described herein, are the preferred ones because they remove any contribution from matrix binding responsible for reducing aptamer functionality. On the other hand, colorimetric assays, like the MBs one, are widely applied in quantitative detection of mycotoxins, especially in the real-sample applications [29,30,32].

Comparison of In-Silico and Experimental Binding Affinity of Aptamers to Mycotoxins
Characterization of aptamer binding is critical for studying aptamer molecular recognition and integrating them into diverse applications. In the current decade, the screening of aptamers through computational docking and molecular dynamics study has greatly increased and the application for computational prediction of the aptamer affinity against its specific target is becoming a commonly used procedure [35]. Indeed, the limitation in terms of costs of the high-throughput screen of aptamers through experimental trials are overcome by in-silico methods that allow to rapidly and reliably screen aptamer candidates.
To compare the in-silico prediction results with the experimental ones for each aptamer-toxin couple, the K D results, calculated by using MST and MBs assays and expressed as nmol/L values, were translated to free energy ones by means of the relation ∆G 0 = −RTln(K D ) formula, where T is held at 298 • K. This, as discussed above, is then compared to an estimate of ∆G 0 obtained as E dock + E conf . Figure 6 is a visual representation of the comparison, overlapping for each aptamer the experimental results from MST and MBs tests and the in-silico prediction. Vertical dashed lines were used to describe those cases where K D and the corresponding binding energy are above the technique measurement limit. The agreement between simulation and experimental data is almost quantitative, however, predictions for FB1, and in particular its binding to FB1_10 deserve a few considerations. In general, such model calculations as those in the present paper, mostly presuppose that the system is free in solution. This is not granted under all respects. For example, as said above, the SimRNA potential is derived from experimental structural data, most of which certainly do not refer to dilute RNA solutions. In most cases, the experimental MST and MBs affinities are similar, making the comparison more straightforward. However, as explained above, for AF_AB3, FB1_10 and FB1_39T3, MBs is not able to detect any binding. Moreover, in the case of AF_APT1, MST gives a negative response. In these limited four out of ten cases, we cannot make any quantitative comparison. However, the calculated prediction for FB1_10 is certainly suspect. In the case of FB1_10 we also relieved the sampling constraints on conformations to extend the study to more improbable ones, but could not detect any effective binding. More in general, a more accurate benchmark is needed for molecules like fumonisins, both to assess the reliability of calculations and to improve the efficiency of the simulation. Basically, when the MST and MB assays give similar or strictly correlated affinity values, as is the case for OTA, one can avoid pondering whether the in-silico data are better compared to the first or the latter set of experimental data, a question which is not easily addressed because of the sources of the SimRNA potential. On the other hand, AFB1 comes with two aptamer samples, one with MB good affinity and very poor MST affinity, the other doing the exact opposite. MBs is not able to detect any binding. Moreover, in the case of AF_APT1, MST gives a negative response. In these limited four out of ten cases, we cannot make any quantitative comparison. However, the calculated prediction for FB1_10 is certainly suspect. In the case of FB1_10 we also relieved the sampling constraints on conformations to extend the study to more improbable ones, but could not detect any effective binding. More in general, a more accurate benchmark is needed for molecules like fumonisins, both to assess the reliability of calculations and to improve the efficiency of the simulation. Basically, when the MST and MB assays give similar or strictly correlated affinity values, as is the case for OTA, one can avoid pondering whether the in-silico data are better compared to the first or the latter set of experimental data, a question which is not easily addressed because of the sources of the SimRNA potential. On the other hand, AFB1 comes with two aptamer samples, one with MB good affinity and very poor MST affinity, the other doing the exact opposite. A possible explanation of the more homogeneous behavior of the five OTA aptamers with the two different experimental approaches could be related to their secondary structure. Indeed, the sequence analysis of OTA aptamers by using the "m-Fold" program revealed in all the five aptamers two conserved consensus sequences. The first sequence (GGGTGTGGG) was located in the steam region, which stabilizes the whole structure due to high GC contents, while the second one (AGGGAGT) was located in the in the single-strand terminal loop. Both these two regions were most likely responsible for the binding with OTA independently from the surrounding in which the recognition takes place [16]. In the case of AFB1 and FB1 aptamers investigated herein, none of these or other conservative sequences were observed in their structure, as well as, a lower content of GC was present. Both these factors could be responsible of the minor stability of these families of aptamers in the tested conditions and then, to the greater lack of homogeneity of the results. A further confirmation of these results comes from the numbers of applications reported in literature for these three families of aptamers. In particular, the most numerous applications to food samples A possible explanation of the more homogeneous behavior of the five OTA aptamers with the two different experimental approaches could be related to their secondary structure. Indeed, the sequence analysis of OTA aptamers by using the "m-Fold" program revealed in all the five aptamers two conserved consensus sequences. The first sequence (GGGTGTGGG) was located in the steam region, which stabilizes the whole structure due to high GC contents, while the second one (AGGGAGT) was located in the in the single-strand terminal loop. Both these two regions were most likely responsible for the binding with OTA independently from the surrounding in which the recognition takes place [16]. In the case of AFB1 and FB1 aptamers investigated herein, none of these or other conservative sequences were observed in their structure, as well as, a lower content of GC was present. Both these factors could be responsible of the minor stability of these families of aptamers in the tested conditions and then, to the greater lack of homogeneity of the results. A further confirmation of these results comes from the numbers of applications reported in literature for these three families of aptamers. In particular, the most numerous applications to food samples contaminated with mycotoxins were described for OTA aptamers followed by AFB1 aptamers and to a lesser extent for FB1 aptamers [29,30,32], thus suggesting the major versatility of OTA and AFB1 aptamers compared to the other ones.

Conclusions
In the present study we propose for the first time an in-silico pipeline approach for the rapid screening of different ssDNA aptamers against three mycotoxins such as fumonisin B1, aflatoxin B1 and ochratoxin A. Furthermore, the comparison of the outputs of the theoretical screening with experimental ones was a useful strategy to validate the theoretical results. The proposed in-silico method for affinity prediction is conveniently cheap, requiring modest computational resources and being easily automatable to cope with large aptamer databases. With the current parameters, in general, good agreement of theoretical results with those obtained with experimental assays was observed and in particular with MST assay in which both aptamer and target mycotoxin are free in solution. Furthermore, it worked far more reliably and efficiently for aflatoxin B1 and ochratoxin A and to a lesser extent with fumonisin B1, and therefore, plausibly, with similar molecules, endowed with a large number of free dihedral angles. It is, however, reasonable that a strict adaptation of the RxDock docking protocol could provide good results also for this class of molecules.
The consistency of the in-silico results with the experimental findings confirmed that this strategy is a promising tool to be used for the screening of aptamer affinity towards small molecules like mycotoxins, even though the selectivity of aptamers observed by computational methods should be confirmed by experimental tools.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4360/12/12/2983/s1, Table S1: In-silico selectivity tests for aflatoxin B1 (AFB1) aptamers towards ochratoxin A (OTA) and comparison with docking calculations for the target AFB1. Table S2: In-silico selectivity tests for ochratoxin A (OTA) aptamers towards aflatoxin B1 (AFB1) and comparison with docking calculations for the target OTA. Table S3: Timings (s). Funding: This research has been supported in part by the MYCOKEY project "Integrated and innovative key actions for mycotoxin management in the food and feed chain" (H2020-Grant Agreement No. 678781).

Conflicts of Interest:
The authors declare no conflict of interest.