Experience-related alterations in white matter structure and gene expression in adult rats

White matter (WM) plasticity during adulthood is a recently described phenomenon by which experience can shape brain structure. It has been observed in humans using diffusion tensor imaging (DTI). However, it remains unclear which mechanisms drive or underlie WM plasticity in adulthood. Here, we combined DTI and mRNA expression analysis and examined the effects of somatosensory experience in adult rats. Somatosensory experience resulted in differences in WM and grey matter structure. C-FOS mRNA expression, a marker of cortical activity, in the barrel cortex correlated with the structural WM measures, suggesting that WM plasticity is activity-dependent. Analysis of myelin-related genes revealed higher myelin basic protein expression in WM, while genome-wide RNA sequencing analysis identified 134 differentially-expressed genes regulating proteins involved in functions related to cell proliferation and differentiation, neuronal activity modulation and regulation of myelination. In conclusion, the macroscale measures of WM differences identified in response to somatosensory experience are supported by molecular evidence, which strongly suggest myelination as, at least, one of the underlying mechanisms.


ABSTRACT (160 words)
White matter (WM) plasticity during adulthood is a recently described phenomenon by which experience can shape brain structure. It has been observed in humans using diffusion tensor imaging (DTI). However, it remains unclear which mechanisms drive or underlie WM plasticity in adulthood. Here, we combined DTI and mRNA expression analysis and examined the effects of somatosensory experience in adult rats. Somatosensory experience resulted in differences in WM and grey matter structure. C-FOS mRNA expression, a marker of cortical activity, in the barrel cortex correlated with the structural WM measures, suggesting that WM plasticity is activity-dependent. Analysis of myelin-related genes revealed higher myelin basic protein expression in WM, while genome-wide RNA sequencing analysis identified 134 differentially-expressed genes regulating proteins involved in functions related to cell proliferation and differentiation, neuronal activity modulation and regulation of myelination. In conclusion, the macroscale measures of WM differences identified in response to somatosensory experience are supported by molecular evidence, which strongly suggest myelination as, at least, one of the underlying mechanisms.

INTRODUCTION (306 words)
There is accumulating evidence that structural changes in white matter (WM) occur in response to changes in experience, even during adulthood (Sampaio-Baptista and Johansen-Berg, 2017). For example, neuroimaging studies have reported experience-induced structural WM plasticity in both humans (Scholz et al., 2009;Hofstetter et al., 2013) and rodents (Blumenfeld-Katzir et al., 2011;Sampaio-Baptista et al., 2013). However, it is unclear how such macroscale changes relate to the underlying molecular mechanisms of experiencedependent WM plasticity during adulthood.
Structural changes in barrel cortex, such as synapse and spine formation (Trachtenberg et al., 2002) have been extensively reported in response to whisker stimulation , deprivation (Holtmaat et al., 2006) and somatosensory learning (Kuhlman et al., 2014), along with corresponding changes in expression levels of genes, such as BDNF (Rocamora et al., 1996) and synaptophysin (Ishibashi, 2002). More recently, somatosensory enrichment in middle-aged rats resulted in an increase in oligodendrocyte numbers and integration in barrel cortex (Hughes et al., 2018). However, it is not clear whether somatosensory experience also induces WM plasticity or whether such structural plasticity is localized only to cortical regions.
Here, we investigated if somatosensory experience in adult rats induces WM plasticity at the macro and the molecular scale, by combining neuroimaging and mRNA expression analysis. Adult rats were trained in a texture detection task (TDT) (von Heimendahl et al., 2007) and WM structure was assessed with diffusion tensor imaging (DTI). We used a recently developed analysis tool, Permutation Analysis of Linear Models (PALM) (Winkler et al., 2016) implemented in the FMRIB software library (FSL) (Smith et al., 2004), which allows inference over multiple modalities or measures of imaging and nonimaging data, to assess the joint and individual contribution of the 4 DTI measures.
Additionally, we analysed mRNA expression to further support the structural findings and identify putative molecular mechanisms underlying experience-dependent WM plasticity.

