Brain-spinal cord interaction in long-term motor sequence learning in human: An fMRI study

The spinal cord is important for sensory guidance and execution of skilled movements. Yet its role in human motor learning is not well understood. Despite evidence revealing an active involvement of spinal circuits in the early phase of motor learning, whether long-term learning engages similar changes in spinal cord activation and functional connectivity remains unknown. Here, we investigated spinal–cerebral functional plasticity associated with learning of a speciﬁc sequence of visually-guided joystick movements (sequence task) over six days of train- ing. On the ﬁrst and last training days, we acquired high-resolution functional images of the brain and cervical cord simultaneously, while participants practiced the sequence or a random task while electromyography was recorded from wrist muscles. After six days of training, the subjects’ motor performance improved in the sequence compared to the control condition. These behavioral changes were associated with decreased co-contractions and increased reciprocal activations between antagonist wrist muscles. Importantly, early learning was characterized by activation in the C8 level, whereas a more rostral activation in the C6-C7 was found during the later learning phase. Motor sequence learning was also supported by increased spinal cord functional connectivity with distinct brain networks, including the motor cortex, superior parietal lobule, and the cerebellum at the early stage, and the angular gyrus and cerebellum at a later stage of learning. Our results suggest that the early vs. late shift in spinal activation from caudal to rostral cervical segments synchronized with distinct brain networks, including parietal and cerebellar regions, is related to progressive changes reﬂecting the increasing ﬁne control of wrist muscles during motor sequence learning.


Introduction
Over the past decades, a plethora of neuroimaging studies has sought to provide insights into the neurofunctional mechanisms mediating mo-Abbreviations: MSL, motor sequence learning; SRTT, serial reaction rime tasks; CNS, central nervous system; SEQ, sequence; RND, random; BOLD, blood oxygenation level dependent; EPI, echo-planar imaging; PNM, physiological noise modelling; FoV, field of view; tSNR, temporal signal to noise ratio; PPI, Psychophysiological Interactions; SCT, spinal cord toolbox; EMG, electromyography; FCR, flexor carpi radialis; ECR, extensor carpi radialis; AAR, automatic artifact removal; IAJ, integrated absolute jerk; RMS, root mean square; ACC, anterior cingulate cortex; SMA, supplementary motor area; SPL, superior parietal lobule; CB, cerebellum; ANG, angular gyrus. movements come to be performed effortlessly as a unitary sequence through repeated practice. Ample evidence indicates that MSL (a form of procedural memory) takes place over different phases (e.g., fast [minutes to hours] and slow [hours to days]), with different patterns of cerebral neuroplasticity associated with each (e.g., ( Dayan and Cohen, 2011 ;Doyon and Benali, 2005 ;Doyon et al., 2018 ;Lohse et al., 2014 )). Specifically, changes in neuronal activity in the cortico-cerebellar and corticostriatal circuits (in conjunction with hippocampus) have been observed from the fast to the slow phase as the motor sequence is being acquired and consolidated ( Lohse et al., 2014 ;Doyon et al., 2009 ;Doyon et al., 2003 ;Lehericy et al., 2005 ).
Yet, an important question, largely overlooked by the neuroimaging community so far, is whether MSL relies on brain plasticity alone, or it involves other parts of the central nervous system (CNS), such as the spinal cord. Traditionally, the spinal cord has been viewed as a mere pathway for relaying signals from higher-order brain regions to the peripheral effectors, and from the sensory cells back to the brain, with a limited role in information processing ( Bizzi et al., 2000 ). Yet, a large body of evidence, primarily from animal neurophysiology research using traditional learning paradigms and lesions studies, has revealed that spinal neurons do support habituation and sensitization learning mechanisms, as well as Pavlovian and instrumental conditioning; all representing proxies of spinal cord neuroplasticity ( Grau, 2014 ;Wolpaw, 2007 ). Similarly, there is also accumulating evidence from human electrophysiology studies showing that spinal motoneurons excitability (as measured via the Hoffmann H-reflex) can change following a single motor skill training session, hence supporting the hypothesis that intrinsic spinal cord plastic changes can be observed during this form of motor learning ( Albuquerque et al., 2018 ;Meunier et al., 2007 ;Lungu et al., 2010 ;Roche et al., 2011 ;Christiansen et al., 2017 ). Indeed, the first direct evidence that MSL relies on neuroplasticity at different levels of the CNS, and that the spinal cord contributes uniquely to this process, has been provided by our group ( Vahdat et al., 2015 ). In the latter study, we used an innovative functional magnetic resonance imaging (fMRI) protocol, during which we scanned the brain and cervical spinal cord simultaneously, while participants performed an MSL task ( Vahdat et al., 2015 ). This simultaneous brain-spinal cord fMRI protocol allowed us to assess, in vivo , the MSL-related activity at multiple CNS levels, and to probe the functional interactions that develop between the spinal and supraspinal levels during motor learning. Importantly, the results from this study, not only revealed an increase in blood oxygenation level dependent (BOLD) spinal cord activity in C6-C8 spinal segments, when comparing the execution of a complex vs. simple sequence of movements, but also for the first time, to demonstrate that this functional change was independent from that of supraspinal sensorimotor structures known to be anatomically connected to the spinal cord or involved in MSL. Finally, our findings also showed that the initially strong functional connectivity between the spinal cord and sensorimotor cortices diminish with learning, while the connectivity between the spinal cord and cerebellum increase over the course of the acquisition process ( Vahdat et al., 2015 ). Altogether, our results strongly suggested that during the fast MSL learning phase, the spinal cord contributes distinctively from the brain to the motor learning process.
In addition to its role during the early motor skill learning phase, there is evidence that the spinal cord may also contribute to the acquisition of a more proficient motor program, following additional practice during the slow learning phase. Several animal-model studies have indeed provided direct support to this long-term spinal cord plasticity hypothesis ( Wolpaw, 2007 ;Levine et al., 2012 ;Peters et al., 2017 ;Ueno et al., 2018 ). In humans, indirect evidence of long-term spinal or corticospinal plasticity has also been provided through several neurophysiological studies that revealed a decrease of amplitude in the H-reflex and motor evoked potentials in motor skill experts like athletes and dancers ( Mynark and Koceja, 1997 ;Nielsen et al., 1993 ;Ryder et Moscatelli et al., 2016 ). Yet, despite these results, there is a lack of conclusive and direct evidence in humans that the spinal cord is involved in the slow learning phase of a new motor skill due primarily to the fact that testing using these electrophysiological paradigms is done offline, and not while subjects are performing the skilled motor behavior of interest. Moreover, unlike fMRI, these electrophysiological techniques do not provide a precise localization of the changes in spinal activity, nor an assessment of these changes and their interaction on a large scale at both spinal and supraspinal levels.
To fill this knowledge gap, we thus employed our pioneering simultaneous brain/spinal cord fMRI data acquisition approach ( Vahdat et al., 2015 ;Cohen-Adad et al., 2010 ;Finsterbusch et al., 2013 ;Vahdat et al., 2020 ) in young healthy subjects in order to: ( Dayan and Cohen, 2011 ) test the hypothesis that the spinal cord plays an active role in the slow acquisition phase of a novel sequence of movements, and ( Doyon and Benali, 2005 ) assess changes in the spinal cord -brain functional connectivity occurring during the early and late phases of motor sequence learning.

