Exploring PAZ/3′-overhang interaction to improve siRNA specificity. A combined experimental and modeling study

A direct connection between the PAZ/3′-overhang binding affinity and the siRNA potency and specificity is defined through complementary experimental and computational results.


RNA synthesis and siRNA preparation
Modified RNA strands were synthesized on the 0.2-μmol scale using LV200 polystyrene supports. 3'-modified RNAs were synthesized on the 1-μmol scale using CPG functionalized with β-L-Thymidine, 1 L-threoninol-Acridine, 2 2-deoxyRibitol 3 and GNA-Thymine 4 units as solid supports. All oligonucleotides were synthesized on an Applied Biosystems 3400 synthesizer using commercially available reagents and 2'-O-TBDMS-5'-O-DMTprotected phosphoramidites (A Bz , G dmf , C Ac and U) in DMT-ON mode. The coupling time was 15 min and the coupling yields of natural and modified phosphoramidites were >97%. Unmodified (wt, 19, SCR and BLT) and phosphorothioate (PS) oligonucleotides (Table S1 and Table S2) were purchased from Sigma-Aldrich. SiRNA sequences described by Terrazas et al. 5 were used to design siRNA duplexes against Renilla gene (Table S2) and non-active scrambled sequences (SCR and BLT) ( Table S3). Solid supports were treated at 55 °C for 1 h with 1.5 mL of NH 3 solution (33%) and 0.5 mL of ethanol. Then, the oligonucleotides were purified using Glen-Pack TM RNA purification cartridges (Glen Research), quantified by absorption at 260 nm and analyzed by HLPC Column: Nucleosil 120-10 C18 column (250 × 4 mm). Solvent A: 5% ACN in 0.1 M aqueous TEAAc (pH = 7) and solvent B: 70% ACN in 0.1 M aqueous TEAA (pH = 7). Flow rate: 3 mL/min. Conditions: 20 min linear gradient from 15% to 80% B and 5 min 80% B. Finally, the oligonucleotides were confirmed by MALDI mass spectrometry (Table S1). SiRNA duplexes were prepared by annealing equimolar ratios of sense and antisense siRNA strands were in siRNA suspension buffer (100 mM KOAc, 30 mM HEPES-KOH, 2 mM MgCl 2 , pH 7.4) to a final concentration of 20 μM, duplexes were then heated at 95 °C for 5 min and slowly cooled to 4 °C. Table S1. RNA sequences used for the formation of siRNA duplexes.

Psicheck2 AS and SS Reporters
Psicheck2 AS and SS vectors have been previously described. 8

Transfection and Luciferase Assay
24 H prior transfection, HeLa cells were seeded in 24-well plates at density of 1x10 5 cells well -1 ; MEF and TEF TRBP -/were plated in 24-well plates at density of 0.8x10 5 cells well -1 . At time of transfection, the growth medium was changed with fresh DMEM +10% FBS without Pen/Strep (500 µl well -1 ). Renilla expression was measured on lysates collected 24 h post-transfection using Dual-luciferase reporter assay system (Promega) and Glomax luminometer (Promega). The ratios of Renilla luciferase (hRluc) to Photinus luciferase (hluc+) protein activities were normalized to mock transfection and the mock activity was set as 100%.

Isolation of RNA and RT-qPCR
Total RNA was extracted from MEF and MEF Ago2 -/with TRIzol reagent (Invitrogen). Then, isolated RNA was quantified by NanoDrop (Thermo Scientific) and 2.5 µg of each RNA sample was treated with DNase I (RNase free) (New England Biolabs) following the manufacturer's instructions. Reverse transcription reaction was performed on 0.5 µg of total RNA using Oligo (dT) Primer and Revertaid H minus RT enzyme (Thermo Scientific) according to the manufacturer's instructions. cDNA was subsequently diluted 4 times UltraPure DNase/RNase-Free distilled water (Invitrogen) and 1 µL of resulting cDNA to the Maxima Sybr Green/ROX qPCR Master Mix (Thermo Scientific). Real-time qPCR was accomplished in a total volume of 20 µL, following the manufacturer's instructions. The reference gene GADPH was used as the internal control. Renilla silencing was calculated and represented as 2 −∆∆Ct method. All primer pairs were purchased from Sigma-Aldrich and Primer-Blast was used as the primer designing tool.

Statistical Analysis
Statistical analysis was performed using GraphPad Prism software (GraphPad, San Diego, CA, USA). IC 50 determination was performed using non-linear regression analysis (log [inhibitor] vs. normalized response). Table S3. Sequences of siRNA duplexes used in competition assay.

Forward
Reverse

MD simulations 1.2.1 Starting structures
The starting structures used to describe the interaction between PAZ domain and 3'-overhang of siRNA were obtained from Protein Data Bank (PDB), selecting (i) the structure of PAZ from human Ago eIF2c1 (Set I) bound to the double stranded 9-mer siRNA, (PDB code: 1SI2) 9 and (ii) the structure of human Ago2 (Set II) bound to 4-mer single stranded miRNA (PDB code: 4F3T). 10 From the latter, the atoms corresponding to PAZ domain and the described four terminal nucleotides belonging to 3' end of siRNA were selected.

