Fast gaze reorientations by combined movements of the eye, head, trunk and lower extremities

Large reorientations of the line of sight, involving combined rotations of the eyes, head, trunk and lower extremities, are executed either as fast single-step or as slow multiple-step gaze transfers. In order to obtain more insight into the mechanisms of gaze and multisegmental movement control, we have investigated time-optimal gaze shifts (i.e. with the instruction to move as fast as possible) during voluntary whole-body rotations to remembered targets up to 180° eccentricity performed by standing healthy humans in darkness. Fast, accurate, single-step movement patterns occurred in approximately 70 % of trials, i.e. considerably more frequently than in previous studies with the instruction to turn at freely chosen speed (30 %). Head-in-space velocity in these cases was significantly higher than during multiple-step transfers and displayed a conspicuously regular bell-shaped profile, increasing smoothly to a peak and then decreasing slowly until realignment with the target. Head-in-space acceleration was on average not different during reorientations to the different target eccentricities. In contrast, head-in-space velocity increased with target eccentricity due to the longer duration of the acceleration phase implemented during trials to more distant targets. Eye saccade amplitude approached the eye-in-orbit mechanical limit and was unrelated to eye/head velocity, duration or target eccentricity. Overall, the combined movement was stereotyped such that the first two principal components accounted for data variance almost up to gaze shift end, suggesting that the three mechanical degrees of freedom under consideration (eye-in-orbit, head-on-trunk and trunk-in-space) are on average reduced to two kinematic degrees of freedom (i.e. eye, head-in-space). Synchronous EMG activity in the anterior tibial and gastrocnemius muscles preceded the onset of eye rotation. Since the magnitude and timing of peak head-in-space velocity were scaled with target eccentricity and because head-on-trunk and trunk-in-space displacements were on average linearly correlated, we propose a separate controller for head-in-space movement, whereas the movement of the eye-in-space may be, in contrast, governed by global, i.e. gaze feedback. The rapid progression of the line of sight can be sustained, and the reactivation of the vestibulo-ocular reflex would be postponed, until gaze error approaches zero only in association with a strong head-in-space neural control signal.

Introduction
Upright human subjects reorient to angularly distant targets (≥90°) by assembling fairly stereotypical, simultaneous, combined movements of the eye, head, trunk and lower extremities (Land 2004;Anastasopoulos et al. 2009). Such large reorientations to remembered targets can be accomplished as accurate, single-step shifts of the visual axis, covering at least 85 % of target eccentricity. Frequently, however, one or more short intervals of stationary gaze (i.e. of zero gaze velocity) emerge before target acquisition. In these cases, the initial segment of the gaze shift falls short of the target and more than 50 % of the target eccentricity is covered by a series of small saccades, presumably nystagmic fast phases, and head-in-space displacement (multiple-step gaze shifts; cat: Bergeron and Guitton 2002;humans: Anastasopoulos et al. 2009).
During the less time-efficient multiple-step gaze transfers, target acquisition takes at least 100-300 ms longer. It has been suggested that the higher the head-in-space velocity during this task, the more likely that a single-step gaze shift pattern will be elicited (Anastasopoulos et al. 2009). For example, patients with Parkinson's disease and spasmodic torticollis generate accurate single-step gaze displacements to large target eccentricities less frequently than normal subjects, most likely because of the low head-inspace velocity which they can attain (Anastasopoulos et al. 2011, but this improves after deep brain stimulation treatment (Lohnes and Earhart 2012). However, the question of whether increasing head (and other body segments) velocity during whole-body gaze transfers will result in a higher proportion of single-step gaze shifts has never been directly investigated. In simple ecological terms is this single-step gaze transfer pattern the one we use when we want to look at an eccentric target fast?
Not only is the effect of an imposed time demand on the task not known; the exact kinematic or neural constraints which have to be fulfilled so that the more time-efficient single-step gaze shifts can be assembled are not known either. The initial segment (up to 90°) of such gaze reorientations is characterized by tight segmental coupling, as the three mechanical degrees of freedom available (i.e. eye-inorbit, head-on-trunk and trunk-in-space) are reduced to two kinematic degrees of freedom (Sklavos et al. 2010). It is not, however, clear how this reduction is related to specific metric/kinematic measures of the coordinated movement.
A limitation of our previous studies is that data were obtained from relatively older subjects (Ss), turning at natural, freely chosen speeds (Anastasopoulos et al. 2009;Sklavos et al. 2010). Also, these studies assessed behaviour in a relatively small number of trials and only during the initial, brief segment of gaze progression. Thus, we do not know whether the reduction in the degrees of freedom prevails until target acquisition, given the increased contribution of the trunk in later stages of gaze shift (Sklavos et al. 2010). This is important because the trunk is the heaviest body segment and biomechanical challenges in controlling its motion accurately might interfere with the tight segmental coupling observed in the initial motion stages. In the current paper, we forced young subjects to execute the task as soon as possible and thus increase head-in-space velocity. In doing so, we expected to induce single-step gaze transfer patterns more frequently, and this should allow us to determine the statistical significance and degree of covariation among eye, head and trunk movements in late stages of gaze transfer, that is, almost until arriving on the target. Also, it can be expected that the dynamic trajectories of the involved body segments in such time-optimal task executions might reveal the time course and switching of the neural control signals of the underlying control mechanism; because of the instruction and the considerable head and trunk inertia, the magnitude of neural control signals to the effectors moving these segments is expected to be either maximal or minimal under these circumstances and switching between acceleration and deceleration sharply defined ('bang-bang control'). Finally, the timing of segmental activation within this motor synergy was further evaluated by recording the initiation of electromyographic (EMG) activity from several muscles involved in stepping. A top-down sequence of segmental movement initiation has been observed when subjects turn on the spot in upright posture with the eye starting to rotate first and the foot last (Hollands et al. 2004), but data show that eye saccade initiation is considerably delayed when compared with data from sitting subjects (Scotto Di Cesare et al. 2013). Therefore, we hypothesized that this saccadic delay is introduced in order to optimize the coordination between the large number of rotating segments and the anticipatory postural adjustments required before the onset of movement. EMG recordings may disclose the timing and pattern of activation of several lower limb muscles in relation to eye rotation onset and spatial prediction.

Subjects
We studied six healthy right-hand/leg dominant adults (three females), mean age 27 ± 1.2 years. No one wore spectacles, and selection ensured excellent physical condition. The study was approved by the Imperial College Riverside Ethics Committee in accordance with the Declaration of Helsinki, and all subjects provided informed consent.

Protocol and apparatus
Details have been given previously (Anastasopoulos et al. 2009). In short, Ss stood in the centre of a circular array (radius 1.2 m) of eight LEDs, placed at 45° intervals at eye level, in darkness. At the beginning of each trial, they fixated and aligned their head, body and feet with the central LED. After a delay of 10 s, the central target was extinguished, thus indicating that another LED in one of the seven eccentric locations (45°, 90° and 135° either right or left of centre as well as at 180°) had been lit. The subject had to fixate, turn and align his body with the lit LED (outbound trials). After an interval of 15 s, the eccentric LED was turned off (Fig. 1, upper schematic), thus cueing subjects to return to the initial, central position (inbound or 'return' direction) and look in the direction of previously illuminated central LED. The central LED was, however, off during returning, and Ss moved in complete darkness (see below). Ss performed eighteen trials (outbound followed by return) to each LED location in random order.
It should be noted that trials to 45° targets are visually driven. For outbound trials to 90° or larger targets, the LED is initially not visible and Ss had no hint as to its location or in which direction they should turn. In inbound or 'return' trials, however, Ss knew where the centrally positioned target was and in this way the combined movement to the centre was exclusively memory driven (the central LED was turned on after additional 2.8 s).
In contrast to all our previous studies with this paradigm, where Ss were given no instruction as to the speed of the movement required, here Ss were instructed to accomplish turns as soon and as fast as possible.
Head-in-space, upper trunk and feet horizontal (yaw plane) movements were recorded using a Polhemus Fastrak Fig. 1 A rightward singlestep gaze shift to the central remembered target at 135° offset in complete darkness. The panels above show a cartoon with the successive target presentations and head/trunk positions adopted before and after the combined movement. The left and middle vertical dotted markers indicate the onset and termination of the rapid displacement of the line of sight. Above and below the EMG recordings are shown position and velocity traces, respectively. Note the reacceleration in the eye velocity trace after the offset of the main eye saccade. Both during the acceleration and deceleration phases of head-in-space movement, changes in head-on-trunk velocity (second trace from the bottom) are complemented by variations of trunk-in-space velocity, such that a regular, approximately bell-shaped head-in-space velocity profile is maintained throughout the movement. Reappearance of the central LED is indicated by the right vertical marker. Synchronous EMG activity in the anterior tibial (tibant) and gastrocnemius (gastr) muscles precedes eye, head and trunk rotations . Markers were placed on a light, adjustable helmet, on the spinous process of the seventh cervical vertebra and dorsally on each foot. Horizontal eye-in-orbit rotations were recorded using bi-temporal DC electrooculography (EOG) with a flat response to 90 Hz. EMG activity was collected by adhesive surface disc electrodes (CareFusion, Middleton, WI, USA, diameter 1 cm) placed at a distance of 3 cm near the motor point of the following muscles on the right side: anterior tibial, gastrocnemius, biceps femoris and external oblique.