Participants
Thirty young, healthy subjects were recruited to participate in this study. We included only right-handed individuals and excluded those with a history of neurological or psychiatric disease, motor-system complications, use of neurological medication or with any MRI-incompatible object or material in their body. We also excluded individuals who previously participated in motor-learning experiments, as well as those with previous training in playing an instrument for more than three consecutive years in the last five years. Of the participants fulfilling all the inclusion and none of the exclusion criteria, we subsequently excluded the data from two subjects who did not come for the training sessions following the first day of scanning, and from another subject who had excessive movements (more than 7 mm in plane) during the neuroimaging session while executing the motor task, hence resulting in failed motion correction. Thus, the final sample considered for analysis consisted of 27 participants (14 females, mean ± SD age = 24.85 ± 2.98 years old). The ethics committee at the centre de Recherche de l'Institut Universitaire de gériatrie de Montréal (CRIUGM) reviewed and approved the study (protocol number: CMER-RNQ 15-16-06). All participants gave written informed consent at the beginning of the experiment and were debriefed and compensated for their participation after having completed the experiment.

General procedure
The study consisted of 6 testing sessions, which occurred daily over a period of 6 days. Participants were scanned twice; at the beginning and at the end of a 6-days motor skill learning regimen. On Day 1, participants were acquainted with the set-up by practicing the MSL task inside a mock MRI scanner. The sequence of movements practiced in the mock scanner was different from the other sequence used during the actual experiment. Following the participant's habituation with the experimental set-up and task, they were then moved into the actual 3.0T MRI scanner, where they were asked to perform both the sequence and random conditions of the motor task (please, see the task description below for more details). For the next four behavioral-only sessions, each taking place on subsequent days, participants were placed inside the mock MRI scanner and asked to perform the sequence learning condition (SEQ, see below). Finally, on Day 6, subjects were again tested (at the same time of the day as the Day 1) and scanned following the same procedures as the one used in the first session, except for the habituation period.

Fig. 1. (A)
Illustration of the sequence of movements that subjects were required to practice during each session. Subjects were asked to perform the task using an MRI compatible joystick that they manipulated with their (dominant) right hand. (B) The shimming and scanning parameters were adjusted for each sub-volume (brain, spinal cord) separately described in ( Finsterbusch et al., 2013 ). The brain sub-volume covered the whole cerebrum (30-33 slices), while the spinal cord sub-volume covered the C3-T1 segments of the spinal cord (8-10 slices).

Motor sequence learning paradigm
Participants performed the motor task using a joystick manipulated by their dominant (right) hand ( Fig. 1 -A). The task began with the visualization of a cursor (i.e., a black circle) in the center of the screen, which showed the location of the joystick corresponding to the neutral position of both the joystick and the wrist, and of 4 static targets located clockwise at 3, 6, 9 and 12 o'clock, corresponding to wrist abduction, flexion, adduction and extension, respectively. Participants were instructed to move the cursor to the location of a target as soon as the circle was filled (by moving the joystick toward the corresponding location). As soon as the cursor touched the target's location, the latter disappeared, and the next target appeared at another location among the 3 remaining targets. Participants were asked to reach targets as quickly as possible, but unbeknown to them, a repeating sequence of 8 locations was repeated ten times in a block of 80 trials. There was a short break (20 s) between two consecutive blocks and each session comprised a total of 15 blocks. Two different conditions of the motor task were administered. In one condition, the targets appeared in a sequential order that was repeated every 8 trials (8-item sequence learning condition; SEQ), whereas in the other condition, the targets appeared in a pseudorandom non-sequential order (random condition; RND). The specific sequence ( Fig. 1 A) was selected such that each target was hit twice, and no transition was repeated twice. We also made sure that, on average, the random condition followed the same distribution of targets and transitions. Same number of blocks and trials used in each session. Overall, the average time of the participants to complete the motor sequence training on day 1 was around 30 min while it reduced to less than 25 min after 6 days of practice.