Molecular systems construction
In the preparation of PAZ, the lacking residues were added to the initial crystal structure and optimized through Modeller9v15. 11 The N and C termini were capped with neutral acetyl and N-methyl group, respectively. Hydrogen atoms were added using Leap module of AMBER14. The overall charge of PAZ domain was set to correspond to pH=7. In comparison to experimental counterpart, the last two residues positioned in the 3'-end terminal of the RNA x-ray crystal structure, were substituted by the natural nucleotide thymidine or by chemically modified nucleotides (see Table S1). The remaining nucleotide sequence match the one used experimentally. Similarly to protein construction, the Leap module was used to modify siRNA molecules. The new modified residues were parameterized by using the R.E.D.D. server 12 which automates the calculation of RESP charges at the HF/6-31G(d) level of theory, consistent with the AMBER14SB force field. Specifically for 2'-O-methyl and phosphorothioate modifications additional force fields have been used. 13,14 Both phosphorothioate residues were simulated considering R p configuration. All crystal water molecules surrounding PAZ/siRNA complex were considered during simulations.

Molecular dynamic simulations
PAZ/siRNA complexes were parameterized using AMBER's ff14SB force field, which collects the recent updates of the backbone and the χ torsion parameters. [15][16][17] The original complex conformations were initially minimized for 2500 cycles using the steepest-descent algorithm in implicit Generalized Born (GB) solvent, and the final geometries were then solvated in an octahedral box. All systems were solvated using a TIP3P waters 18 maintaining a minimum distance between the complex and the box wall of approximately 12 Å. For net charge neutralization sodium ions were added. In order to ensure conditions close to the experimental ones, Na + and Clions were added to the complex to attain a concentration of approximately 150 mM, and their positions were randomized with Amber's Ptraj program. 19 The number of solvent molecules added to each system was made to be identical within related crystal structures, allowing the energetic comparisons between siRNA modifications. All the solvated systems were initially minimized in two steps. During initial step, the atom positions of the siRNA and PAZ molecules were restrained applying 25 Kcal mol -1 Å -2 of positional restraints, using 1000 steps for the steepest-descent algorithm, followed by a conjugate-gradient minimization involving another 1000 steps. In the second step, the entire system was minimized without restraints using the same algorithms as the previous step but with 2500 steps for each. Subsequently, the systems were heated from 0 to 298 K using Langevin thermostat over 1 ns with collision frequency of 2 ps -1 and 5 Kcal mol -1 Å -2 of positional restraints on all of the solute atoms. Final equilibration at constant pressure and at 300 K was carried out for 2 ns using Langevin thermostat, and using 0.5 Kcal mol -1 Å -2 of positional restraints on the solute atoms. Finally, unrestrained systems were simulated under NPT ensemble conditions using periodic boundary conditions and particle mesh Ewald (PME) 20 for long range charge interactions. For all heating, equilibration and production runs, the time step was set to 2 fs, and short-range interactions were set to a cutoff of 10 Å. Bonds involving hydrogen were held fixed using the SHAKE algorithm. 21 System coordinates were recorded every 2 picoseconds. Six short simulations were produced from the crystal structure with different initial velocities, via random number generator seeds using the option ig = -1 of AMBER. Production simulations were carried out using CPU version of PMEMD program in AMBER14 22 on Navigator computing cluster of LCA-UC. 23 VMD was used for visualization. 24 RSMD, RMSF and H-bond analysis were carried out in AMBER's Ptraj program. Data from H-bond analysis were compiled and processed using an in-house script developed by the authors using Fortran 77 and R programming. The hydrogen bonds were considered relevant when the distance between acceptor-donor is equal or less than 3.5 Å and the angle cutoff of equal or higher than 135°. The last 15 ns of the trajectory were used for RMSF and H-bonds sampling.

Energy analysis using MM-GBSA
Generalized Born and surface area continuum solvation (MM-GBSA) method 25 is an end-state calculation frequently employed to estimate the binding free energy in a noncovalently bound protein-ligand complex. To obtain a detailed description of the PAZ/siRNA interaction, MM-GBSA was used for implicit solvation using Debye-Hückel screening. Following the MM-GBSA approach, the binding free energies were calculated according to: where each term is estimated from snapshots, denoted as ( ), taken from MD trajectories. For each snapshot, the free energy of a state is estimated by: where the first two terms, ( ) and / ( ) are calculated from internal (bond, angle, dihedrals), electrostatic and van der Waals energies. The third term of the equation, ( ) comprises the contribution of polar and non-polar terms in which the polar contribution is typically obtained by using the generalized Born (GB) model whereas the non-polar tern is estimated from the solvent-accessible surface area (SASA) determined with the LCPO method. 26 Finally, the last term, ( ) corresponds to the conformational entropy to the binding. MM-GBSA analyses were made with MMPBSA.py program 27 in AMBER14 using the last 5000 frames from the equilibrated trajectory with a single trajectory approach. Nmode calculations are particularly computationally demanding for large systems, such as protein-RNA, and the analysis was conducted by analyzing using 8 to 15 snapshots. For MM-GBSA, a level-set-based dielectric model was used with an ionic strength equal 150 mM, and a solvent probe of 1.4 Å radii. Using the same frames as those used in the binding free energy calculations; the interaction components for PAZ and siRNA complex were also determined for evaluating the individual binding free energy per-residue at complex interface. For multiple short simulations, it was considered the average values from the 6 independent simulations, uncorrelated data points. SEM values were simply computed, dividing the sample standard deviation by the square root of the total number of simulations.