ANIMALS
All behavioural experiments were conducted at Radboud University Nijmegen (The Netherlands). The experiments were approved by the Animal Ethics Committee of the Radboud University Nijmegen (The Netherlands), according to Dutch legislation and all procedures were performed according to the project and personal licenses held by the experimenters. 48 animals (3 months old, male Long Evans rats (250 -450 g) (Harlan, Bicester, UK)) were housed in standard laboratory conditions under a 12-h light/12-h dark cycle at 20 ºC temperature and 40 -70 % humidity. The animals were housed individually for more precise control of their general welfare and because group housing may interfere with the task experience. All animals were given appropriate time to acclimate after delivery (1 week minimum) and had ad libitum access to food and water. After this period, they were handled daily for one week before the start of the task. Before the task, animals were exposed to the testing arena for 10 minutes each day for 1-2 days under dim visible light.
The animals were given no access to water for a period of 24h before the first session.
From here on, they received water during the task (0.1 ml of water per correct trial) and water was also made available ad libitum for 30 min after the session. The delay between the end of the task and the time period when water was freely available varied between 30 minutes to 2 hours in order to prevent the animals from learning that water would be available after the testing period.

TEXTURE DETECTION TASK
The texture detection task (TDT) is based on a previously described task (von Heimendahl et al., 2007). Rats were trained to distinguish between a smooth and a textured surface using operant conditioning as described below. Training was performed in the dark to avoid the influence of visual cues on performance. Potential olfactory cues were removed from textures by washing them at least once every individual animal session, and by using different sets of identical textures that were interchanged randomly between animals and sessions.
Rats were tested individually. During testing, the animal was placed on a 30 cm elevated platform with two water dispensers on each side. Under this platform was a small bridge where the animal could place its front paws for a short period of time in order to reach the stimulus presented in front of the platform. The stimulus consisted of a series of rectangle shapes with patterns that could be varied depending on the animal's performance in the task.
During the shaping period the animal was placed on the apparatus, and every appropriate response was rewarded. First, water was randomly delivered in order for the animal to learn where the water was placed. After the animal had learned this, it was rewarded for leaning on the edge of the platform and reaching the stimulus. Finally, water was delivered when the animal touched the texture with its whiskers.
During the training phase of the TDT, the stimulus was either a smooth texture (reference texture) or a positive copy of sandpaper on a resin material. Each animal was trained with a fixed association (e.g., turn left on rough, right on smooth). Only if it approached the correct drinking spout, the animal was given a water reward (0.1 ml per reward); for an incorrect choice, it received no water. The next trial started with a delay of 5 s. Between trials, the texture's stand was turned about its vertical axis by a computercontrolled stepping motor, which allowed for quick, randomized, and automated switching between textures. Each session lasted for about 30 min, during which the animal performed between 60-100 trials.

EXPERIMENTAL DESIGN
Animals were randomly assigned to the TDT and control groups, balancing for weight to obtain equal weight averages between the groups.
TDT group: For the TDT, animals (total n=28) were trained to distinguish between the reference texture (smooth) and a P100 texture (162 um average particle diameter). Rats were trained 5 days a week. Individual rats were sacrificed the day after they reached criterion (2 sessions at > 80% accuracy) on the P100 texture, in order to have comparable performance levels between animals. This resulted in a variation in the number of task exposure days which was then controlled for in subsequent analyses as described below. After reaching performance criterion in the P100 texture, a subgroup of rats (n=8) continued training with finer textures in a stepwise manner: P150 (100 um average particle diameter), P220 (68 um), P280 (52.2 um), P360 (40.5 um), P400 (35 um), P500 (30.2 um) and P600 (25.8 um). When the animals performed above criterion (> 80% accuracy) for a given texture, the rough texture was changed to a finer one on the following training day. Rats were sacrificed the day after they reached > 80% accuracy on the P600 texture.
Negative control experiments were performed in this subgroup of rats to demonstrate that the animals distinguished smooth vs rough textures and not other sensory attributes of the task (e.g. noises or odours). To do that, animals were presented with the same texture (P400 versus P400) and their performance was assessed. When rats were presented with the same texture their performance accuracy dropped to chance levels ( Supplementary Fig. 1). Active Control group: 12 rats were matched to an individual in the TDT group. Rats were water restricted and exposed to the TDT task for the same period of time as the matched animal. However, these animals were rewarded randomly, and not in relation to textureresponse contingencies. They received a similar number of rewards as the matched animal throughout the entire training period and were sacrificed after the same number of training sessions as the matched animal in the learning group.
Passive Control group: Caged controls (n = 20) were handled and weighed daily; their body weight served as a reference body weight with respect to the other group.