MRI data acquisition
To investigate reciprocal influences of the brain and the spinal cord during motor learning, we used a specific MRI pulse sequence that enabled us to acquire BOLD data simultaneously from the cervical spinal cord and brain ( Finsterbusch et al., 2013 ;Tinnermann et al., 2017 ). Specifically, this sequence allowed us to acquire data from 2 sub-volumes of images within the same TR, one covering the cervical spinal cord (C3-T1 spinal segments), and one covering the entire brain, with different field of view and spatial resolutions ( Fig. 1 -B). At the beginning of the scan session, five sample echo-planar imaging (EPI) volumes were acquired, and the acquisition parameters for each of these 2 sub-volumes were optimized to ensure optimal quality data for both the brain and cervical spinal cord. For this purpose, shimming properties and resonance frequencies were adjusted individually for each of the 2 sub-volumes prior to the start of the test runs in each scanning session and were then updated during BOLD signal recording in the brain and the spinal cord. Details of this shimming procedure are discussed elsewhere ( Tinnermann et al., 2017 ). In addition and based on those pre-scan acquisitions, a specific "z-shim" slice-specific approach was used to maximize the signal intensity within the cervical spinal cord ( Finsterbusch et al., 2012 ). Due to the size of the area, and to reduce the likelihood of displacement and artifacts related to magnetic susceptibility differences between bone and soft tissue, a small rectangular shim box (red box, Fig. 1 B) was placed around the cord.
A 3T TIM Trio Siemens MR scanner (Siemens Healthcare, Erlangen, Germany) equipped with a 12-channel head coil paired with a 4-channel neck coil (both receive only) was used for the current study. Participant's head and neck were supported with extra foam pads to minimize movements during scanning. The lower edge of the head coil (along with the upper lip of the participant, approximately C2-C3 vertebrae) were aligned to the magnet isocenter for optimal coverage.
For functional imaging, a total number of 40-43 EPI slices (depending on the participant's head size) were acquired in descending-order for both sub-volumes (interleaved within each volume). The upper subvolume covering the brain contained 30-33 transversal slices. The whole brain of most participants fitted within this volume, except for a few subjects for whom the box was tilted upward approximately 10°(around the X-axis), thus cutting a small part of the upper motor cortex, specifically the part that corresponds to the lower limbs, which were not involved in the current task. Slices in the brain sub-volume had a field of view (FoV) of 220 ×220 mm 2 , with 2 × 2 * 3.5 mm 3 voxel size and a 3.5 mm gap between slices. The flip angle was 90°, the bandwidth per pixel was 1515 Hz yielding an echo spacing 0.74 ms. We used a 7/8 partial Fourier encoding and parallel imaging (GRAPPA) with an acceleration factor of R = 2 and 48 reference lines to minimize the echo time, hence resulting in an echo train length of 35.6 ms and an echo time of 30 ms for the brain. Parallel imaging also helped reducing image distortions.
The lower sub-volume included 8-10 slices covering the cervical spinal cord from C3 to T1 spinal segmental levels and was oriented parallel to the spinal cord at the C6 level. Slices in the lower sub-volume had an FoV of 132 ×132 mm2, with 1.2 × 1.2 × 5 mm3 voxel size and 9 mm gap between slices. The flip angle was 90°, the bandwidth per pixel was 1263 Hz yielding an echo spacing of 0.90 ms. 7/8 partial Fourier and parallel imaging ( R = 2, 48 reference lines) was again used, resulting in an echo train length of 43.3 ms and an echo time of 33 ms. Finally, the TR for the measurement of all slices (brain + spinal cord) was 3140 ms, except for 3 subjects who had a TR of 3050 ms or 3200 ms (depending on each participant coverage within the field of view). Both brain and spinal cord sequences included a fat saturation pulse.
Furthermore, two saturation pulses were applied in a V-shaped configuration to minimize ghosting and inflow artifacts related to blood flow in the major cervical vessels. To reduce the level of noise during the functional acquisition runs, signal from the head coil was only used for the brain images, while signal from the neck coil was only used for the spinal cord images, using specific reconstruction procedure described elsewhere ( Finsterbusch et al., 2013 ). The average ( ± SD) temporal signal to noise ratio (tSNR) within the brain and spinal cord masks were 48.1 ± 13.3 and 12.5 ± 4.1 respectively.
For structural scanning, we employed two imaging sequences. First, we used a 3D-MPRAGE T1-weighted sequence with the following characteristics: 175 sagittal slices covering entire head and neck down to the C8/T1 vertebrae, TR = 2.3 s, TE = 3.45 ms, FA = 9°, TI = 1.1 s, FoV = 192 ×240 ×320, with a voxel resolution of 1 × 1 * 1 mm 3 . In addition, we used a T2-weighted MEDIC sequence with the following parameters: 12 transversal slices, slice thickness = 5 mm (no gap), in-plane resolution = 0.6 × 0.6 mm 2 , TR = 464 ms, TE = 22 ms, FA = 20°, bandwidth 180 Hz per pixel, FoV = 184 ×256 mm 2 , GRAPPA acceleration factor = 2, reference lines = 48. The latter slices were placed at the same position to the ones used for the spinal fMRI sequence. The T2 image enabled us to acquire a proper GM/WM contrast in the spinal cord and to improve our segmentation and registration outcomes. We also recorded the participants' heart rate (using pulse oximeter) and breathing rhythms during each, the latter measures being used for physiological noise modeling for the spinal cord fMRI data analyses using the Physiological Noise Modelling (PNM) toolbox from the FSL platform (see supplementary information for more details).

Scanning procedure
After moving to the real scanner and following the localizer and the initial shim process, the first of the 2 structural volumes (T1) was acquired, followed by 2 runs of fMRI data. The order of the 2 runs (RND or SEQ conditions) was counterbalanced between subjects. Subjects were not informed explicitly about the experimental conditions to prevent possible expectation bias interference. The functional runs were then followed by the second structural imaging (T2s) scan. The procedure used for the second scanning session (Day 6) was the same as the one employed in the first scan session, except that there was no additional practice of the motor sequence inside the mock-scanner.

MRI data processing
Image processing was done separately for the brain and the cervical spinal cord using the spinal cord toolbox (SCT-v3.1; ( De Leener et al., 2017 )), FSL-6.0 ( Jenkinson et al., 2012 )) and in-house Matlab-2016a programs. Details about Image processing and the analysis of MRI data is provided in supplementary information. Briefly, spinal cord fMRI processing was started with the removal of the first two volumes, and then motion correction was done. The spinal cord was segmented in the motion-corrected images and then registered to the MNI-Poly-AMU template by first registering the image to T2 and then back to the template. After smoothing and high-pass filter, the data was used along with Physiological Noise Modelling (PNM) regressor for the analysis. Brain fMRI data was started with the removal of the first two volumes, then slice timing correction and motion correction. Motion-corrected images were registered to the MNI-152 template using the high resolution T1 images. High-pass filter and smoothing were done before the analysis.