Data acquisition and analysis
The on-off LED signals, EOG, body position markers and EMG were sampled at 500 Hz and stored for offline analysis. A smoothing Savitzky-Golay filter was used for offline velocity measurements of the EOG position trace (MAT-LAB™ function sgolayfilt), while head and trunk velocities were calculated after subjecting the position signals to moving mean linear filtering (16 + 1, not weighted). Movement onset/offset for the various segments was calculated from the time derivative. Gaze (eye-in-space) displacement was obtained by adding EOG and head signals after appropriate scaling. Subtracting the trunk from the head signal resulted in an estimate of head-on-trunk movement. Biphasic deflections in the eye velocity trace were thought to represent blinks (Gandhi 2012); these trials were not further analysed. Raw EMG data were low-passed filtered at 10 Hz. Onset of EMG activity in relation to the onset of gaze shift (i.e. of the eye saccade) was defined visually (Fig. 1). To improve consistency, only rightward trials to initially non-visible targets (i.e. target eccentricity ≥90°, outbound and inbound) starting with the right foot were analysed for EMG onset.

'Within channels' principal component (PC) analysis
We applied PC analysis to examine and quantify statistically the relationships among a set of 3-covarying variables (channels), i.e. eye-in-orbit, trunk-in-space and head-ontrunk, of n-observations by transforming the original data to a new set of 3-orthogonal variables, i.e. the PCs. The three PCs are derived in a decreasing order of importance (statistical significance), so that the first and the last PCs, respectively, account for maximal and minimal amount of variance in the original data. The PCs combine linearly with the original (elemental) variables, and if the three elemental variables are totally uncovaried, three statistically significant PCs will account for data set variability. In contrast, if two of the variables covary, then only two significant PCs account for data set variability and the space defined by the three elemental variables is reduced by one dimension, i.e. variation in the data set under consideration can be described by a plane (i.e. a two-dimensional coordinate system). Similarly, if all three variables covary, then only one significant PC accounts for data set variability and the space defined by the three variables is reduced by one further dimension (i.e. to a straight line). Thus, the number of significant PCs and the percentage of data variance they account for stand for the degree of covariation among the variables. Percentages of the variation in the original set of variables, which the principal components account for, equal the normalized eigenvalues ˆ j = j / 3 j=1 j ,ˆ 1 <ˆ 2 <ˆ 3 of the correlation matrix Σ = Tr(X) · X (Tr is for transpose) where is the matrix of the data set with n rows (samples) and 3 columns (variables). E (m) represents eye, T (m) trunk and H (m) T head-on-trunk angular displacement recorded as the mth sample of the trial, m = 1, …, n. For the purpose of PC analysis, raw position data were used (i.e. before filtering), and EOG was resampled at 30 Hz. Increasing time intervals, starting always with the onset of gaze, was analysed ('sliding' approach). Data centring and auto-scaling were carried out by subtracting the arithmetic means µ(E (m) ), µ(T (m) ), µ(H (m) T ) and dividing by the standard deviations σ (E (m) ), σ (T (m) ), σ (H (m) T ), respectively. By multiplying the array of the data set consisting of the matrix X n×3 with the array of eigenvectors (loadings, v 3×3 ) of Σ, we getP = X · v, where P n×3 is the array of the scores of the PCs. These are organized in columns of an increasing order of importance, such that the third, second and first columns contain the first, second and third principal components, respectively, accounting for ˆ 3 >ˆ 2 >ˆ 1 of the variance of the original variables, respectively. Then, the contribution of the X (i) original variable to the P (j) principal component (weight) is v ij . The loadings of the three variables define the new coordinate system. Eigenvalues ( j ), j = 1, 2, 3 and their corresponding eigenvectors (v j ), j = 1, 2, 3 were computed by using the MATLAB™ function eig. For each data set (trial), the eigenvalue (k) j and weights v (k) ij , k = 1, …, N, are calculated. Further, by averaging across all trials, the mean eigenvalue

PCs' statistics
Bartlett's test of sphericity was used to validate the hypothesis that the first principal component was significantly different from all those of lower importance, i.e.
both second and third. This has to be accepted if the test statis- is greater than the critical value of Chi-square distribution with 1 2 (j − 2)(j + 1) degrees of freedom. As Bartlett's test is insufficient for testing the significance of the last two principal components, i.e. those of the lowest importance (j = 1, 2), a test for proportion of total variance was applied post hoc. A desired proportion of total variance at 99.0 % was selected and excluded all principal components, which accounted for <1 % of the total variance; one-sample independent t test was performed for the arithmetic means of variance proportion of second and third PCs and tested whether they were significantly >1 %. When normality failed according to a one-sample Kolmogorov-Smirnov nonparametric test, a binomial nonparametric test for the median value was carried out instead.