BRAIN PREPARATION
Rats were sacrificed by direct decapitation without anaesthesia on the day after they reached criterion. On the day of sacrifice, animals were trained on their respective task for 15 minutes, then placed back in their home cage for 15 minutes, after which they were sacrificed and the brains were removed. The passive control group was handled for 15 min prior to the sacrifice. The right hemisphere was frozen on dry ice and kept at -80 °C for molecular analysis and the left hemisphere was immersed in 4% PFA for DTI acquisition.
For DTI acquisition, all left brain hemispheres (n=48) were placed into falcon tubes (50ml) in pairs (one from each group), one hemisphere above the other, and embedded in 2 % agarose gel (Sigma-Aldrich) (Sampaio-Baptista et al., 2013). The hemispheres were aligned to each other along the posterior -anterior axis.

MRI ACQUISITION
All 48 ex-vivo left brain hemispheres were scanned in pairs overnight with a 4.7T

TISSUE DISSECTION
A subset of right brain hemispheres (n=24; 15 from the TDT group and 9 from the passive control group) were dissected for qPCR and RNA-SEQ. The experimenter was blind to the group for all the following procedures. All procedures were performed under RNase-free conditions. The brain hemispheres were sliced into 300 µm coronal sections using a cryotome (Leica GmbH, Germany) at -15 °C and mounted on glass slides. Cytochrome oxidase-stained reference sections were used as a template to locate the barrel cortex, following stereotactic coordinates (Paxinos and Watson, 1998). Punches of the barrel cortex and in WM directly underneath it of the right hemisphere were taken using a 1.20-mm micropunch (Harris Inc., UK) and stored at -80 °C before RNA isolation took place.

RNA ISOLATION
The experimenter was blind to the group in all the following procedures and the groups were randomly distributed to ensure equal distribution of groups to avoid any technical bias.
Samples were homogenized with a TissueLyser (Retsch GmbH, Germany) in TRIzol® Reagent (Invitrogen Co., USA). RNA was isolated with TRIzol® Reagent (Invitrogen), according to the manufacturers' protocol. The procedure was modified for small amounts of tissue by using 800 µl of TRIzol® Reagent and adding 1 µl of glycogen (Fermentas Inc., USA). RNA concentration and quality was determined with a NanodropTM ND-1000 spectrophotometer (Thermo Fisher Scientific Inc., USA) and 1% agarose gel electrophoresis, respectively. The samples were kept at -80 °C until further analysis.

SYNAPTONEUROSOME PREPARATION
Synaptoneurosomes were prepared by the method described by (Williams et al., 2009), with some modifications. Brain tissue punches were homogenized with a Teflonhomogenizer (12-14 strokes at 1000 rpm) in 4 mL of homogenization buffer, containing 0.35 M sucrose pH 7.4, 10 mM 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES), 1 mM ethylenediaminetetraacetic acid (EDTA), 0.25 mM dithiothreitol, 8 U/ml RNAse inhibitor and a protease inhibitor cocktail (Roche). Cell debris and nuclei were removed by centrifugation at 1000 g for 10 min at 4°C yielding pellet P1 and supernatant S1. The S1 fraction was passed sequentially through a series of filters with decreasing pore sizes of 80, 40 and 10 µm (Millipore). The final filtrate was centrifuged at 2000g for 15min at 4°C yielding pellet P2 and supernatant S2. Pellet P2 containing synaptoneurosomes was resuspended in 200 µL of homogenization buffer. Enrichment of synaptic components in the synaptoneurosomal fraction was assessed by western blot in control experiments. To assess cfos expression in synaptoneurosomes, RNA was isolated using the Trizol method, followed by downstream QPCR.

QUANTITATIVE PCR (qPCR)
Prior to cDNA-synthesis, 0.5 µg of each RNA sample was treated with 2 U DNase (Fermentas Inc., USA), in the presence of RiboLockTM RNase Inhibitor (20 U/µl) (Fermentas Inc., USA). For cDNA synthesis, through random priming, the RevertAidTM H Minus First Strand cDNA Synthesis kit (Fermentas Inc., USA) was used, following the manufacturer's guidelines. Prior to analysis, each cDNA sample was diluted 1/15 with MilliQ water. qPCR reactions were performed with the Rotor-Gene 6000 Series (Corbett Life Science Pty. Ltd., Australia). For each reaction, 2.5 µL of each diluted sample of cDNA was added to a mix containing 6.25 µL 2X MaximaTM SYBR Green qPCR Master Mix (Fermentas Inc., USA), 1 µL of each primer (5 µM) and 1.75 µL MilliQ water. Primers were designed using NCBI Primer-Blast (www.ncbi.nlm.nih.gov/tools/primer-blast/) and synthesized at Sigma-Aldrich (UK). Cycling conditions were 10 min 95°C followed by 40 cycles of 15 sec at 95°C, 30 sec at 60°C and 30 sec at 72°C. After cycling, a melting protocol was performed, from 72 °C to 95 °C, measuring fluorescence every 1 °C, to control for product specificity.
For the candidate gene analysis of the WM, the following genes were selected for their role in myelin and WM plasticity (PLP1, OLIG1, NOGO-A, MBP, MOG, MAG, PDGFRA, NKX2.2, NKX2.6, SOX10, MOBP, ARC). The c-FOS gene was selected to confirm Barrel Cortex neuronal activation in response to the TDT task. Relative expression of the selected genes of interest was calculated after obtaining the corresponding Ct values and correcting for unequal sample input using geNorm (Vandesompele et al., 2002), which identifies the two most stably expressed housekeeping genes from a set of three tested candidate genes (ACTB, YWHAZ and CYCA) previously reported to be stably expressed in the brain (Bonefeld et al., 2008) to calculate a normalization factor for each sample. This normalization factor was then used to obtain the relative differences between the samples for each primer pair.