EMG recording and preprocessing
EMG signals were recorded using MRI compatible surface electrodes from two muscles, flexor carpi radialis (FCR) and extensor carpi radialis (ECR), while subjects were in the scanner performing the task. EMG was recorded using a Biopac EMG100C-MRI smart amplifier with the sampling rate equal to 2500 Hz (gain = 500). Online noise filtering was done using 2 analogue filters for the FCR and ECR muscle activity (Fig. S1A). Offline preprocessing of the EMG signal and automatic artifact removal (AAR) were carried out by estimating the periodicity of MRI gradient artifacts on the raw EMG time series, and then using a sliding average template of these artifacts to subtract them from the signal during each run.
The raw time series were first up-sampled ten times, from 2.5 to 25 kHz, to facilitate the temporal interpolation of the signal. To compensate for temporal drift of the artifact over time, we then detected the local maximum for each TR, and computed the peak intensity in the time series that corresponded to the slice excitation. We performed a temporal interpolation of the signal such that the number of data points between 2 peaks were constant and equal to the real TR. Once the time series had the exact same number of points for each TR, we were then able to compute an artifact template. To account for the small fluctuations in MRI artifacts over time, we used a sliding average of 21 TRs to compute local templates instead of calculating a global template of the artifact and subtracted these from the signal. Finally, we down-sampled the resulting time series from 25 kHz back to its original sampling frequency (i.e., 2.5 kHz). Figure S1 presents snapshots of the raw EMG signal, as well as of cleaned (denoised) signal.

Performance and EMG analysis
In order to use the rich kinematics/trajectory data recorded from joystick and assess learning-related changes in movement kinematics, we used a measure of jerkiness instead of simply calculating the reaction times for hitting targets ( Balasubramanian et al., 2015 ). To assess the participants' performance in each block, a measure of jerkiness of movement for reaching the next target was calculated from the third time derivative of pointer displacement (or the second derivative of velocity). The absolute value of the jerk at each time point was calculated and then integrated over time throughout the task to calculate the integrated absolute jerk (IAJ) ( Goldvasser et al., 2001 ). The mean IAJ scores in the first 7 blocks and the last 7 blocks were calculated separately to assess learning within the run of 15 blocks, with the hypothesis that the movements will become smoother in the SEQ condition as a result of learning.
For the EMG data preprocessing, full-wave rectification was applied on the AAR-cleaned EMG signal, followed by low-pass filtering (Butterworth filter, fco = 10 Hz) to obtain a measure corresponding to the envelope of EMG signal. For each muscle and condition, the envelope signal was then divided by a normalization factor to account for EMG electrode placement differences across days. This factor was obtained by taking the median value of the root mean square (RMS) of the EMG in a 500-ms sliding window during RND blocks in each session, given that the movement velocity in the RND condition was not different across sessions. We then calculated measures of the "wasted ", "effective " and "total " signals between the pair of antagonistic muscles (FCR and ECR) using their normalized envelope EMG signal at each time point t, as follows ( Banks et al., 2017 ): Equation 1 .
We then estimated the level of co-contraction between FCR and ECR ( Banks et al., 2017 ) and their reciprocal activity ( Flanders and Cordo, 1987 ), as follows: Equation 2 .
where RMS is the root mean square operator calculated during each movement element. Then, for each condition and session, the cocontraction and reciprocal activation values were averaged across 4 movement elements along the direction of wrist flexion/extension, where the FCR and ECR muscles mainly act as antagonists. Furthermore, as EMG activity, and hence EMG power, is closely related to movement speed ( Hof et al., 2002 ), EMG signals were additionally normalized using the mean movement speed in each movement block in order to account for differences in speed across blocks and conditions. Hence, we measured the EMG power in each trial by calculating the RMS of the normalized envelope signal during each individual movement element (from one target to the next one). For each muscle and condition, the EMG power was then averaged across all movement elements in that condition. The average EMG power between FCR and ECR muscles were also calculated. Note that normalization through the use of the movement speed does not affect the calculation of co-contraction and reciprocal activation, as the speed normalization factor would be canceled out in the numerator and denominator, as indicated by the formulas in Equation 2.

Statistical analysis
For the purpose of presentation, the mean ± standard error of the mean (SEM) is reported. A repeated measure ANOVA with Blocks (1:7 vs 9:15) * Day (Day 1 vs. Day 6) * Conditions (SEQ vs. RND) as the withinsubject factors was run on the mean IAJ. Paired sample t-tests were used to compare the first seven blocks with the last 7 blocks of each run within each session. Normality and equality of variance were checked for each analysis and the degree of freedom was adjusted when needed. Alpha level was set to 0.05 for all analyses. we noticed that it could be explained by the following two-way interactions: that is between block and condition [F(1,26) = 25.66; p < 0.001; n2p = 0.50], and between day and condition [F(1,26) = 10.39; p = 0.003; n2p = 0.75], indicating that the learning effect (operationalized as significant performance changes in the Sequence, but not in the Random condition) was evident within and between each practice day.

Behavioral results
Paired sample t-tests comparing the first and the last 7 blocks of each day and condition (see Fig. S2) showed that subjects' trajectories for reaching a target were significantly less jerky in the last 7, as compared to the first 7 blocks of practice on each day during the SEQ condition [SEQ-D1: t(26) = 4.63, p < 0.001; SEQ-D7: t(26) = 4.67, p < 0.001], whereas there was no such significant difference in the RND condition [RND-D1: t(26) = − 0.23, p = 0.82; RND-D7: t(26) = 1.52, p = 0.14]. This indicates that learning-related changes (i.e., reducing jerkiness with practice) were evident in the SEQ condition in each training session, where participants could make use of the sequential pattern to perform smoother movements, but not in the RND condition when stimuli appeared at unpredictable locations. This interpretation is further supported by the fact that, on Day 1, participants' trajectories for reaching a target in the SEQ condition were significantly less jerky than those in the RND condition [t(26) = 9.49, p < 0.001]. Similarly, the trajectories for reaching a target in the SEQ condition on Day 6 were also significantly smoother than those in the RND task [t(26) = 14.85, p < 0.001]. In addition, there was a significant improvement in jerkiness (a reduction in mean IAJ) when comparing the subjects' trajectories in the SEQ condition across days (less jerky on Day 6 than Day 1) [t(26) = 3.25, p = 0.003], whereas there was no such difference between days in the RND condition [t(26) = 1.10, p = 0.28].