Multilinear regression
) displacement trajectories across all N-trials, separately for various target eccentricities and time intervals from gaze onset (Table 1), were calculated. Any point of trunk displacement trajectory T (a) where ε i are fitting standard errors. This can be written as a matrix equation by introducing vectors The estimate of b which minimizes the sum square error Tr(ε) · ε, as in multiple linear regression, is then given by b = (Tr(A) · A) −1 · (Tr(A) · T (a) ). The new estimate of trunk displacement trajectory is T (a) = A ·b and, in principle, the 3D vector (E (a) , H (a) is an estimate of the goodness of fit of the surface plane to the vector of displacement trajectories. Finally, are the angles of orientation (in radians) of the surface plane with respect to trunk-in-space/head-on-trunk and head-on-trunk/eye-inorbit planes, respectively.

Patterns of voluntary reorientations
Details of the movement pattern in normal subjects moving at freely chosen speed have been given and illustrated extensively in earlier publications (Figs. 1, 2, 3, Table 1 Covariation of eye-in-orbit, head-on trunk and trunk-in-space displacements Segmental coupling is quantified by the percentage (%) of variance accounted for by principal components. The first and second PCs were overall significant (p < 0.001). The third PC accounted for <1 % of displacement variance and was non-significant (p > 0.1). T time interval analysed during single-step gaze transfers separately to 90°, 135° and 180° targets, starting always from gaze onset; N, number of trials analysed; G, mean amplitude ± SD in degrees, of gaze displacement during the time interval under consideration. The orientation of the covariation plane is indicated by the angles ϕ (in degrees, with respect to the plane defined by the axes eye-in-orbit/head-on-trunk) and θ (with respect to the plane defined by the axes trunk-in-space/head-on-trunk). The coefficient of determination R 2 was always >96 % In the outbound condition, particularly of ≥90° amplitude and, hence, of unpredictable target location, the prevailing movement pattern involved termination of the primary gaze shift well before target acquisition. Compensatory, slow-phase eye movements would rotate the eye centrally, and one or more intervals of stationary gaze, due to the vestibulo-ocular reflex, would appear (multiplestep gaze shifts as defined in our previous publications).
Outbound turns and visually guided gaze shifts (i.e. to ±45° visual targets) will not be considered further here.
Inbound or 'return' reorientations to targets ≥90° consisted predominantly of single-step displacements of the visual axis covering at least 85 % of target eccentricity (Fig. 1). Here, we concentrate mainly on the description of these responses to the centrally located, remembered target and try to identify the kinematic constraints allowing for the tight segmental coupling and the execution of such quick and accurate gaze saccades. Single-step movement patterns occurred in 82, 78 and 56 % of inbound trials to 90°, 135° and 180° target eccentricities, respectively, i.e. considerably more frequently than in our previous studies with the instruction to turn at freely chosen speed (60, 35 and 25 %, respectively).
Metrics and kinematics of time-optimal, single-step gaze transfers As already reported in the literature for large saccades with the head free to move (reviewed by Freedman 2008), peak eye saccade velocity (448 ± 135°/s, mean ± SD) 1 s

Head-in-space velocity
Trunk-in-space velocity Head-on-trunk velocity Fig. 2 A family of eight superimposed head-in-space velocity trajectories from single-step gaze shifts. Displacements to 90° (left) and 135° (right) eccentricity targets (upper rows) and the corresponding trajectories of trunk-in-space (middle rows) and head-on-trunk velocity (lower rows) from the same trials. All of the traces are aligned on the peak of head-in-space velocity and normalized such that headin-space velocity has in all cases the same arbitrary value. They are ranked according to the displacement magnitude of peak head-ontrunk deflection. Note how the cooperative or opposing motion of head-on-trunk and trunk-in-space produces in all cases a regular, approximately bell-shaped velocity profile of head-in-space movement Kinematics of head-in-space movement. Though quite variable, peak head-in-space velocity and acceleration time interval increase significantly with target eccentricity and duration (261 ± 87 ms) were unrelated to eye saccade amplitude (46.8 ± 5.6° on average, approaching thus the mechanical eye-in-orbit limit). An example executed under the current 'as fast as possible' instruction is shown in Fig. 1. In most instances, the eye saccade velocity profiles had a single peak. In several occasions, however, a second, considerably smaller, peak was discernible shortly before eye saccade offset. In other instances, one or two tiny eye saccades (of <4° amplitude) appeared before gaze shift end (Fig. 1). These reaccelerations occurred almost always after the velocity of head-in-space motion had begun to decrease but, still, before gaze end (see below).
Gaze transfer was terminated by the emergence of VOR (at mean amplitude 84.6 ± 5.0°, 131.6 ± 8.1° and 172.1 ± 9.8°; unsigned errors 6.4°, 7.1° and 9.4° for 90°, 135° and 180° target eccentricity, respectively), and the line of sight was stabilized, in darkness, quite near the memorized centrally located target. Ss sometimes made a corrective small saccade after the LED was turned on.
Head-in-space movement resulted from the sum of headon-trunk and trunk rotation in space. The dynamic trajectory of head-in-space velocity was conspicuously monophasic, increasing smoothly to a peak and then decreasing slowly until realignment with the memorized gaze direction (deceleration phase usually somewhat longer, Figs. 1 and 2). Head-in-space peak velocity (Fig. 3a) and duration increased thereby as a function of target eccentricity (factorial ANOVA, F = 26.1 p < 0.0001 and F = 25.1 p < 0.0001, respectively). In contrast to head-in-space, the dynamic trajectories of trunk-in-space, and most particularly headon-trunk motion, were more variable (Fig. 2). Just as peak head-in-space velocity and duration, the head-in-space acceleration phase was also significantly dependent on target eccentricity, i.e. it was shorter when moving to 90° targets and longer when reorienting to 180° targets (ANOVA, F = 18.0, p < 0.0001, Fig. 3b; compare left and right records Fig. 2). Head-in-space acceleration during this time interval was not different among the various target eccentricities (on average 1362 ± 742°/s 2 ), confirming that the Ss generally tried to comply with the instruction 'go as fast as you can'. Noteworthy, in reorientations made by means of the multiple-step pattern (Fig. 4), mean head-in-space acceleration was significantly low (on average 593 ± 310°/ s 2 , i.e. less than half the acceleration observed when assembling single-step gaze shifts, ANOVA, F = 49.4, p < 0.0001). In many instances, the premature termination of primary gaze shift long before reaching the position of the centrally located remembered target was signalled by an abrupt deceleration break in the, up to this point, monotonically increasing head-in-space velocity trace (arrows in Fig. 4). As a result of these modifications, peak velocity was significantly lower during multiple-step than during single-step gaze shifts (172 ± 52 vs. 279 ± 80, 211 ± 64 vs. 363 ± 110 and 266 ± 54 vs. 428 ± 127°/s to targets of 90°, 135° and 180° eccentricity, respectively, F = 29.7, p < 0.0001). However, even if interrupted by variable periods of stationary gaze, the displacement of the visual axis was terminated again very close to the remembered target.
As already reported in the literature (reviewed by Leigh and Zee 2006) for smaller head-free gaze displacements, gaze saccade duration, though quite variable, correlated significantly and increased linearly with gaze amplitude in single-step gaze transfers (n = 181, R = 0.7, slope 3.2 ms/°). Single-step transfers were approximately 100-300 ms faster than multiple-step gaze displacements (315 ± 75, 449 ± 107, 599 ± 135 vs. 506 ± 78, 611 ± 132, 791 ± 138 ms to targets of 90°, 135° and 180°  Fig. 4 A rightward multiple-step gaze shift to the central remembered target at 135° offset in complete darkness. Target presentations and head/trunk positions adopted before and after the combined movement are as shown in the cartoon of Fig. 1. Extinguishment of the eccentric LED ('go' signal) and reappearance of the central LED are indicated by the left and right vertical markers, respectively. The premature termination of the initial gaze shift long before target acquisition (arrowhead, bottom trace) coincides with the emergence of a plateau of the up to this point increasing head-in-space velocity trace (arrow) eccentricity, respectively, F = 82.8, p < 0.0001). Finally, head-in-space oscillations (quantified by the absolute ratio of the first negative peak to peak head velocity, on average 0.06 ± 0.03) did not increase with the amplitude of headon-trunk movement (R 2 = 0.24).
In order to examine whether head-in-space control is autonomous or driven by an ultimate goal-oriented function (in this case gaze reorientation), we examined how the timing of peak head-in-space velocity was associated with the end of the gaze shift in single-step movement patterns (Fig. 5). If head-in-space movement was solely controlled by gaze signals, it might be expected that, because of the forced instruction, head-in-space would continue to accelerate until arriving at the desired gaze position (final goal, vertical middle marker in Fig. 1) and then start to decelerate, as gaze shift end does not coincide with head-in-pace alignment with the remembered target. Note that in Fig. 5, however, peak head-in-space velocity does not coincide with gaze shift end; it occurs earlier and earlier as gaze shift duration and amplitude increase, suggesting that the control of head-in-space motion does not exclusively aim at conveying the line of sight as soon as possible to the target. In contrast, the duration of the 'acceleration phase time' increases with target eccentricity, implying that the latter is one of the factors dictating the timing of peak head-inspace velocity.
Covariance of eye, head and trunk motion during single-step gaze transfers (principal component analysis, PCA) Mean ± standard deviations of percentages of variance accounted for by individual principal components (PCs) are shown in Table 1, separately for single-step gaze shifts to 90°, 135° and 180° targets. PCA was applied whenever sufficient numbers (N > 15) of trials for statistical analysis were available. The first PC was significant (Bartlett's test, p < 0.001) and accounted for almost 90 % of variance in all analysis time windows. The second PC was also significant overall and accounted approximately for about 10 % of variance (p < 0.001). In contrast, the third PC was non-significant in all time windows and target eccentricities, accounting for <1 % of variance. Variance up to half a second after movement initiation could be thus accounted for by essentially the first two PCs (i.e. the sum of their percentages was approximately 100 %), suggesting that the three mechanical degrees of freedom under consideration are reduced to two kinematic degrees of freedom. Visualization of these relationships is exemplified in Fig. 6. Mean eye, head and trunk displacements from 53 trials to 135° targets up to 333 ms (i.e. 10 data samples) are plotted in the three-dimensional position space defined by the variables (eye-in-orbit; head-on-trunk and trunk-in-space). Time progresses from the origin of the coordinate system to the right, upward and backward (Fig. 6a). Two significant PCs imply that the data approximately vary in only two orthogonal directions defined by the first and second PCs, i.e. in a plane. The projection of data samples in this plane is illustrated in Fig. 6b (x and y axes give mean PC scores, units are normalized as a consequence of autoscaling, c.f. "Materials and methods"). PC1 increases with time (from −1.7 to 2.6), while PC2 first increases (from −0.4 to 0.5) and then decreases to −0.5. All variables have positive mean PC1 weights (i.e. show the trend represented by PC1). Eye mean PC2 weight is again positive, while those of head and trunk are both negative. In consequence, head and trunk displacements covary in both directions described by the significant PCs.
An approximation of the orientation of the covariation plane, calculated from displacement data by means of multilinear regression, is depicted in Fig. 6a and described by the angles θ and ϕ. Note from Table 1 that θ is often close to 90° while ϕ is mostly larger than 45°. Because θ is mostly close to 90°, the covariation plane is approximately perpendicular (orthogonal) to the plane defined by the axes head-on-trunk and trunk-in-space. In consequence, headon-trunk and trunk-in-space displacements are linearly correlated, i.e. they lay approximately upon a regression line and can be replaced by one variable, for example, headin-space. The slope of this relationship is larger than unity Head-in-space attains peak velocity well before gaze shift end (vertical marker). The time span between these events increases with gaze amplitude. This is exemplified for 90°, 135° and 180° reorientations (lower, middle and upper traces, respectively) on the right. The shaded area beneath the head-in-space velocity traces represents the contribution of head-in-space displacement to gaze amplitude as ϕ is mostly >than 45°. In this way, trunk-in-space is on average larger than head-on-trunk displacement throughout the gaze transfers. Of note, these relationships do not apply during 'multistep' (nystagmic) transfers. Here, three derived significant PCs account for the data variance (Sklavos et al. 2010).

Onset of EMG activity in lower limb and trunk muscles
Ss activated the right anterior tibial and gastrocnemius muscles approximately 30 ms on average before starting to rotate gaze (i.e. the eye) rightwards to the centrally located remembered target (Fig. 7, black quadrangles and circles; also Fig. 1). Biceps femoris and lateral oblique (black triangles) were activated approximately 40-50 ms later than eye movement onset. Note the larger scatter of values of onset of EMG activity in these proximal muscles. When the location of the target was not known beforehand, i.e. in the unpredictable outbound trials, both distal and proximal muscles were activated mostly after the initiation of eye rotation (grey quadrangles, circles and triangles, p < 0.003 Mann-Whitney U test, inbound vs. outbound). There was not any statistical difference in the onset of EMG activation from the right tibial and gastrocnemius muscles, both for inbound and for outbound reorientations, i.e. they were activated on average synchronously. Not unexpectedly, EMG activity in these muscles decreased as the right foot started to rotate rightwards.

Discussion
More than 70 % of reorientations to the centrally located remembered target in this study were accomplished by rapid, accurate single-step movement patterns (in contrast to our previous studies at free speed where single-step gaze shifts were observed in approximately 30 % of cases; Anastasopoulos et al. 2009). Ss were instructed to move as fast as possible and were thus able to attain high head and trunk velocities. In the remaining 30 % of trials, however, the progression of the line of sight was interrupted by one or more variable intervals of stable gaze (multiple-step, nystagmic or 'scanning' pattern). All the same, gaze displacement ended again close to the remembered target. Multiple-step gaze shifts were characterized by significantly low head-in-space velocity. High head-in-space velocity during task execution thus increases the probability of eliciting  . 7 Latency of EMG activity relative to gaze shift onset. Median and interquartile ranges in the right tibialis anterior (tibant), gastrocnemius (gastr), biceps femoris (bifem) and obliquus externus (oblex) muscles. Gaze onset is indicated by the horizontal dotted line. Only those rightward rotations are included during which stepping started with the right foot. While median EMG onset precedes or coincides with eye rotation when the location of the target is known (inbound or return trials, black quadrangle, circle and triangles), it follows the rotation onset in outbound, unpredictable trials (same notation in grey) single-step gaze shifts and quickly acquiring a predetermined direction of the visual axis, as hypothesized in previous studies (Anastasopoulos et al. 2009(Anastasopoulos et al. , 2011. The slower, multiple-step gaze shift pattern, incorporating sections of VOR mediated gaze stability, looks suitable for visually scanning the environment during lower velocity roaming head movements. Subjects confronted with the task to acquire the target as soon as possible quickly increased the velocity of the head-in-space to a peak before decelerating up to the target. Statistical evaluation of segmental spatial covariation by PCA confirmed that eye, head and trunk displacements are coupled at the kinematic level such that their motion can be controlled until gaze shift end with a two degrees of freedom system (i.e. eye-and head-in-space). This presumably presents the CNS with a simpler control task.
Separate control of head-in-space movement A conspicuous, central feature of large, single-step reorientations to remembered targets was the high peaked monophasic head-in-space velocity traces (Fig. 2), suggesting that head-in-space rotation is related to control signals switched from maximum 'on' to maximum 'off' ('bang-bang' double pulse). This was a consequence of the forced instruction to the Ss and indicates that inertial characteristics of the head/trunk mechanical system predominate. Head-on-trunk and trunk-in-space movements were, instead, more variable both during the accelerating and decelerating phase ( Fig. 2; also before and after gaze shift end). Further evidence suggestive of a separate or independent head motion controller is the fact that time to peak velocity and the magnitude of peak velocity for the 'composite' head-in-space variable increased with target eccentricity (Fig. 3). Similar 'main sequence' relationships have been earlier described for trunk restrained primates (Tomlinson and Bahra 1986) or cats (Guitton et al. 1990), i.e. when head-on-trunk equals head-in-space movement. Zangemeister et al. (1982) studied time-optimal head movements in sitting humans and found that the duration of head motion was independent of movement amplitude. In contrast to the results of Zangemeister et al., head movement in our study appears to be controlled by the duration of force generation rather than by the amount of force (Fig. 3), probably because of the large target eccentricities and the large trunk inertia involved. The head control system is, in the current experiments, clearly functioning close to the upper limit of force generation, and this strict demand may be responsible for the fact that Ss fail to acquire in all occasions distant targets with single-step gaze shifts (see below). Under the requirement of maximal effort in executing a heavy task, it is conceivable that Ss had to prolong the application of force and consequently failed more frequently with increasing target eccentricity.
A further salient finding of our study in favour of separate control for head-in-space is that command timing does not relate to gaze displacement but, instead, head-in-space displacement (Fig. 5). If the head was exclusively under gaze control, a continual increase in head-in-space velocity up until gaze shift end would be expected. Our findings contribute to the ongoing debate whether it is gaze or its individual eye and head components that are controlled during orienting behaviours (Guitton et al. 2003;Choi and Guitton 2006;Freedman 2001;Sylvestre and Cullen 2006) with evidence supporting the view that the controller for gaze and head-in-space movement may not be identical. The neural mechanism underlying the control of the large number of muscles and joints of the neck and trunk, such that a solution appropriate for the task in hand (i.e. timeoptimal displacement of head-in-space) can be achieved, remains to be elucidated. Notably, single-step displacements are characterized on average by linear correlation of head-on-trunk and trunk-in-space displacements throughout the gaze transfer. Moreover, trunk displacement is on average larger than head-on-trunk displacement throughout visual axis transfer. Multiple-step movements do not fulfil these kinematic constraints as they do not demonstrate reductions in significant PCs and covariation among the segments (Sklavos et al. 2010). The present results may facilitate future investigations on the control mechanism. They have to be interpreted, however, with caution as Ss were required to reestablish the same body position, and this may have encouraged coordinated, en-bloc movements, reducing variability. The study of more complex or ecological tasks in the future may refine our results and reveal further parameters of segment covariation.
Why did Ss frequently choose to move head-on-trunk far below 60-70° (i.e. the limit of the so-called neutral zone, McClure et al. 1998) though they thus increase the duration of head-in-space movement up to the target (Fig. 2)? Our results show that this strategy does not reduce unwanted head-on-trunk oscillations (Saglam et al. 2011). It possibly represents an anticipatory mechanism aiming at reducing the amount of control energy required to rotate headon-trunk to more eccentric positions due to concomitant increase in the applied force (Kardamakis and Moschovakis 2009). Contrary to these feed-forward models, we suggest, a substantial contribution from feedback control, based mainly on vestibular input, because the movement of the support (i.e. the feet) is extremely variable and subject to many possible disturbances (Sklavos et al. 2010). It may also be easier to coordinate the many segments and achieve a well scaled and timed head in space velocity, as our Ss did managed to do, if afferent online velocity information (=vestibular input) is available. This question could be answered by testing patients with bilateral vestibular failure during whole-body gaze transfer task, and preliminary data in two such patients support this view (ongoing project).
The control of eye movement and constraints on gaze shift duration Separate control signals for the eye and head have been proposed by several researchers because of the variability in the relative timing of the start of eye and head movements and the effect of initial eye-in-orbit position on the contribution of the eyes and head to gaze displacement (for a review, see Freedman 2008). Separate control for the most distal and proximal parts of linked systems has been modelled recently by Daye et al. (2014). According to this scheme, the most distal subpart (here the eye) of such a 'hierarchical' structure is governed only by global goal feedback (here gaze). Indeed, the relationships between amplitude, velocity and duration of the eye-in-orbit movement of large single-step gaze displacements were not stereotyped but depended on how the head moved. For example, reaccelerations of the eye-in-orbit motion after the peak of head-in-space velocity may have been generated to correct deviations of the progressing gaze trajectory from the efference copy of ongoing motor commands and complement the separately progressing head motion before the impending gaze shift termination (Fig. 1). Remarkable accuracy in darkness (>90 %) was achieved in spite of the variable head-in-space contribution even for multiple-step displacements. We suggest that a mechanism monitors continually overall visual axis progression and modifies accordingly the eye component, given that head-in-space movement has to obey bangbang control appropriately for a certain target eccentricity.
The fact that high head-in-space velocity and single-step gaze shifts are interconnected makes sense biologically when target acquisition occurs under 'as soon as you can' conditions, but how are these two variables associated mechanistically? It is clear from our findings that reemergence of the VOR and gaze stabilization does not depend solely on gaze error, as most gaze models propose (cf. Galiana and Guitton 1992), but they are related to the strength of the descending signal moving head-in-space, i.e. the command sent to agonists moving 'head-in-space'. In multiple-step gaze shifts, characterized by premature termination of 'saccadic mode' in spite of the instruction, average head-in-space acceleration was less than half of that during single-step reorientations (cf. "Results", compare Fig. 1 with Fig. 4). We argue that high head-in-space acceleration results both in adequately fast head-in-space movement, contributing quickly to approximately 45°, 90° and 135° of displacement for targets at 90°, 135° and 180°, respectively (grey area under the velocity traces, given that eye rotation covers utmost 45-50° on average, Fig. 5) and in prolongation of the time interval available for 'saccadic mode'. The neural mechanism that determines the duration of 'saccadic mode' and the termination of gaze or eye saccades is largely unknown. It has been hypothesized that descending signals from the rostral superior colliculus and/or frontal eye fields (Bergeron and Guitton 2002;Hanes et al. 1998;Segraves 1992) gate omnipause neuron activity (which resumes at saccade end) and by this disinhibit gaze bursters in the pons allowing for gaze displacement. Electrical stimulation of the omnipause neurons brainstem region in primates stops ongoing gaze shifts in mid-flight (Gandhi and Sparks 2007). During the stimulation, interval gaze position remains constant ('stabilization mode') and, while head movement continues along its trajectory, the eye counter-rotates with a roughly equal velocity in the opposite direction (i.e. the VOR is reactivated). Still, mechanisms other than OPN discharge may contribute to saccadic termination. The oculomotor cerebellum has been also implicated in the termination of rapid gaze displacements (Fuchs et al. 2010;Rucker et al. 2011). We suggest a straightforward, mechanistic explanation of our findings, i.e. that the gating period of this putative mechanism is too short if the descending signals accelerating head-in-space are weak. The stronger demands imposed on the duration of head-in-space movement due to the large eccentricities and large segmental inertia revealed the conditional ability of the gaze control system to delay beyond a top limit switching from 'saccadic mode' to 'stabilization mode'. The mechanism is thus not uniquely controlled by a gaze-motor-error signal (cf. Paré and Guitton 1998). Still, the experimentally observed relationship between the strength of the 'head-in-space' neural control signal and the duration of the 'saccadic mode' interval may be not a causal one; the reverse may be also true, i.e. head-inspace velocity is simply low when multi-step displacements are planned, e.g. when a visual scanning pattern is called for.
Stabilization of the ankle joint precedes eye rotation to the target Foot rotations in this experimental paradigm are extremely variable and start mostly considerably later than movements of the eye, head and trunk (on average 350 ms later than the eye, Anastasopoulos et al. 2009). As foot latency is significantly reduced when target location is known (i.e. in return trials), it has been previously hypothesized that proprioceptive input and local mechanisms (i.e. from hip stretch) are overruled by supraspinal inputs and that stepping initiation is thereby controlled more rostrally. On the other hand, segmental movement variability analysis has previously indicated that that lower extremities contribute initially to postural rather than gaze displacement control (Sklavos et al. 2008). Indeed, our current results suggest that the turning synergy to predictable targets is initiated with the synchronous activation of antagonistic leg muscles. Stiffness of the ankle joint thus increases, and we hypothesize that forces involved in trunk rotation will less likely threaten balance. More detailed measurements are, however, required in this direction. Also, EMG activity in muscles that rotate the foot may be expected to ensue, but this was not recorded and cannot be answered with certainty. Such ankle stabilization is likely part of an anticipatory postural adjustment in preparation of a complex movement sequence (Assaiante et al. 2000), in this case ultimately aiming to shift the eyes in space. Given that in conventional head-fixed saccadic experiments, anticipatory adjustments and assemblage of multisegmental synergies are not required, our EMG findings thus agree with the view that the long eye latencies observed in this paradigm relate to active delaying of the eye saccades until all anticipatory postural adjustments and segments are ready to go (Scotto Di Cesare et al. 2013).


extfieldsofstudy: Biology Medicine

text:
Molecular Modeling of Myrosinase from Brassica oleracea: A Structural Investigation of Sinigrin Interaction

Myrosinase, which is present in cruciferous plant species, plays an important role in the hydrolysis of glycosides such as glucosinolates and is involved in plant defense. Brassicaceae myrosinases are diverse although they share common ancestry, and structural knowledge about myrosinases from cabbage (Brassica oleracea) was needed. To address this, we constructed a three-dimensional model structure of myrosinase based on Sinapis alba structures using Iterative Threading ASSEmbly Refinement server (I-TASSER) webserver, and refined model coordinates were evaluated with ProQ and Verify3D. The resulting model was predicted with β/α fold, ten conserved N-glycosylation sites, and three disulfide bridges. In addition, this model shared features with the known Sinapis alba myrosinase structure. To obtain a better understanding of myrosinase–sinigrin interaction, the refined model was docked using Autodock Vina with crucial key amino acids. The key nucleophile residues GLN207 and GLU427 were found to interact with sinigrin to form a hydrogen bond. Further, 20-ns molecular dynamics simulation was performed to examine myrosinase–sinigrin complex stability, revealing that residue GLU207 maintained its hydrogen bond stability throughout the entire simulation and structural orientation was similar to that of the docked state. This conceptual model should be useful for understanding the structural features of myrosinase and their binding orientation with sinigrin.

Introduction
Myrosinase (also termed thioglucoside glucohydrolase, sinigrinase, or sinigrase EC: 3.2.1.147), is an S-glycosidase involved in breakdown of thioglucosides such as glucosinolates (GSL), a large group of secondary metabolites especially abundant in cruciferous plants. Myrosinase is localized in specialized GSL-free idioblast cells termed myrosin cells [1]. GSLs, classified based on their structure as indolic or aliphatic, are biologically inactive in their native state; myrosinase catalyzes the hydrolysis of these thioglucosides into D-glucose and a thiol (sulfate) such as aglycone (thiohydroximate-O-sulfonate). The thiols produced are unstable intermediates that differ mainly based on side chain structures, myrosinase substrate, plant species, and other reaction conditions like pH or availability of ferrous ions and environmental factors [2][3][4]. Compounds spontaneously formed from aglycone in this process are biologically intact, active and sometimes toxic, and include isothiocyanates and nitrile, epithionitriles, thiocyanates, and oxazoldine-thione [5], which are involved in plant defense against herbivory, pathogen invasion and tissue damage [1,3]. Concentrations of GSLs vary among different parts of the plant [6,7]. External factors such as tissue damage and pathogen attack trigger the release of GSLs, which results in interaction with myrosinase and eventually degradation of the glycosides. In plants, myrosinase activities vary between species, cultivars, and parts of the plant examined. Myrosinase activity is highest in seeds and seedling stages as indolic GSLs are broken down to form indoleacetic acid (IAA), a phytohormone essential for developing plants [8]. In cruciferous species, the interaction between myrosinase and GSLs for defense creates a so-called binary mustard oil bomb. These myrosinase are activated by ascorbic acid and the three-dimensional structures of several binary and ternary complexes of the enzyme with different transition state analogues, with ascorbic acid and with inhibitors mimicking the substrate have been reported. In addition, mechanistic details already known include the position of a water molecule placed perfectly for activation by ascorbate and for nucleophilic attack on the covalent intermediate of the reaction [9]. Biochemical, genetic and phylogenetic data, [10] reported that the myrosinase-GSL system evolved from cyanogenic glucosides and O-ȕ-glucosidases. Products of GSLs possess toxic and pharmacological effects on a wide range of organisms [11,12]. GSLs induce Phase I and Phase II detoxifying enzymes during cancer infection and alter steroid hormone levels during their metabolism, in turn aiding in protecting the cell from oxidative damage [13]. Cruciferous species are known for the presence of GSLs, conferring their chemo-preventive effects [14]. Broccoli (Brassica oleracea var. italica) was reported to act against and provide protection from liver cancer in animals [15]. In humans, glycosides are involved in reducing DNA damage in lymphocytes [16,17]. There was significant anti-nutritional activity in animals when they were fed with GSLs [14]. GSLs from plants have become an interesting area of work due to their benefits to animal and human health such as anti-microbial [18], anti-nematode [19] insecticidal, weed-inhibiting [20], herbicidal [10,21], allelopathic [22], and anticarcinogenic [15,23] activities. To date, more than 130 GSLs have been identified and reported in plants, of which more than 30 are found in the Brassicaceae family [24]. Lenman et al [25] reported that myrosinase shows differential expression throughout development and also in various parts of the plant. Molecular breeders have used traditional breeding methods to increase the content of GSLs due to their beneficial activities. This study aims to improve understanding of myrosinase-sinigrin interaction at the molecular level. To date, one conceptual model was constructed from Brassica juncia and described their myrosinase structural features [26]. Homology modeling is a powerful method to predict three-dimensional structures from amino acid sequences [27]. This approach can be used to obtain structural knowledge of proteins, enzyme-substrate interactions, and protein-ligand interactions [28,29]. In this study, we have applied this technique to predict the structure of myrosinase from Brassica oleracea (B. oleracea, wild cabbage) and examine its interaction with sinigrin, one of the possible substrates.

Template Identification and Homology Modeling
The amino acid sequence of myrosinase (MYR) from cabbage (B. oleracea, wild cabbage) was obtained from the NCBI protein sequence database (Accession number: ABS30827) [30]. The primary sequence of MYR protein from cabbage (B. oleracea) consists of 546 amino acids, which was used for modeling. Homology modeling was performed using the Iterative Threading ASSEmbly Refinement (I-TASSER) server based on the ab initio/threading method [31]. To predict the three-dimensional structure of cabbage (B. oleracea) MYR, 1MYR_A (Sinapis alba (S. alba) MYR) [32], 1E4M_M (S. alba MYR) [9], 1DWF_M (S. alba MYR) [33], 1VO2_A (Sorghum bicolor Dhurrinase) [34], and 4JHO_A (Oryza sativa beta-glucosidase) [35] structures were chosen as the templates. These template structures were selected by multiple threading approach Local Meta-Threading-Server (LOMETS) based on sequence similarity to cabbage (B. oleracea) MYR. Alignments of the query sequence with each template were generated using the ConSurf server. This program provides multiple sequence alignment with conserved features, which provides information about specific amino acid features in their local environment [36]. A total of five three-dimensional (3D) models were generated by I-TASSER; among them, the best model was identified based on confidence score (C-Score) [37].

Molecular Dynamics Simulation
Simulations of molecular dynamics (MD) were performed for the modeled MYR protein using the GROMOS96 43a1 force field implemented in GROningen MAchine for Chemical Simulation (GROMACS) version 4.6 [38,39]. A cubic box was constructed around the 3D model with a minimum of 10 Å between the box edges and the model surface. Appropriate counter-ions (Na + or Cl í ) were added for neutralizing this system. Further it was solvated with simple point charge (SPC) water molecules and minimized using the steepest descent method to remove spurious contacts. In addition, this system was equilibrated 100 ps under NVT (canonical ensemble) and NPT (isothermal-isobaric ensemble) conditions. The temperature of the system was set to 300 K after applying position restrains to the protein. This temperature and pressure was maintained using the Berendsen weak-coupling method. A cut-off radius of 0.9 nm was used in the simulation for van der Waals and electrostatic interactions [40]. The LINCS algorithm was used to reduce close contact between the complex system and bond constraints [41]. Finally, a 20-ns production was carried out for the MYR model using the particle mesh Ewald (PME) electrostatics method under NPT conditions. For analysis, final coordinates were saved every 2 ps.
The final MYR trajectory coordinates of the 20-ns simulation were selected, subjected to energy minimization. Further, the refined model was validated using Protein Quality Predictor (ProQ) [42] and VERIFY3D plot [43] and then used for molecular docking studies.

Analysis of MYR Active Site and N-Glycosylation Sites
The refined model was analyzed using Catalytic Site Atlas (CSA) database [16] and Q-SiteFinder [44] for active site identification. Those sites were compared with template structures for structural evaluation. Further N-glycosylation sites were identified using the NetNGlyc1.0 server [45].

MYR-Sinigrin Interaction
The two-dimensional (2D) structure of sinigrin (2-propenylglucosinlate) was retrieved from the NCBI PubChem database with the following accession number CID: 5464493. This molecule was drawn using ChemSketch (Advanced Chemistry Development, Inc. Toronto, ON, Canada), saved in a 2D (.mol) format, and then converted into a 3D (.pdb) format by importing into Discovery Studio 3.5 visualizer (DS 3.5) (DS, Accelrys, Inc. SanDiego, CA, USA). This molecule was optimized using the Conjugate Gradients method [46] followed by Steepest Descent in 200 steps using the PyRx program [47]. Universal Force Field (UFF) was used for minimization of the sinigrin molecule [48]. Molecular docking analysis was performed using AutoDock Vina [49]. Active site residues were carefully assigned to the binding pocket of the MYR protein, and docking steps were followed from previous studies [50]. The best docking orientation was identified from binding affinity score and hydrogen bond interaction to the active site based on visual inspection. The 2D graphical visualization of MYR and sinigrin interaction was generated by DS 3.5. The X-Score program [51] was used to calculate binding affinity scores of sinigrin and non-bonded interactions.

Complex Stability Refinement by Molecular Dynamics Simulation
The MYR-sinigrin docked-complex file obtained from docking studies was subjected to MD using GROMACS 4.6. The sinigrin structural topology file was generated using the Dundee PRODRG2 server [52] with chirality, full charges and energy minimization. In addition, MYR topology was created using the GROMACS pdb2gmx utility. The above-mentioned dynamics parameter was used to produce a 20-ns simulation for the modeled MYR and MYR-sinigrin complex. The detailed MD procedure was implemented according to our previous study [51]. In this analysis, root mean square deviation (RMSD), and root mean square fluctuation (RMSF) calculations were computed using the GROMACS utilities of g_rms, and g_rmsf, respectively. For this study, Chimera [53], and DS 3.5 were used for graphical visualization, and graphs were plotted using Microsoft Excel.

Homology Modeling
The full-length cabbage (B. oleracea) MYR protein consists of 546 amino acids including a glycoside hydrolase family 1 (Glyco_hydro_1) domain (41-522 amino acids). A search for homologous proteins for modeling cabbage (B. oleracea) MYR prompted us to use 1MYR_A, 1E4M_M, 1DWF_M, 1VO2_A and 4JHO_A protein structures as templates. Multiple sequence alignment was performed with these template proteins to identify sequence alignments of structural motifs (Figure 1). The alignment results revealed good conservation of cabbage (B. oleracea) MYR with the S. alba proteins, with 91% query coverage and 72% sequence identity. Cabbage (B. oleracea) MYR shared 40% identity with the protein from Sorghum bicolor (1VO2_A) and 39% identity with that of Oryza sativa (4JHO_A). In addition, three potential S. alba templates were found to share their common ancestor with cabbage (B. oleracea) MYR. The best model was selected based on C-Score value, and this score was used for estimating the quality of predicted models by I-TASSER. Then, the selected best model was subjected to MD Simulations. Figure 1. Sequence alignment of myrosinase with their template structures generated by using Chimera visualization tool. Gaps and conservation residues were denoted in dotted lines and gray color, respectively. In addition, block levels represented as identical (tallest blocks), conserved (Intermediate blocks) and mild conserved (Small blocks) and not identical (no blocks).

Model Refinement and Structure of MYR
The stability and MD properties of the constructed model were evaluated using MD simulations. To identify the energetically favorable structure for this model, MD simulation was observed up to 20 ns. The RMSD values were examined throughout the entire simulation time to check the structural stability of the tested model. RMSD values of MYR backbone atoms attained equilibration state in approximately 7.1 ns with an average of 0.43Å and remained stable throughout the course of the simulation (Figure 2A). RMSF analyses of backbone atoms of each amino acid are shown in Figure 2B; much structural variation were found in secondary structure. Numerous reports have confirmed that O-glycosidase family enzymes are linked to a single (ȕ/Į) 8 barrel [54]. A final snapshot from the 20-ns MD simulation was subjected to further evaluation and refinement. The stereo-chemical quality of the enhanced model was validated by Veryfy3D and ProQ. The Verify3D evaluation server reported that our optimized model scored over 80% of amino acids in the 3D/1D profile, which indicates residues are in favorable positions (Figure 3). ProQ revealed an LGscore of 4.808, which indicates that our conceptual MYR structure was an extremely good model. From these evaluation results, we confirmed that the developed model was satisfactory, and we therefore used it for subsequent analyses. Furthermore, this refined MYR protein was compared with template structures to identify sequence conservation. The alignment results revealed that our modeled Cabbage (B. oleracea) MYR shared with sequence conservation of 72.14%, 71.54%, 71.54%, 40.29%, and 39.13% with 1MYR_A, 1DWF_M, 1E4M_M, 1VO2_A, and 4JHO_A structures, respectively. In addition, the three templates from S.alba (1MYR_A, 1DWF_M, and 1E4M_M) were shown high similarity with refined MYR including Glyco_hydro_1 domain. The GOR4 analysis has revealed that our cabbage (B. oleracea) MYR model contains secondary structural elements of 22.71% helix, 20.51% strand, and 56.78% coil. Three disulfide bridges were found, in the Cys2-Cys456, Cys34-Cys454, and Cys226-Cys234 regions. The known myrosinase crystal structure from S. alba has shown similar structural features [26]. Diverse myrosinase proteins are involved in plant defense mechanisms, and although their origins vary, they are associated with common ancestors. Moreover, N-glycosylation sites are important to maintain the molecular stability and solubility of myrosinase in the dehydrated environment. N-glycosylation sites were predicted for the optimized MYR model using the NetNGlyc1.0 server. A total of ten Asn-Xaa-Ser/Thr sites were observed, at residues 110, 152, 177, 238, 329, 366, 380, 497, 518, and 533. Among these potential N-glycosylation sites, those with a threshold value over 0.5 were considered as potential sites and the residues were 110, 177, 329, 497, and 533. In comparison, ten sites were predicted in S. alba (six potential) and seven were reported from a theoretical model of Brassica juncea (B. juncea) [26]. The number of predicted N-glycosylation sites from cabbage (B. oleracea) was close to S. alba, in accordance with good amino acid conservations including potential sites (Supplementary Figure S1).

Molecular Interaction Studies
To better understand MYR from cabbage (B. oleracea), an interaction study between the substrate sinigrin and the enhanced MYR protein active sites was carried out using the Autodock Vina program. For docking studies, potential active site information was examined from the S. alba myrosinase crystal structure and confirmed with the CSA database and Q-site finder. Residue ARG115 was involved in the catalytic target site, and GLN207 and GLU427 had catalytic nucleophile activities. In addition, THR210, ASN348, and TYR350 were present in the binding pocket of MYR for substrate binding (Figure 4). These key amino acids were assigned as active sites in docking to identify productive binding orientations. Three-dimensional docked representations of the sinigrin molecule in the MYR active site and interacted 2D plot are shown in Figure 5. Further, myrosinase-sinigrin interaction was evaluated by binding affinity scores and hydrogen bond formation in the MYR active sites. In the docked state, the substrate sinigrin formed hydrogen bonds with the catalytic target site and two catalytic nucleophile residues, along with other residues of SER117, ASP221, and LYS285 ( Table 1). The docking program predicted that MYR-sinigrin binding affinity was í7.2 kcal/mol in the docked state. The nucleophile residue GLN207 was involved in the formation of one hydrogen bond to the HE21-O sinigrin group with hydrogen bond length of 2.10Å. A second nucleophilic residue, GLU427, formed two separate hydrogen bonds with the hydrogen atoms present in the sinigrin (H-OE1, and H-OE2) with hydrogen bond length of 1.85Å and 2.34Å respectively. In addition, one hydrogen bond was observed at the residue SER117 HG atom (HG-O) with 1.7Å hydrogen bond length. We also observed two hydrogen bonds with LYS485 (HZ2-O and HZ2-O) with lengths of 2.43Å and 2.41Å, respectively. Moreover, sinigrin contains seven rotatable bonds for enzyme binding, including one nitrogen (N) atom and sulfur group. This specific sulfur group was used to recognize the MYR binding and complex stability. In the docked state, residues GLU427 and LYS485 interacted with N atoms and the sulfur group respectively. In addition, the following hydrophobic interaction residues also contributed to molecular stability: TRP475, THR426, ASN206, PHE483, GLU482, SER56, ALA57, PRO161, TRP71, TYR58, and PRO223.

Stability Evolution of the MYR-Sinigrin Complex
An MD simulation was conducted to examine the MYR-sinigrin complex stability, and differences were monitored by RMSD with (holo) and without (apo) sinigrin. Accordingly, two individual MD simulations were carried out; the first contained only MYR as discussed above, the second consisted of the MYR-sinigrin complex file that was obtained from docking analysis. The protein mobility and structural changes were measured by RMSD plot of MYR backbone atoms throughout the 20-ns simulation. RMSD analysis of the MD trajectory noted that fluctuations were attained to 6.20 ns with an average of 0.26 Å, and then retained stability throughout simulation ( Figure 6A). Apo and holo forms of MYR attained equilibration at over 6 ns. RMSD values of the apo form ranged between 0.0005 and 0.4711 Å, which was higher than those of the holo form, which ranged between 0.0005 Å and 0.2918 Å during 20 ns trajectories. Further, RMSF values were calculated for the apo and holo forms of MYR to estimate local protein mobility. Analysis of both RMSF values against backbone atoms indicated that a few atoms highly fluctuated in the N-terminal and C-terminal ends. In addition, the backbone atoms of RMSF values were compared with MYR (apo) and MYR-sinigrin complex (holo) ( Figure 6B). Analyzing the trajectory files, the crucial residues had RMSF values ranging from 0 to 0.11 Å. In addition, compare to apo form of MYR, nucleophilic and catalytic residues namely GLN207, GLU427, and ARG115 were seen with slight fluctuations. The results from both apo and holo form of MYR showed that the residues in the sinigrin binding regions stabilized upon ligand recognition. Further, higher fluctuations in RMSF values were noted at N-terminal and C-terminal than other regions. The superimposition of initial and final structure revealed that sinigrin-binding orientation was well preserved and slight structural changes were observed in the active sites. Moreover the RMSF values of predicted ten N-glycosylation sites were ranged between 0.08 Å and 0.28 Å over 20-ns simulation time (Supplementary Figure S2). Analysis of potential N-glycosylation sites revealed that slight and higher fluctuations were observed at the residues of 110, 117, 329, and 497 and 533, respectively, when compared to apo and holo form of MYR. The final trajectory of the 20-ns simulation was examined to predict binding affinity and hydrogen bond stability. Hydrogen bond stability analysis showed that the nucleophile GLN207 maintained the hydrogen bond throughout the entire 20-ns simulation (Supplementary Figure S3). This nucleophilic residue is responsible for enzyme-substrate interaction [55,56]. Another active site residue, THR210, formed a hydrogen bond in the hydrophobic state. By contrast, the hydrogen bonds of residues SER117, THR221, and LYS485 disappeared compared to the docking state in long dynamic simulations. We calculated the binding affinity of MYR-sinigrin to be í6.98 kcal/mol using the X-Score program ( Table 2). In the case of non-bonded interaction, we observed negative logarithm dissociation constants of hydrophobic pair score (5.12), hydrophobic match score (5.13), hydrophobic surface score (5.11) and predicted mean binding affinity score (5.12). Both docking and MD simulation results suggests that the substrate sinigrin successfully recognized the MYR binding pocket at crucial active sites. In general, conceptual model has been developed based on homology principle and help to characterize their putative functions in preliminary way. In addition, this structure based homology prediction method has limitations at the structural level, including large protein, multi-complex, and multi-domain proteins [57]. The developed model will be helpful to characterize the myrosinase function and sinigrin interaction at molecular level but further experimental confirmations are required.

Conclusions
In this study, we constructed 3D model of MYR from cabbage (B. oleracea) based on sequence and structural similarity. The constructed model, which includes the Glyco_hydro_1 domain, was used for MD simulation to evaluate protein structural stability. The stereo-chemical quality of the optimized model confirmed that all residues were in favorable positions (Veryfy3D) and that the model was extremely good (ProQ). Our refined model included three disulfide bridges at the following positions: Cys2-Cys456, Cys34-Cys454, and Cys226-Cys234. In addition, ten conserved Asn-Xaa-Ser/Thr N-glycosylation sites were detected, among them five residues were potential sites, having a threshold value over 0.5. Notably, N-glycosylation sites, active sites and sequence conservation of our conceptual model showed similarity to the S. alba MYR crystal structure. Further, molecular interaction analysis showed two nucleophile (GLN207, and GLU427) and one catalytic (ARG115) residues were involved in the formation of hydrogen bonds along with hydrophobic interactions. The specific sulfur and nitrogen atoms present in sinigrin were well recognized by MYR active sites in the docked state. MD simulation results also suggested that strong hydrogen bond formation was maintained at residue GLU207 throughout the entire simulation, and binding orientation are the same as observed in the initial MYR-sinigrin complex structure. Our proposed 3D model of MYR and sinigrin binding orientation will be helpful to understand the MYR structural activities.


text:
Genetic Variants in MiRNA Processing Genes and Pre-MiRNAs Are Associated with the Risk of Chronic Lymphocytic Leukemia

Genome wide association studies (GWAS) have identified several low-penetrance susceptibility alleles in chronic lymphocytic leukemia (CLL). Nevertheless, these studies scarcely study regions that are implicated in non-coding molecules such as microRNAs (miRNAs). Abnormalities in miRNAs, as altered expression patterns and mutations, have been described in CLL, suggesting their implication in the development of the disease. Genetic variations in miRNAs can affect levels of miRNA expression if present in pre-miRNAs and in miRNA biogenesis genes or alter miRNA function if present in both target mRNA and miRNA sequences. Therefore, the present study aimed to evaluate whether polymorphisms in pre-miRNAs, and/or miRNA processing genes contribute to predisposition for CLL. A total of 91 SNPs in 107 CLL patients and 350 cancer-free controls were successfully analyzed using TaqMan Open Array technology. We found nine statistically significant associations with CLL risk after FDR correction, seven in miRNA processing genes (rs3805500 and rs6877842 in DROSHA, rs1057035 in DICER1, rs17676986 in SND1, rs9611280 in TNRC6B, rs784567 in TRBP and rs11866002 in CNOT1) and two in pre-miRNAs (rs11614913 in miR196a2 and rs2114358 in miR1206). These findings suggest that polymorphisms in genes involved in miRNAs biogenesis pathway as well as in pre-miRNAs contribute to the risk of CLL. Large-scale studies are needed to validate the current findings.

Introduction
Chronic lymphocytic leukemia (CLL) is the most frequent leukemia of adults in Western countries. It is characterized by the accumulation of monomorphic small round B lymphocytes in peripheral blood, bone marrow, spleen and/or lymph nodes [1]. CLL has been reported to have the highest genetic predisposition of all hematologic neoplasias [2], with an~8.5-fold increased relative risk in first-degree relatives [3]. However, previous genetic linkage studies failed to provide evidence that rare, high-penetrance genes contribute significantly to the familial risk [4], suggesting that additional loci of modest effects could be involved. Actually, recent genomewide association studies (GWAS) have identified several low-penetrance susceptibility alleles associated with CLL [5][6][7][8][9], associations that have been thoroughly validated in multiple independent series [10][11][12][13]. Nevertheless, these studies have focused their effort on identifying genes involved in the risk of CLL and scarcely study regions that are implicated in non-coding molecules, such as microRNAs (miRNAs).
MiRNAs are a group of small, functional, non-coding RNAs of approximately 20 nucleotides that regulate gene expression at the post-transcriptional level by binding to the 3´untranslated region (UTR) of a target gene [14]. This can lead to inhibition of translation or enhanced degradation of their target mRNAs. Through this mechanism, miRNAs can regulate a large number of genes. Consequently, variations in miRNA function can affect many targets with functional consequences.
To date several reports have described altered expression of miRNAs in CLL, suggesting their implication in the development of the disease [15][16][17]. MiRNA expression levels can be altered by genetic polymorphisms in pre-miRNAs and/or also in miRNA biogenesis genes. Additionally, genetic polymorphisms in miRNAs or in mRNA targets altering the mRNA-miRNA binding may be also relevant, affecting the miRNA function. Indeed, mutations in miRNA transcripts have been seen to be common in CLL [17,18], suggesting that CLL predisposition involves genetic alterations in miRNAs. Supporting this idea, a recent study has found that certain miR-SNPs show different frequency in CLL patients compared with a control population [19].
In addition, several studies have hitherto described polymorphisms in miRNAs associated with the susceptibility to different malignancies. For example, rs7372209 in miR26a-1 has been associated with the risk of throat and esophageal cancer and rs11614913 in pre-miR196a2 with breast, lung and gastric cancers, etc [14]. However, despite accumulating evidence that inherited genetic variation in miRNA genes and miRNA processing genes can contribute to the predisposition for cancer, the role of these in the risk of CLL has not been extensively studied. Therefore, the aim of this study was to evaluate whether polymorphisms in pre-miRNAs, and/ or miRNA processing genes, contribute to the risk of CLL.

Study population
The study population included a total of 523 Spanish individuals, 132 CLL patients previously reported [20] and 391 cancer-free controls. Both, cases and controls, were selected from the same residing area of the north of Spain and had West European ancestry. A brief description of the clinical characteristics of CLL population is described in S1

Genes and polymorphisms selection
A total of 118 SNPs were selected for the study. Of them, 72 SNPs were located in 21 genes involved in the miRNAs biogenesis and processing (S2 Table) and 46 SNPs were situated in 44 miRNAs (S3 Table). The SNPs in miRNA processing genes were selected after literature review and by using Patrocles database (http://www.patrocles.org). In each gene, using tagSNPs, we covered almost all the SNPs with potentially functional effects using bioinformatics tools (F-SNP, Fast-SNP, polymirTS [21,22] Patrocles [23], Hapmap, Haploview). We considered functional effects those causing amino acid changes, alternative splicing, those located in the promoter region in putative transcription factor binding sites, or disrupting/creating miRNAs targets. We also selected SNPs from the literature. Only SNPs with a reported minor allele frequency greater than 5% (MAF0.05) in European/Caucasoid populations were selected.
For the selection of SNPs in miRNAs, we had two factors in consideration. A) miRNAs can regulate a wide range of genes that are not completely defined. Therefore, any miRNA could be implicated in the regulation of genes affecting CLL risk. B) The number of miRNAs with polymorphisms described in the databases at the moment of selection with a MAF > 0.01 was affordable. Consequently, we selected all the SNPs in pre-miRNAs with a MAF>0.01 in European/Caucasoid populations, using Patrocles and Ensembl databases (http://www. ensembl.org/) (Welcome Trust Genome Campus, Cambridge, UK) and literature review. The preliminary list of SNPs was filtered using suitability for the Taqman OpenArray platform as criteria.

Genotyping
Genomic DNA was extracted from peripheral blood from controls and from granulocytes (normal cells) from CLL patients as previously described [20]. SNP genotyping was performed at the General Research Services (SGIker) of the University of the Basque Country using Taq-Man Open Array technology (Applied Biosystems), following the manufacturer's protocol. Data were analyzed with Taqman Genotyper software for genotype clustering and calling. As a genotyping control, 20 duplicate samples were placed across the plates.

Statistical analysis
Haploview 4.2 was used to search for any deviation of Hardy-Weinberg equilibrium in controls, determine haplotype block structure, calculate the genotyping success rate for each SNP and compare haplotype frequencies and clinical parameters between cases and controls. Significance of the best result was corrected for multiple testing by 1000 permutations. Genotype frequencies in the cases and controls were compared using a χ2 test. The effect sizes of the associations were estimated by the OR's from univariate logistic regression. The most significant test among the different genetic models was used to determine the statistical significance of each SNP. The results were adjusted for multiple comparisons by the False Discovery Rate (FDR) [24]. In all cases the significance level was set at 5%. Analyses were performed by using R v2.11 software.

miRNAs secondary structures prediction
The RNAfold web tool (http://rna.tbi.univie.ac.at) [25] was used to predict the most stable secondary structures of the miRNAs showing significant SNPs. The analyzed sequences included 50 bp upstream and 50 bp downstream of the pre-miRNAs.

Genotype association study
A total of 28 SNPs showed statistically significant association with CLL risk (S5 Table). After FDR correction, 9 SNPs remained significant. Among them, seven were located in six processing genes ( Table 1, Fig. 1) and two in pre-miRNAs ( Table 2). Two of the significant SNPs, rs3805500 and rs6877842, were located in DROSHA gene. Rs3805500 showed the most significant finding of the analyses under the log-additive genetic model (AA vs AG vs GG). The G allele showed a 0.54-fold decreased risk of CLL (95% CI: 0.38-0.77; P = 0.032). Other 5 SNPs in processing genes, rs1057035 in DICER1, rs17676986 in SND1, rs9611280 in TNRC6B, rs784567 in TRBP and rs11866002 in CNOT1, showed statistically significant results (P<0.05).
We also observed that two SNPs in pre-miRNAs were significantly associated with CLL risk, rs11614913 and rs2114358 in the miR196a2 and miR1206, respectively. The genotype TT rs11614913 was associated with a 2.3-fold increased risk under the recessive model (95% CI: 1.38-3.95; P = 0.002) while the genotypes AA and AG rs2114358 showed 0.52-fold protective risk to CLL under the dominant model (95% CI: 0.33-0.82; P = 0.005).
When the correlation between SNP alleles and clinical parameters were analyzed (Rai stage, morphology, levels of serum markers such as β2-microglobulin or lactate dehydrogenase, presence of genomic aberrations, marrow infiltration, presence of adenopathy, lymphocyte doubling time, residual minimum disease, and treatment) no significant association was found in CLL patients (S6 Table).

Haplotype Association study
We conducted haplotype analyses for the miRNA-related genes with at least two statistically significant SNPs in this study before FDR correction, including DROSHA, DICER1, DGCR8 and TNRC6B. Two haplotype blocks, one in DROSHA and another one in DGCR8, showed statistically significant differences in their frequency in cases versus controls (χ2 test; adjusted P-value<0.05 after permutation testing) ( Table 3). The significant haplotype of DROSHA gene, formed by rs3805500-rs7735863-rs6884823 (AGG), was associated with increased risk (OR = 1.63; 95%

miRNAs secondary structures
In order to check if SNPs in miR196a2 and miR1206 affected their secondary structure, we compared the minimum free energy for optimal secondary structures of both variants for each miRNA. The SNP rs2114358 in pre-miR1206 had a dramatic effect on its secondary structure, whereas the SNP rs11614913 in pre-miR196a2 had none (Fig. 2).

Discussion
In the current study, we identified genetic variants in miRNA processing genes and pre-miRNAs that were associated with the risk of CLL. We examined a total of 91 SNPs, finding nine significantly associated after FDR correction, seven in miRNA processing genes and two in pre-miRNAs. Interestingly, among the significant SNPs in miRNA processing genes, two were located in DROSHA; rs3805500 in the intronic region and rs6877842 upstream. Moreover, a haplotype including the SNPs rs3805500-rs7735863-rs6884823 (AGG) was also found significant in this gene. DROSHA encodes a type III Rnase enzyme, essential for the maturation of pri-miRNAs into pre-miRNAs [26]. Aberrant expression of DROSHA at mRNA or protein level was shown to be associated with patient survival and tumor progression in several types of cancer, including pleomorphic adenomas [27], esophageal cancer [28], cervical cancer [29] and ovarian cancer [30]. Moreover, the SNP rs3805500 which was the most significant finding in our study, was found by other authors to be associated with head and neck cancer recurrence [31]. In addition, a recent study in lung cancer has demonstrated that the polymorphism rs640831 in high LD with rs3805500 (r 2 >0.8) was associated with reduced DROSHA mRNA expression and with expression changes in 56 miRNAs out of 199 analyzed [32], some of which could also regulate genes involved in the development of CLL. Therefore, in the same way, the two polymorphisms detected in the present study could alter the expression of miRNAs involved in the regulation of genes involved in the risk of CLL. Additionally, five extra miRNA biogenesis pathway SNPs were also found significant. One of them, rs1057035, was located within the 3' UTR of DICER1, an important region for RNA stability. Moreover, it was located at the putative miRNA binding site of miR574-3p, which could affect the binding of this miRNA to DICER. DICER plays an important role in the cleavage of pre-miRNAs into their mature form. A recent study showed that rs1057035 C allele was associated with a lower expression of DICER and a decreased risk of oral cancer, which seemed to be due to a higher inhibition of miR574-3p on DICER mRNA [33]. In our study, in contrast, the variant genotype CC was associated with a 2.7-fold increased risk of CLL. These differences between studies suggest that DICER may have different roles in the tumorigenesis process depending on the cancer type. Supporting this idea, some studies have associated reduced expression of DICER with the development of lung cancer [34], while others have seen overexpression of DICER in oral cancer [35]. Therefore, it seems plausible that rs1057035 is affecting the binding of miR574-3p which in turn could be altering the regulation of DICER and the risk of CLL.
Another significant SNP, rs784567, was located in the 5' near gene region of TRBP, an integral component of a DICER-containing complex important to the cytoplasmic processing of pre-miRNAs into mature miRNAs. Previous studies also found rs784567 to be borderline significant in bladder cancer [36] and oral premalignant lesions [37]. Moreover, depletion of TRBP was found to result in significantly impaired miRNA biogenesis [38]. Consequently, genetic variation affecting the expression of this gene could lead to global changes in the miR-NAome and subsequent deregulation of genes involved in CLL origin.
On the other hand, rs9611280 was situated in TNRC6B, gene that mediates miRNA-guided mRNA cleavage [39]. According to in silico analyses (F-SNP), this missense SNP was not predicted to abort protein function, but instead affected the regulation at splicing level. This change could perturb the maturation of the mRNA, changing the sequence of the protein, and therefore contribute to carcinogenesis in CLL.
Finally, the last two significant SNPs identified were rs17676986 at SND1 and rs1186002 at CNOT1, genes involved in RNA-induced silencing complex (RISC). Genetic variants in both genes have been previously associated with autism [40] and ventricular arrhythmias [41] respectively, but this is the first time that are associated with cancer.
In addition to SNPs from miRNA processing genes, there were also two SNPs in miR196a2 and miR1206 significantly associated with CLL risk after FDR correction. miR196a2 rs11614913 variant homozygote TT showed an increased risk of CLL. Different rs11614913 genotype variants were previously reported to be associated with cancer risk; for example, the variant CC was associated with lung [42], hepatocellular [43], glioma [44], and gastric [45] cancer, while the variants TT or TT + CT were associated with breast [46], prostate [47], and head and neck cancer [48]. Despite these differences among studies, probably due to the different populations used in the analyses [49], it seems probable that rs11614913 plays an important role in cancer development. rs11614913 is located in the 3p mature miRNA region, so it could affect maturation of 5p and 3p miRNAs, as well as the expression and function of these miR-NAs. Actually, Hoffman et al already demonstrated that rs11614913 C allele led to less efficient processing of the miRNA precursor to its mature form as well as diminished capacity to regulate target genes in breast cancer [46]. Interestingly, potential targets of miR196a2 include genes with relevance in CLL, such as CDKN1B or HMOX1. Therefore, alterations in the sequence of miR196a2 which may affect miRNA function, could contribute to CLL risk due to the deregulation of its targets.
Finally, miR1206 rs2114358 AG-GG genotypes showed a reduced risk of CLL. Remarkably, miR1206 is located at 8q24.2 a frequent site of chromosomal breaks in cancer and a signature of genomic instability [50]. Moreover, miR1206 is of particular interest as it is part of long non-coding RNA (lncRNA) transcript of the PVT1 gene, locus of both tumorigenic translocations and retroviral insertions [51], which has been also shown to play a role in cell proliferation and apoptosis [52]. Overexpression of PVT1 transcript has been documented in a variety of tumor tissues [53] and miR1206 expression has been found in increased levels in B cell tumors [54]. The allelic effect of miR1206 rs2114358 has been further analyzed in different human cell lines but without concluding results, suggesting that the effect of this SNP was cell line and cell type specific [55]. We found that miR1206 rs2114358 AG-GG genotypes showed a reduced risk of CLL, so it is likely to expect that AA genotype was associated with an increased risk. Moreover, we found by in silico analysis (RNA fold), that A allele altered the miRNA secondary structure, suggesting that this variant alter the maturation of miR1206.
A recent study found that rs56103835 in pre-miR323b was associated with CLL when compared with European samples from the 1000 genome project [19]. We also analyzed this SNP in our population, finding no significant association (S5 Table), neither when CLL population was compared with the European samples from the 1000 genome project (0.143 versus 0.187), nor when our control population was used (0.143 versus 0.183).
In conclusion, the findings of the present study indicated that SNPs located at DROSHA, DICER1, SND1, TNRC6B, TRBP and CNOT1 genes, involved in miRNAs biogenesis pathway, as well as in miR196a2 and miR1206 may contribute to the risk of CLL. To our knowledge, this is the first study analyzing specifically miRNA-related SNPs in CLL, which opens a promising approach to search for new susceptibility markers in CLL. New large-scale studies including functional analyses will help to validate our findings.
Supporting Information S1

Author Contributions
Conceived and designed the experiments: AGO IMG ELL. Performed the experiments: IMG AGC NBA MPG. Analyzed the data: AGO IMG AGC. Contributed reagents/materials/analysis tools: MA. Wrote the paper: AGO IMG.