RNA SEQUENCING (RNA-SEQ)
A subgroup of samples (n=13; 5 PC, 8 TDT animals and 6 AC) was used for RNA-Seq, four pools of total RNA (with equal input from each individual sample) were made after Trizol extraction, and were further purified and DNase-treated using Qiagen colums (RNeasy Plus Microkit, Quiagen). The yield of the purified RNA ranged between 1.5 and 2 ug total RNA per pool. The five pools were as follows: (1)  that the animals needed to reach criterion in the behavioural experiment: LE: 9-10 days, ME: 6-7 days, SE: 5 days.
The purified RNA pools were sent for further quality control and RNA-Seq analysis to the Genomic Services Lab of the HudsonAlpha Institute for Biotechnology (AL, USA); http://gsl.hudsonalpha.org/. RNA-Seq with ribosomal RNA (rRNA) reduction was used using standard protocols (depth >45M pair-end reads per sample), and the resulted raw data were received. Data analyses (alignment and statistical analysis) were performed with GeneSifter (Geospiza, PerkinElmer Inc). We compared the control group with all three TDT groups pools (LS, LM, LF). For each comparison a p value (likelihood ratio test) and fold change (FC) was obtained. To obtain the list of differentially expressed (DE) genes in the control versus TDT groups, the following cut-offs were applied: p value ≤5E-5 in all 3 comparisons, FC>|1.25| in at least 2 out of 3 comparisons (with a minimum FC>|1.2|), and expression levels (reads per kilobase million, RPKM)>6 in at least 1 out of the 4 groups compared.

MRI STATISITICAL ANALYSIS
The data were pre-processed according to standard procedures in FSL (Smith et al., 2004). Tract Based Spatial Statistics (TBSS) (Smith et al., 2006) was applied to the preprocessed data. Images were then analysed as described elsewhere (Sampaio-Baptista et al., 2013). Briefly, all FA maps were aligned with linear and non-linear transformations to the study specific template and averaged to generate the mean FA image, from which the WM skeleton was extracted. The skeleton was thresholded at an FA value of 0.36 to contain only the major tracts (Sampaio-Baptista et al., 2013). Finally, the FA values of the tract centres were projected onto the skeleton for each rat brain and fed into statistical analysis. MD, RD and PD skeleton maps were created with the same method, using the FA registrations and skeleton projections as implemented in TBSS for non-FA images (Smith et al., 2006).
We used Permutation Analysis of Linear Models (PALM) (Winkler et al., 2016) for multi-measures analysis. PALM is a tool that allows inference over multiple modalities, including non-imaging data, using non-parametric permutation methods, similarly to the randomise tool in FSL (Winkler et al., 2014), although offering a number of features not available in other analysis software, such as the ability for joint inference over multiple modalities, or multiple contrasts, or both together, while correcting FWER or FDR across modalities and contrasts (Winkler et al., 2016).
We used PALM to assess the joint and individual contribution of the 4 DTI measures while simultaneously correcting across the tests. Non-Parametric Combination (NPC), as implemented in PALM, was used for joint inference over the 4 DTI measures (FA, MD, RD and PD). NPC works by combining test statistics or p-values of separate (even if not independent) analyses into a single, joint statistic, the significance of which is assessed through synchronized permutations for each of the separate tests. The synchronized permutations for the separate tests accommodate, implicitly, any eventual lack of independence among them. Such a joint analysis can be interpreted as a more powerful, permutation-based version of the classical multivariate analysis of covariance (MANCOVA); differently than MANCOVA, however, NPC allows investigation of the direction of joint effects.
Here we used NPC with Fisher's combining function, testing for effects with concordant directions across the 4 DTI measures. A cluster-forming threshold of t > 1.7 and 5000 permutations were used to determine p-values FWER-corrected for multiple comparisons (across all voxels and the 4 DTI measures). The chosen cluster-forming t threshold was based on the degrees of freedom of the sample. Clusters with a corrected significance of p < 0.05 were deemed significant.
We performed two statistical tests in the WM analysis. First we tested for differences between groups and included the total number of exposure days per animal as a covariate.
Second, we tested for correlations between performance rate and the 4 DTI measures. This was calculated by fitting a logarithmic model and extracting the slope of the percentage of correct trials curve for each individual animal (curves are illustrated in Fig. 1A).
Further, we tested for GM differences using MD only as this measure can indicate changes in tissue density regardless of the structure orientation. We tested for group differences and included the total number of exposure days per animal as a covariate. We performed non-parametric permutation testing with the Randomise tool as implemented in FSL (Smith et al., 2004), with a cluster-forming threshold of t > 1.7 and 5000 permutations were used to determine corrected p-values. Clusters with a corrected significance of p < 0.05 were deemed significant.