Neuroimaging results
We had 2 main goals in the current study: first, to examine spinal cord's involvement in motor sequence learning and second, to characterize the changes in functional connectivity between the brain and spinal cord associated with the short vs. long-term phases of motor learning.

Neural correlates of short-term and long-term motor sequence learning: within session effects
As expected, an analysis contrasting the results of the SEQ vs. RND conditions on Day 1, hence reflecting within-session short term MSL functional changes (Fig. S8A), yielded clusters of activation in the C8 segment of the spinal cord. The peak was located ipsilateral to the hand used in the motor sequence learning task. Results of the same contrast at the brain revealed clusters bilaterally in the anterior insula and the SMA, as well as the right anterior cingulate cortex (ACC) and lobule V of the right cerebellum (Table S5-A).
The same contrast assessing the within-session changes specific to long-term MSL on Day 6 ( Figure S8-B) revealed 3 clusters of increased activity more rostrally than the one found on Day 1, that is between C6-C8 segments of the spinal cord, with peaks located both ipsilateral and contralateral to the hand used for the task. The same contrast at the brain level, yielded activation bilaterally in the anterior cingulate cortex, the right cerebellum lobule IV, and the left cerebellar lobule V (Table S5-B). This pattern of results suggests a switch of activity from more caudal regions of the spinal cord on Day 1 to more rostral regions of the spinal cord on Day 6 of learning. At the cerebral level, the Day 1 to Day 6 transition resulted in more activity in cerebellar regions such as ACC and cerebellum and less activity in superior parietal lobule (SPL) and Insula.

Neural correlates of short-term and long-term motor sequence learning: between session effects
Comparing the results of the short vs. long-term MSL (SEQ > RND and Day 1 > Day 6; Fig. 3: blue), the analysis yielded an activation cluster in the spinal cord at the level of the C8 segment with the peak on the right (ipsilateral) side. The same contrast at the brain level resulted  in activations in the superior parietal lobule bilaterally (peak on the left), the anterior cingulate cortex, the right superior parietal lobule, and the right lateral occipital cortex (Table S6-A).
The reverse comparison (SEQ > RND, Day 6 > Day 1; Fig. 3: red) in the spinal cord resulted in 2 activation clusters: one in the C5 and on the right ventral side and one in C6 in the center with the peak on the left side on the cord. The same contrast applied at the brain level produced activation clusters in the right putamen, the right posterior cingulate, the right occipital pole, and the right cerebellum lobule IX and vermis lobule VIII.

Neural correlates of long-term sequence practice
We contrasted the activation across days (Day 6 vs. Day 1) for the SEQ condition only (SEQ D6 > SEQ D1; Fig. S9: Red-Yellow) to examine the neural correlates of sequence practice only. This analysis resulted in clusters of activity in the right posterior cingulate and cuneal cortices, as well as in the inferior part of the right superior parietal lobule. The same contrast in the spinal cord yielded activation in the rostral part of the C7 segment (on the left side, dorsally). The reverse contrast for the SEQ condition (SEQ D1 > SEQ D6; Fig. S9: Blue-Light Blue) resulted in significant activation clusters in the right superior parietal lobule, the globus pallidus, and the right inferior temporal gyrus at the brain level, whereas it yielded activation in the rostral part of the C8 segment of the spinal cord (on the right side and ventrally). In line with the findings of earlier analyses, here we also noticed changes in the activation foci within the spinal cord from the caudal to more rostral segments as the motor sequence was practiced over many days. We computed correlations between behavioral learning-related changes in motor learning (i.e., reduction in jerkiness of the movement) and activation values in the brain and spinal clusters in order to provide further evidence that these functional neurophysiological changes are related to motor learning. To this end, we calculated the correlation of the between-session difference in mean jerkiness with the beta values extracted from the clusters of activation in the brain and spinal cord, resulting from the contrast between the sequence practice on Day 1 and Day 6 (figure S9). In the brain, there was a positive correlation between the beta values extracted from SPL and the improvement in performance on the sequence task on Day 1 (P(27) = 0.73, p < 0.001). By contrast, in the spinal cord, there was a single positive correlation between the beta values extracted from the cluster surviving the Day6 > Day1 contrast and the improvement in performance for the sequence condition on Day 6 (P(27) = 0.81, p < 0.001).

Learning dependent interaction between the spinal cord and the brain
Furthermore, we investigated the changes in functional interaction between the local spinal network that was activated during the course of learning and the cerebral network involved in motor sequence learning using PPI analysis (see the description in methods). Here we report the results of the PPI analysis for the learning effect (SEQ > RND) on Day 1 and Day 6 separately. Our results revealed that on Day 1, activity found in the C8 segment of the spinal cord was correlated with that of the left precentral gyrus, the left superior parietal lobule, the brainstem, and right cerebellum (crus I and crus II) in the SEQ compared to the RND condition ( Fig. 4 , yellow-red; GRF corrected, p < 0.01). Furthermore, spinal cord activity on Day 6 (levels C6-7) was correlated with that of the right angular gyrus and left cerebellum (crus I) in the SEQ compared to the RND condition ( Fig. 4 , blue-light blue; GRF corrected, p < 0.01).