STATISTICAL ANALYSIS OF mRNA EXPRESSION
Statistical analysis of WM qPCR data was also performed with PALM (Winkler et al., 2016). We tested for group differences with non-parametric permutation testing with a between groups contrast (with total number of exposure days as a covariate). A p-value of < 0.05 was deemed significant, corrected for multiple comparisons across the 12 genes of interest.
For the qPCR analysis of the c-FOS gene of the Barrel cortex, a one-way analysis of covariance was used to test for differences between groups (with total number of exposure days as a covariate) with SPSS. A p-value of < 0.05 was deemed significant.

QPCR VALIDATION
The control versus TDT comparison led to the identification of 134 differentially- Gene Ontology (GO) enrichment analysis of the differentially expressed genes was performed using the web-based gene ontology tool from the Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.7 (http://david.ncifcrf.gov) (Dennis et al., 2003;Huang da et al., 2009a, b). For the enrichment analysis (Functional Annotation Chart tool), default software settings were used, and GO terms with a corrected p-value ≤ 0.01 (Benjamini correction) and at least 5 genes represented in the GO term were considered to be overrepresented (enriched).
Ingenuity Pathways Analysis (IPA) (Ingenuity Systems Inc., USA), was used to perform pathway, network and upstream regulator analyses to explore relationships between genes on the basis of curated information present in the IPA database. For pathway and interaction network analyses, a score was obtained (calculated as the -log of the associated Fisher's exact test p value). This score indicates the likelihood that the assembly of a set of focus genes in a network could be explained by random chance alone; networks with scores of 2 or higher have at least a 99% confidence of not being generated by random chance alone.
Upstream regulator analysis generated a list of putative upstream regulators of the DE genes, and indicate, for each putative upstream regulator, a predicted activation state, activation zscore, p value of overlap and list of putative target genes of the DE dataset.

RESULTS
Three months old rats (n=28) were trained on a TDT (von Heimendahl et al., 2007) that required a texture identity (rough or smooth) to be associated with a reward side (e.g., turn left on rough, right on smooth). After the initial shaping period, it took the trained animals between 5 and 17 days to reach criterion performance of 2 sessions > 80% accuracy ( Fig.1A, C). A subgroup of rats (n = 8) were further trained to detect increasingly more finegrained textures. After the rats had associated the correct reward side with the first texture, increasing the difficulty of the texture discrimination did not alter their accuracy (Fig. 1B, C). We first tested for between-group differences across all voxels in the WM skeleton. A significant cluster was found (p < 0.05, fully corrected across all voxels, Fig. 2A) covering a large area of WM under prefrontal and sensorimotor regions ( Fig. 2A). These are relevant WM regions for the texture detection task, in the light of the barrel cortex's role in somatosensory information processing (for review see (Feldmeyer et al., 2013)). The NPC partial tests showed that, compared to controls, the TDT group had significantly higher FA (p < 0.05, fully corrected across all voxels and the 4 measures) (Fig. 2B), lower MD (p < 0.05, fully corrected across all voxels and the 4 measures) (Fig. 2C), and lower RD (p < 0.05, fully corrected across all voxels and the 4 measures) (Fig. 2D) across similar areas of WM. PD was not significantly different between groups (p = 0.842, fully corrected across all voxels and the 4 measures).
Secondly, to assess the relationship between task performance and the neuroimaging structural measures, we again used a NPC Fisher's joint inference analysis to test for voxelwise correlations between 4 DTI measures and performance rate (slope of the individual curves) across individual animals (n = 28). We did not find a significant correlation between performance rate and the joint 4 DTI measures. However, there was a significant correlation between performance rate and FA (p < 0.05, fully corrected across all voxels and the 4 measures, Fig. 2E) and a trend for RD (p = 0.09, fully corrected across all voxels and the 4 measures, not shown) in a similar area of WM (90% overlap) to that showing group differences (depicted in Fig. 2A). Animals with higher FA (and lower RD) tended to show steeper slopes (i.e. they reached criterion performance with fewer exposure days), suggesting that WM microstructure is related to TDT performance. This correlation could either reflect experience-dependent changes that occurred with task performance, or pre-existing structural differences that give rise to performance variation, or a combination of both.

Grey Matter (GM) MD Is Lower In The TDT Group
To assess effects of experience on GM microstructure we tested for MD differences, as this measure reflects water restriction regardless of the structure orientation and can potentially indicate changes in tissue density and/or water content. This analysis revealed clusters with significantly lower MD (p < 0.05, corrected) in the TDT group (n=28) compared to the PC group (n=20) (Fig. 3). The significant clusters included both frontal and sensory cortex, hippocampus, and subcortical structures such as striatum and thalamic nuclei.

Learning Versus Experience
Our main analysis above compared the TDT group to a passive control group in order to identify general experience-dependent effects. To test for specific effects of associative learning versus experience in WM microstructure, we compared a subgroup (sg) of TDT rats (TDTsg, n = 12) to an active control group (AC) (n = 12). The AC group was matched for number of training days to a subgroup of TDT rats and were exposed to the same texture discrimination apparatus but were provided with rewards that were not contingent on their response. This allowed the AC group to experience the same textures and similar levels of rewards but without the requirement to differentiate between rough and smooth textures in order to gain the rewards. Rats in TDTsg and AC groups spent the exact same time in the task (mean = 8.83 days, S.D. = 3.97).
We used the same statistical approach as above and performed a non-parametric combination (NPC) for joint inference analysis as implemented by PALM, over the 4 DTI measures across all voxels in the WM skeleton, using the Fisher's combining function (Winkler et al., 2016). We tested for differences between the TDTsg (n=12), AC (n = 12) and PC (n=20) groups.
The Fisher test did not reveal significant differences between these groups when considering all measures together. However, the partial tests revealed several trends for FA.
A trend was seen for higher FA in the AC group compared to the PC group in a cluster underlying barrel cortex (p = 0.1, fully corrected, Fig. 4A, C) and in a very similar cluster for the TDTsg group compared to the PC group (p = 0.09, fully corrected, Fig. 4B, C). There were no significant differences nor trends between the TDTsg and AC groups (p = 0.71, fully corrected). These results suggest that learning to distinguish between smooth versus rough textures is not necessary for the detected structural differences and that exposure to the task apparatus, whisker stimulation, and rewards, is sufficient to elicit similar structural changes in both active groups compared to a caged group. This indicates that the requirement to differentiate between rough and smooth textures is not necessary to elicit structural WM changes.

Synaptic C-Fos mRNA Expression In The Barrel Cortex Is Higher In The TDT Group
We assessed synaptic C-FOS expression as an indirect marker of cell activity in the barrel cortex to confirm activation of this area in response to the task. As expected, c-FOS mRNA expression was found to be higher in the TDT group (n = 16) compared to the PC (n =10, p < 0.01) (Fig. 5A).

Measures Of WM Microstructure
We further tested for correlations between c-FOS mRNA and the mean FA, MD and RD values of the significant cluster identified in Fig. 2A. We used Bonferroni correction accepting a p-value smaller than 0.0083 as significant.
There were significant negative correlations (Fig. 5B

Myelin basic protein (MBP) mRNA Expression Is Higher In The WM Co-Localised With The Barrel Cortex Of The TDT Group
To assess effects of experience on mRNA expression in the WM, the WM underlying the barrel cortex was dissected from a subset of right-brain hemispheres (n=24; 15 from the TDT group and 9 from the control group) and qPCR was performed on this tissue. correction. MBP mRNA expression was found to be significantly higher in the TDT group (p < 0.05, corrected for multiple comparisons, Fig. 5D).
We tested for correlations between mRNA and the mean FA, MD and RD values of the cluster found to be significant in Fig. 2A with Pearson correlation coefficient. No significant correlations were found between the mRNA expression of the candidate list and the DTI measures.

Genome-wide RNA sequencing
To gain further insight into the molecular mechanisms underlying the observed structural WM and candidate gene expression differences, we performed an unbiased genome-wide analysis of mRNA expression in WM underlying the barrel cortex, by means of RNA-sequencing in a subgroup of samples (n=19; including 5 passive controls, 8 TDT animals and 6 active controls).
The TDT versus PC comparison led to the identification of 134 differentially expressed (DE) genes, of which 65 were up and 69 downregulated (likelihood ratio test, p ≤ 5E-6, FC cut-off |1.25|, RPKM cut-off 6) (Fig. 6A). From this list of 134 DE genes, 124 genes were also differentially expressed between PC and AC, in the same direction (up or down) and with similar significance as in the PC versus TDT comparison. The TDT versus AC comparison lead to only 6 DE genes (Notch3, Tns1, Zbtb16, Nxn, Yap1, Cfap43), all of them downregulated.

Gene ontology and Ingenuity Pathway analysis
In order to interpret the biological significance of the differentially expressed genes, gene ontology (GO) analysis and Ingenuity Pathways Analysis (IPA) were performed. The list of 6 DE genes (TDT versus AC comparison), led to no significant findings at the designated threshold (see Methods). In contrast, the list of 134 DE genes (TDT versus PC comparison) yielded several significant findings that are reported in detail below.  Table 1). This indicates that somatosensory experience induces activation of the transcriptional machinery in WM. For example, MAPK signalling pathway has been implicated in cell proliferation, differentiation and development (Zhang and Liu, 2002). genes of this list were also significantly differentially expressed between PC and AC. Warm

IPA analysis: Upstream Regulators and Networks
To identify molecules upstream of the genes that potentially explain the observed 134 DE genes, an IPA 'Upstream Regulators' analysis was performed (see Methods) (Supplementary Table 2). In simple terms, this analysis uses the IPA database to identify upstream regulators that match the direction of regulation of the downstream differentiallyexpressed molecules from the RNA sequencing dataset (Fig. 6A).
This revealed five predicted upstream regulators showing a high degree of concordance between predicted (by the IPA database) and actual direction of regulation (CREB1, CREM, GnRH-A, dalfampridine and bicuculline) (relationships illustrated in Fig.   6B), all of them with a predicted activated state, an activation z-score > 2.5 and a p-value < 5E-15, with at least 12 target molecules of the DE list dataset.
CREB1 and CREM had overlapping target molecules (29 target molecules for CREB1, of which 16 were also CREM targets), and had a high degree of overlap with GnRH-A targets (12 targets, of which 9 overlapping with CREB1/CREM targets) (Fig. 6B bottom). CREB and CREM are transcription factors activated by phosphorylation in response to cAMP and other signals (for review see (Mayr and Montminy, 2001)).
Additionally, two drugs, dalfampridine and bicuculline, shared the same 16 target molecules (Fig. 6B top). Dalfampridine is a broad-spectrum voltage-gated potassium channel blocker that broadens the action potential. It shows beneficial effects in multiple sclerosis patients, probably through restoration of axonal conduction (Dunn and Blight, 2011) and has been shown to promote remyelination after acute nerve injury (Tseng et al., 2016).
Bicuculline blocks GABA A -mediated inhibition, thereby increasing neuronal activity. This shows that genes involved in neuronal activity modulation are differentially expressed in our sample, possibly indicating increases or functional changes in axonal activity.
Next we sought to identify functionally related networks of genes and important regulatory hubs (Fig. 7). There were two networks with a score >30 containing at least 17 molecules of the DE dataset. In Network 1, two main 'hub' molecules (i.e. where most relationships converge to), Akt and Creb (Fig. 7A), were visually identified, while in Network 2 the main 'hub' molecules were Erk1/2 (Fig. 7B)

DISCUSSION (1013 words)
Somatosensory experience results in structural changes including increases in dendritic spines  and cortical myelinogenesis (Hughes et al., 2018). Here we have demonstrated that somatosensory experience also results in structural WM plasticity and identify some molecular correlates that provide candidate mechanisms underlying these findings. We did not find structural or genome-wide mRNA expression differences between the TDT group and an active control group, suggesting that learning the associative task is not necessary for the detected plastic changes and that mere exposure to somatosensory stimulation is sufficient.
The WM structural diffusion metrics were found to correlate with barrel cortex synaptic c-fos expression, suggesting that molecular correlates of cortical activity relate to macroscale measures of WM structural plasticity. Synaptoneurosomes are enriched in synaptic terminals (pre and post) and might also include other cellular (for instance astrocytes) components. Although c-fos is often regarded as an immediate early gene with exclusive expression in neurons, its expression has also been shown in astrocytes under certain conditions, so we cannot completely exclude non-neuronal contribution to the synaptic c-fos expression (Herrera and Robertson, 1996). There is increasing evidence from in vitro (Demerens et al., 1996) and in vivo (Mensch et al., 2015;Etxeberria et al., 2016) studies that neuronal activity modulates myelination, even in adulthood (Gibson et al., 2014).
In line with this, our mRNA expression analysis indicated that myelination-related genes (e.g. MBP, Akt and Erk1/Erk2) (Ishii et al., 2012;Jeffries et al., 2016) are differentially expressed in response to somatosensory experience. Further, the DTI analysis revealed higher FA and lower RD in the TDT group indicating that water diffusion is more hindered across WM tracts after somatosensory experience. Additionally, lower MD was found diffusion in both WM and GM, indicating higher overall restriction of water which could potentially be related to greater tissue density in these areas. While definitive biological interpretation of DTI changes is challenging (Sampaio-Baptista and Johansen-Berg, 2017), this pattern of DTI differences in WM is consistent with cellular mechanisms such as higher myelin thickness or internode length (Etxeberria et al., 2016). Accordingly, we found higher MBP mRNA expression suggesting that myelination has been triggered by somatosensory experience.
Increases in myelin are consistent with the higher FA and decreased RD found in the current study. These findings are congruent with previous neuroimaging studies that have found higher FA, along with higher MBP immunostaining intensity in response to complex motor and cognitive learning (Blumenfeld-Katzir et al., 2011;Sampaio-Baptista et al., 2013).
However, additional mechanisms may also contribute. For example, changes in axon diameter (Sinclair et al., 2017) We have also detected expression differences in a set of genes involved in neuronal modulation and another set of genes involved in regulation of myelin formation or remodelling, which together could underlie, at least in part, the structural WM findings detected here via DTI. The genome-wide mRNA analysis identified 134 differentially expressed genes that were found to be or part of gene networks or regulated by molecules, involved in functions related to neuronal plasticity (CREB1, CREM), memory (CREB1, CREM), neuronal activity modulation (Dalfampridine and Bicuculline), myelin-sheath thickness regulation (Erk1/Erk2), and proliferation control, differentiation and protein synthesis in OPCs (CREB1, CREM, Akt and Erk1/Erk2). All these indicate that WM has undergone functional and structural plasticity in response to somatosensory experience. In particular, CREB and CREM are involved in regulating the transcription of several genes (cfos, BDNF, etc) and have been implicated in neuronal plasticity and memory (Lonze and Ginty, 2002;Benito and Barco, 2010). More recently, CREB has been connected with the transcriptional control of MBP (Meffre et al., 2015), also found to be differentially expressed in our study. Several studies strongly support the role of Akt and mTOR signalling in promoting myelination (reviewed in (Norrmen and Suter, 2013)). Erk1/Erk2, together with Akt/mTOR, have been proposed as important signalling pathways for the control of proliferation, differentiation and protein synthesis in oligodendrocyte precursor cells (OPCs) (Cui and Almazan, 2007;Bibollet-Bahena and Almazan, 2009;Dai et al., 2014). Erk1/Erk2, from the mitogen-activated protein kinase (MAPK) pathway (also significantly enriched in our GO analysis), is an important regulator of myelin-sheath thickness in the CNS (Ishii et al., 2012;Jeffries et al., 2016). Conditional upregulation of Erk1/Erk2 results in global increases in myelin thickness by preexisting oligodendrocytes of adult mice, faster nerve conduction velocity and behavioural changes (Jeffries et al., 2016). Furthermore, Erk2 has been described to have an important role in oligodendrocytes in the translational control of MBP (Michel et al., 2015), which is in line with our findings. These results suggest that somatosensory experience may potentially trigger both de novo myelin formation as well as thickness increases of myelin sheaths.
In conclusion, we report structural WM and GM differences and identify genes that are differentially expressed after somatosensory experience. Additionally, macroscale WM structural plasticity was found to relate to cortical activity as measured by c-fos mRNA, suggesting that cortical experience-dependent mechanisms trigger WM changes. Future studies could also examine the specific effects of the identified genes on MRI measures by combining genetic (McKenzie et al., 2014;Jeffries et al., 2016) or pharmacological manipulation in rodents with imaging read-outs. This would allow to precisely identify the molecular and cellular mechanisms which underlie changes in MRI measures of plasticity and could offer important clues to the biological changes underlying imaging signals recorded in humans.