Learning-related changes in EMG activity
As changes in the spinal cord activation might be closely related to/manifested by EMG activity, we computed 3 different scores to assess changes in FCR and ECR muscle activities (the main muscles involved in the task) as a result of motor learning. First, we calculated the strength of reciprocal activation between FCR and ECR muscle pair during flexion and extension movements, that is when they act as antagonistic muscles (see Equation 2, Methods). This measure of reciprocal activation was differently modulated in the SEQ and RND conditions depending on the practice day (significant interaction, day (Day 1, Day 6) by condition (SEQ, RND) repeated measures ANOVA, F (1,26) = 5.99, p = 0.02) ( Fig. 5 A). Pairwise comparisons revealed higher reciprocal activation in the SEQ condition compared to RND condition on Day 6 ( t (1,26) = 3.5, p = 0.002), as well as in the SEQ condition on Day 6 compared to Day 1 ( t (1,26) = 2.47, p = 0.02). By contrast, the strength of reciprocal activation did not differ between the SEQ and RND conditions on Day 1 ( p > 0.4), nor between the RND conditions on Day 1 and Day 6 ( p > 0.5).
Next, we calculated the measure of co-contraction of FCR and ECR muscle pair during flexion and extension movements (see Equation 2, Methods). As for the measure of reciprocal activation, the amount of co-contraction was also differently modulated in the SEQ and RND conditions depending on the practice day (significant interaction, day (Day 1, Day 6) by condition (SEQ, RND) repeated measures ANOVA, F (1,26) = 6.75, p = 0.01) ( Fig. 5 B). Pairwise comparisons further revealed decreased levels of co-contraction in the SEQ condition compared to RND condition on Day 6 ( t (1,26) = 3.76, p < 0.001), as well as in the SEQ condition on Day 6 compared to Day 1 ( t (1,26) = 2.6, p = 0.01). Again however, the level of co-contraction between the SEQ and RND conditions on Day 1 ( p > 0.2), or between RND conditions on Day 1 and Day 6 ( p > 0.5) did not differ.
Finally, we calculated the average EMG signal power of FCR and ECR muscles during task performance (see Equation 2, Methods). For this analysis, EMG was further normalized by movement speed to account for differences in speed across days and conditions (see methods). A two-way repeated measures ANOVA with the normalized EMG power as dependent variable showed no significant interaction between Day and Condition ( p > 0.6) (Fig. S10). Also, the normalized EMG power did not differ between Day 1 and Day 6 (no main effect of Day, p > 0.9), or between SEQ and RND conditions (no main effect of Condition, p > 0.2). Overall, our EMG analyses showed that long-term motor learning through multiple days of practice changes the activation pattern in antagonistic muscles, while the average EMG power per unit of movement velocity remains unchanged.

Discussion
Here, we investigated the neural correlates mediating the early (fast) and late (slow) motor sequence learning through the simultaneous acquisition of fMRI data in the brain and cervical spinal cord. Analyses of the behavioral data revealed that subjects learned the new sequence of movements as they showed improvements in performance within and between training sessions, hence reflecting both short-and long-term motor learning, respectively. The between-sessions behavioral changes were also paralleled by similar changes in EMG metrics, with reduced co-contractions and increased reciprocal EMG activity during sequential movements from the early-to the late-learning phases. Importantly, fMRI results at the spinal cord and brain levels revealed both intra-and inter-sessions learning-related changes in neural activity. As predicted, early learning related activations in the spinal cord were present in a small area located caudally in the C8 segment, ipsilateral to the hand executing the task. By contrast, performance of the learned sequence in the late-learning session was accompanied by widespread and more rostral activations in 3 clusters across C6 and C7 spinal segments. Parallelly, increased cerebral activity was observed in inferior frontal and temporal regions as well as in the insula, caudate and cerebellum (Crus Co-Contraction. For both variables, post-hoc t-tests revealed a significant difference across days only for the SEQ condition and a significant difference between the Sequence (SEQ) and Random (RND) conditions only on Day 6. ( * p < 0.05; * * p < 0.005). Error bars represent 95% Confidence Interval. I) within the early learning phase, while cerebellar activations in lobules IV and V as well as in the anterior cingulate cortex, cuneus and amygdala were observed in the late-learning session. Between-sessions contrasts revealed increased activations in early vs. late learning in the superior and medial frontal regions, precuneus, inferior parietal, as well as left caudate nucleus. Yet as expected, the opposite contrast (late vs. early) yielded learning-related increased activations in the right putamen, the right posterior cingulate, the right occipital pole, and the right cerebellum lobule IX and vermis lobule VIII. Finally, the connectivity analysis showed that in earlier stages of learning (Day 1), activity in the right cerebellum, left precuneus, and left precentral increased their functional connectivity with the spinal cord ROI (C8). By contrast, during the later stage of learning (Day 6), the increased connectivity with the late-learning related spinal ROI (C6) was observed in the angular gyrus and cerebellum.

Behavioral findings
Consistent with the literature using variants of Serial Reaction Time tasks (SRTT, ( Lungu et al., 2010 ;Nissen and Bullemer, 1987 ), our behavioral results demonstrate that participants' performance improved both within-as well as across training sessions. This suggests that such improvement was due to learning of the pattern of sequential movements, above and beyond the mere general effect of practicing move-ments needed to reach targets with a joystick. Across sessions, we also observed further improvements that were specific to the SEQ condition only, hence showing evidence that participants also reached the late MSL phase. In our joystick reaching task, both RND and SEQ performances became faster from Day1 to Day 6 practice, however, the most significant change was manifested in the kinematics data, as calculated by IAJ scores, which evaluates how smoothly the trajectories were performed. Thus in this task, long-term changes in performance due to learning appeared to be better characterized using kinematic measures, rather than through the subjects' performance speed.

Electromyographic findings
The changes in EMG measurements, which paralleled those seen in the behavioral output, are likely to reflect plasticity in the corticospinal circuits and within the spinal cord. This assertion is based on results from a previous electrophysiological study by our group ( Lungu et al., 2010 ), in which we used the same joystick-based SRTT task and H-reflex measurements (assessing activity in the spinal circuitry) acquired pre-, during and post-learning of a new sequence of finger movements. In the latter study, the H-reflex amplitude decreased significantly more in the SEQ, compared with the RND condition, until subjects reached asymptotic performance (i.e., a plateau), indicating that this was related to the sequential content of the task, and not simply a consequence of execut-ing wrist movements. Given that we employed the same task in the current study, it is appropriate to assume that the reduced co-contractions and increased reciprocal EMG activity, are indicators of changes of activity in the spinal and corticospinal circuits as well. Indeed, previous electrophysiology studies that assessed the corticospinal excitability through H-reflex conditioned by transcranial magnetic stimulation over the primary motor cortex (M1), during and after motor skill acquisition, have reported increases in corticospinal excitability ( Wiegel and Leukel, 2020 ). Studies that have assessed only the local spinal changes through H-reflexes of muscles involved in executing a motor skill to be learned also found decreases in unconditioned reflex amplitude immediately after learning, both in the upper ( Giboin et al., 2020 ;Kitano et al., 2009 ), as well as the lower limbs ( Mazzocchio et al., 2006 ). Interestingly, and in line with the long-term changes in EMG in our study, the decreases in unconditioned H-reflex amplitude reported elsewhere were found to be persistent 24 h after learning ( Giboin et al., 2020 ). Finally, studies that used conditioned H-reflexes have found reduced homosynaptic H-reflex depression ( Winkler et al., 2012 ) and increased presynaptic inhibition ( Kitano et al., 2009 ), both findings providing evidence for local spinal cord plasticity following motor skill learning. Altogether these findings, and those reported here, suggest that the spinal cord plays a role in the stabilization and facilitation of movement execution, which is achieved through modification of muscle and joint stiffness as the motor skill is being practiced.

Cerebral correlates of short-and long-term motor sequence learning
There is a large body of evidence that MSL is associated with different brain activation patterns following short-and long-term practice ( Doyon et al., 2018 ;Lohse et al., 2014 ;Doyon et al., 2009 ;Lehericy et al., 2005 ;Doyon, 2016 ;Albouy et al., 2013 ;Ungerleider et al., 2002 ). We found a decrease of BOLD activity in the caudate, but an increase of activity in the right putamen from the early to the late learning session; a pattern that has been described in a previous meta-analysis of brain activations associated with long-term motor learning paradigms ( Lohse et al., 2014 ) as well in numerous studies from our laboratory ( Doyon and Benali, 2005 ;Doyon et al., 2009 ;Doyon et al., 2003 ). Short-vs. long-term learning has also been associated with an increase of activity in the precuneus, which is thought to be involved in the integration of visuospatial information when attention is required for completing a cue-based task ( Nardo et al., 2011 ), the superior parietal lobule, known to be involved in the spatial integration and motor learning ( Coynel et al., 2010 ), and lateral occipital cortex which is linked to the observation of targets when learning assessed through the use of a SRTT paradigm is required ( Rothkirch et al., 2018 ). These findings are in accordance with the model of motor skill acquisition which proposes a 'specialization' of cortico-striatal and cortico-cerebellar networks as a function of the type of motor skill acquired ( Doyon et al., 2009 ). Yet, contrary to findings of some previous studies ( Lohse et al., 2014 ), we observed an increase in activity in the cerebellum when comparing late vs. early learning sessions. Although conjectural, this apparent discrepancy can be explained by the fact that movements needed to control the joystick have a strong motor adaptation component, in addition to the sequential aspect of the task. The cerebellum is known for its involvement in motor coordination and our previous study has also suggested increased connectivity between the spinal cord and the cerebellum ( Vahdat et al., 2015 ). Likewise, our results show an increase in hippocampal activity in late vs. early learning sessions, in apparent contrast with the view proposing that this structure is thought to be particularly implicated in the early, rather than late phase of MSL ( Albouy et al., 2013 ). We believe, however, that this inconsistency may be due again to the spatial nature of the task, which is expected to elicit this structure.

Spinal correlates of short-and long-term motor sequence learning
A plethora of studies have reported evidence of learning-dependent plasticity in the spinal cord following the acquisition of new motor skills ( Wolpaw, 2007 ). Most of these studies employed electrophysiological techniques (e.g., H-reflex) to show changes following the practice of simple motor control tasks or motor learning paradigms ( Lungu et al., 2010 ;Evatt et al., 1989 ;Carp et al., 2006 ). Previously, we have provided evidence of the spinal cord's independent contribution to MSL using fMRI and a finger tapping task ( Vahdat et al., 2015 ). Despite differences in the task used for measuring MSL, our current results replicate the findings in our previous study, as we show that short-term motor sequence learning measured using a target-reaching task is associated with an increase of activity in the C8 segment of the spinal cord, an area very similar to the one that was activated during the finger MSL task and one that is known be involved in the control of fine finger and hand movements ( Vahdat et al., 2015 ).
More importantly, however, in addition to replicating our previous results, the current findings expand them in two critical ways. First, thanks to the spatial resolution achieved with our fMRI pulse sequence, we were now able to determine the lateralization of the activity (ipsilateral to the hand used for the task), and to demonstrate that the peak of activation was in the ventral horn of the spinal cord, where motoneuron cell bodies are located. Additional activation at the center and contralateral side relative to the hand may suggest activation of local spinal interneurons at the commissural level ( Hanna-Boutros et al., 2014 ). Future studies with higher resolutions can help us to delineate this observation. Second, the present pattern of results shows that similar to the brain, the activation pattern reflecting long-term motor learning in the (C5/C6) spinal cord is distinguishable from that of the short-term motor learning (C7/C8). The C5-C8 spinal segments correspond to myotomes related to the innervation of upper limb muscles ( Martini et al., 2018 ). Interestingly, the EMG data suggest a decrease in co-contractions between FCR and ECR muscles. This is thus concomitant with the increase in the fine motor control exerted by the wrist flexors and extensors as learning progresses across multiple sessions. The C5-6 is a known hub for most ECR/FCR motoneurons while C7-8 mostly innervates extrinsic finger muscles. A more rostral activation associated with long-term learning is in line with studies suggesting that in motor learning paired with visual feedback, the motor commands involves ECR/FCR muscles during wrist flexion and extension than other muscles like those involving the fingers ( Danion et al., 2003 ;Ambike et al., 2016 ;Hirose et al., 2020 ). Accordingly, these changes suggest that with more practice, participants develop a more refined motor control by changing the synergy between muscles of the hand and wrist. The optimal feedback theory assumes the existence of redundant anatomical (muscles and joints), kinematic (velocity and trajectory of movement) and neurophysiological (multiple motoneurons innervating different fibers in the same muscle and motoneurons of different muscles receiving common inputs) alternatives to complete a specific action ( Pierrot-Deseilligny and Burke, 2012 ;Scott, 2004 ). Accordingly, the completion of the same series of actions in a task can be done by changing the synergy between different groups of muscles that could differ substantially from one to another, but still, achieve the task requirements ( Sohn et al., 2013 ). The effects of longterm practice of the motor sequence may have helped the system to elicit the most cost-effective activation pattern of the hand and wrist muscles, combined. The changes in activation from caudal to more rostral segments might thus be related to progressive changes of muscle synergies during learning, with more focal and fine motor control of wrist muscles (motoneurons located at C5-C6 segmental levels) with practice, involving spinal interneurons (spinal cord intermediate zone), and less co-activation of hand muscles (motoneurons located at C7-C8).

Brain-spinal cord interaction in early and late motor sequence learning
Another main contribution of the current study was the assessment of brain-spinal cord interaction during early and late MSL, made possible by the simultaneous fMRI of the brain and cervical spinal cord. The brain regions that were synchronized with those of the spinal cord involved during short-term learning (C8) were the precuneus, thought to be involved in the visuospatial coordination of motor action ( Wenderoth et al., 2005 ), and the motor cortex in the contralateral hemisphere, as early in the learning of a motor skill the motor cortex is actively involved in the control of every single motor units at the execution site ( Wymbs et al., 2012 ). By the end of the learning on Day 1, synchronization between the spinal cord, the motor cortex and the precuneus was no longer significant, and was replaced by an increase in synchronization of activity in the spinal cord and the angular gyrus, a region known to be involved in spatial attention ( Chen et al., 2012 ). On Day 6, however, we then observed an increase in the synchronization between the spinal cord and angular gyrus, which may suggest that after 5 sessions of practice, attentional mechanisms become more important for the completion of a motor sequence, while for the random practice, we did not see much change on Day 6. Finally, at the end of day 6, activity in the spinal cord and the cerebellum became more synchronized suggestive of a better performance in the SEQ condition, which is a rhythmic practice, synchronization between the cerebellum and the spinal cord is also required ( Arshavsky et al., 1986 ).
Research on the brain-spinal cord interaction in the acquisition of motor skills is mostly driven by studies in animal models. Those studies suggest the existence of differential connections between the interneurons located in the posterior and anterior horn of the spinal cord and the sensory and motor cortices ( Levine et al., 2012 ). Using a combination of optogenetics, electrophysiology and behavioral assessments in primates, it has been shown that this organization and the interaction between the two subsystems is crucial for skilled motor behavior: blocking the activity of neurons in the motor cortex and a subgroup of spinal interneurons resulting in impaired movement reaching, while damage to the sensory cortex and another subset of spinal interneurons impaired the release of movements in a goal-oriented task ( Ueno et al., 2018 ). The only previous study which used fMRI to investigate motor learning and the interaction between the brain and spinal cord has supported the idea behind the unique contribution of the human spinal cord to motor learning ( Vahdat et al., 2015 ). Yet it should be noted that our previous study focused only on the early stage of the learning and showed similar pattern of connectivity increase between the spinal cord and cerebellum as observed on Day 1 of the current study. In current study, we provide behavioral, electrophysiological and neuroimaging evidence for the long-term spinal cord plasticity related to motor skill acquisition in humans. However, given the organization of the central nervous system and the anatomical hierarchy between the brain and cord, it would be important to study the causal relationship between these structures during different phases of learning. The PPI analysis does not examine effective connectivity (one approach that allows the statistical estimation of causality) between the brain and the spinal cord. Other methods for estimating effective connectivity in fMRI data might be helpful in deciphering the pattern of causal relationship between the brain and spinal cord signals.
Furthermore, there are a few limitations in the current study that need to be considered when designing future studies. Similar to most of neuroimaging studies of motor learning, our sample consists of only right-handed participants, and thus this may have an impact on the lateralisation observed in our results. In addition, because of the nature of block-design GLM analysis on whole spinal cord data, a segmentation of the gray matter was not employed in the current analysis. Besides, we were only able to record EMG from two groups of muscles. Future studies that include more recorded muscle activity or use high-density EMG recording might give better understanding of changes in muscles synergies and their relation to spinal cord activities from early to late motor learning.

Conclusion
To our knowledge, this study is the first to use fMRI to investigate the spinal cord plasticity and its interaction with the brain associated with short-and long-term learning of a novel motor skill. We demonstrate that, long-term motor learning in the human spinal cord is distinguishable from short-term motor learning, and that the changes in the spinal cord interaction with the brain is related to the requirements of the task at-hand as revealed by the behavioral indicators of motor learning (jerkiness of movements). In a joystick task that relied on the integration of information from visual and motor networks, the increases in movement fluency associated specifically with the long-term acquisition of the motor sequence were associated, not only with specific changes within the spinal cord activity from caudal to rostral segments, but also with a shift in its cerebral connectivity from the sensorimotor network initially, to the attentional network and other brain areas involved in the control of rhythmic motor actions.

Author statement
Ali Khatibi was involved in data collection, data processing, analysis, and manuscript preparation.
Shahabeddin Vahdat was involved in design, data collection, analysis, manuscript preparation.
Ovidiu Lungu was involved in design, data collection, analysis supervision, manuscript preparation.
Jurgen Finsterbusch was involved in MRI sequence development/optimization and manuscript preparation.
Christian Büchel was involved in MRI sequence development /optimization and manuscript preparation.
Julien Cohen-Adad was involved in MRI sequence optimization, spinal cord MRI data processing optimization, and manuscript preparation.
Veronique Marchand-Pauvert was involved in design, EMG data processing and analysis, and manuscript preparation.
Julien Doyon was involved in design, supervision (data collection and analysis), and manuscript preparation.