Transcriptional mapping of the primary somatosensory cortex upon sensory deprivation.

Experience-dependent plasticity (EDP) is essential for anatomical and functional maturation of sensory circuits during development. Although the principal synaptic and circuit mechanisms of EDP are increasingly well studied experimentally and computationally, its molecular mechanisms remain largely elusive. EDP can be readily studied in the rodent barrel cortex, where each "barrel column" preferentially represents deflections of its own principal whisker. Depriving select whiskers while sparing their neighbours introduces competition between barrel columns, ultimately leading to weakening of intracortical, translaminar (i.e., cortical layer (L)4-to-L2/3) feed-forward excitatory projections in the deprived columns. The same synapses are potentiated in the neighbouring spared columns. These experience-dependent alterations of synaptic strength are thought to underlie somatosensory map plasticity. We used RNA sequencing in this model system to uncover cortical-column and -layer specific changes on the transcriptome level that are induced by altered sensory experience. Column- and layer-specific barrel cortical tissues were collected from juvenile mice with all whiskers intact and mice that received 11-12 days of long whisker (C-row) deprivation before high-quality RNA was purified and sequenced. The current dataset entails an average of 50 million paired-end reads per sample, 75 base pairs in length. On average, 90.15% of reads could be uniquely mapped to the mm10 reference mouse genome. The current data reveal the transcriptional changes in gene expression in the barrel cortex upon altered sensory experience in juvenile mice and will help to molecularly map the mechanisms of cortical plasticity.

Data Description

Context

Sensory experience powerfully shapes neural circuits. Changes due to sensory organ deprivation such as eye closure, digit amputation, and whisker trimming provide powerful means for studying mechanisms of experience-dependent cortical plasticity.

In the whisker system, experience-dependent plasticity is most commonly studied in the barrel cortex subfield of the primary somatosensory cortex where neural representations of whiskers change in response to altered patterns of incoming sensory information. As originally shown in the barrel cortex [1], sensory deprivation induced by transient whisker trimming is sufficient to perturb neural receptive fields both during development and in adulthood. Previous work has also shown that the cellular basis of deprivation-induced decreases in whisker-evoked representations are primarily due to a reduction of synaptic strength in monosynaptically connected feed-forward neuronal networks in behaving animals [2, 3]. Conversely, whisker-sparing-induced enhancement in whisker representation is mediated at least in part by the long-term synaptic facilitation expressed along the L4 projections in vivo [4]. Identification of the molecular events that mediate these bidirectional changes in synaptic connectivity will benefit from systematic analysis of the gene transcription. Therefore, we performed RNA sequencing in the barrel cortex with or without sensory deprivation across cortical layers 2–4. This database will assist molecular and cellular neurobiologists in addressing the molecular mechanisms associated with experience-dependent plasticity and will enable statistical approaches to determine the dynamics of the coupled changes across molecular pathways as cortical circuits undergo plastic changes in their organization.

Methods

Animals

All experiments were performed in accordance with the Animal Ethics Committee of the Radboud University in Nijmegen, the Netherlands. Pregnant wild-type mice (Charles River: Wilmington, Massachusetts, United States, stock number 000664; RRID:NCBITaxon_10090) were kept at a 12-hour light/dark cycle with access to food ad libitum. Cages were checked for birth daily. To induce experience-dependent plasticity, pups underwent bilateral plucking of their C-row whiskers under isoflurane anaesthesia at P12 (Fig. 1). Control animals were not plucked but anaesthetized and handled similarly. After recovery, pups were returned to their home cage. Every other day, pups were checked for whisker regrowth, and whiskers were plucked if present. At P23–P24, pups were randomly selected from their litter for slice preparation and tissue collection. For each experimental condition (i.e., whisker deprived or control), 4 female pups were used; thus each group consisted of 4 independent biological samples (also known as biological replicates). Samples from cortical layer (L) 4 and L2/3 were treated independently with their own corresponding groups of control, deprived, first-order spared, and second-order spared columns, as detailed in Fig. 1.

Figure 1:

Overview of the experimental design, sample collection, and data organization. (A) Pups were bilaterally spared or deprived of their C-row whiskers between P12 and P23–P24, when acute slices were made and column- and layer-specific tissues were excised. (B) RNA was isolated, checked for integrity and purity, and subsequently sequenced. (C) Organization of the database. Colour codes denote experimental groups. Same denominations are used in the read counts matrix file (see the Supplementary Data).

Slice preparation and sample collection

Pups were anaesthetized using isoflurane and then perfused with ice-cold carbogenated slicing medium (108 mM ChCl, 3 mM KCl, 26 mM NaHCO3, 1.25 mM NaH2PO4, 25 mM glucose, 1 mM CalCl2, 6 mM MgSO4, and 3 mM Na-pyruvate). Next, pups were decapitated and the brain was quickly dissected out, and 400 μm thalamocortical slices from each hemisphere were prepared as described before [2, 3]. Slices were transferred to 37°C carbogenated artificial cerebrospinal fluid (ACSF) (120 mM NaCl, 3.5 mM KCl, 10 mM glucose, 2.5 mM CaCl2, 1.3 mM MgSO4, 25 mM NaHCO3, and 1.25 mM NaH2PO4) where they were kept for 30 minutes and recovered at room temperature for another 30 minutes until tissue collection.

After incubation, slices were placed under a Nikon: Eclipse FN1 microscope (Nikon: Minato, Tokyo, Japan). The holding chamber was continuously perfused with room-temperature carbogenated ACSF. Due to the 55° cut, slices were obtained in which S1 barrels from specific rows (A–E) could be identified [2]. A thin, long glass pipette was pulled using a Sutter instruments P-2000 pipette puller, which was used to make intercolumnar incisions from L1 to the bottom of L4, after which the slice was placed under a binocular dissection microscope, where the location of specific barrel columns could now be readily identified by eye. A sterile 32G needle was then used to cut out L2/3 and L4 separately from each column. Tissue from columns A/E and B/D were pooled as they both constitute second- and first-order spared whiskers, respectively. Immediately after dissection, tissue samples were snap-frozen in liquid nitrogen and stored at −80°C until further use. All tools that came into direct contact with brain tissue were treated using RNAseZap (Thermo Fisher Scientific: Waltham, Massachusetts, United States, #AM9780) in order to minimize RNAse contamination.

RNA solation and quality control

Tissue samples originating from the same rows and layers were pooled within each animal. From control animals, only the C column tissues were used (also see the Re-use potential section). Tissue was quickly dissolved in Qiazol (Qiagen: Hilden, Germany, #79306), after which RNA was isolated using the miRNeasy Mini kit (Qiagen: Hilden, Germany, #217004), DNAse treated (Thermo Fisher Scientific: Waltham, Massachusetts, United States, #EN0521), and cleaned up using the RNeasy MinElute kit (Qiagen: Hilden, Germany, #74204), all following the manufacturer's instructions. Samples were then stored at –80°C until further processing.

RNA sample integrity was determined using Agilent Tapestation (High Sensitivity RNA Screentape). Sample RINs ranged from 7.1 to 8.8. To further assess RNA purity and integrity, RNA samples were used in reverse transcription polymerase chain reaction (RT-PCR) to confirm that cDNA could be produced and that a large (∼1000 bp) amplicon could be obtained. To produce cDNA, SuperScript II Reverse Transcriptase (Thermo Fisher Scientific: Waltham, Massachusetts, United States, #18064014) and random hexamer primers (Roche: Basel, Switzerland, #11034731001) were used. The resulting cDNA was then added to a PCR reaction mix, which further consisted of Jumpstart Ready Mix (Sigma P2893) and exon-exon junction-spanning CamKII primers (FW TCCAACATTGTACGCCTCCAT; RV TGTTGGTGCTGTCGGAAGAT). From all cDNA samples, a fragment of the expected size could be amplified, suggesting that the RNA samples contained pure RNA of sufficient integrity. All RNA samples thus passed our quality control criteria and were subjected to RNA sequencing.

RNA sequencing

RNA sequencing was conducted at the Genomics Core Facility of the EMBL, Heidelberg, Germany (RRID:SCR_004473). The cDNA library was generated using the non-stranded NEBNext Ultra RNA Library Preparation Kit for Illumina (NEB: Ipswich, Massachusetts, United States, catalogue #E7530), which includes oligo-dT bead selection of mRNA. For library enrichment, 13–14 PCR cycles were performed. Pooled libraries were sequenced on the Illumina: NextSeq 500 instrument (Illumina: San Diego, California, United States of America) (RRID:SCR_014983) in a 75-bp paired-end mode using high-output flow cells.

Data validation and quality control

Sequencing read quality was assessed using FastQC (Babraham Bioinformatics: Babraham, England; RRID:SCR_014583), the results of which were merged using MultiQC (RRID:SCR_014982) [5]. The results are displayed in Fig. 2. Per base quality phred scores range from 34.80 to 35.15, indicating base call accuracies of >99.9% (Fig. 2A). Overall, 91.48–94.03% of reads had a mean phred score of 30 or above (Fig. 2B). In line with these scores, per base N content (i.e., percentage of bases that could not be confidently called) was very low, with a maximum value of 0.053%.

Figure 2:

FastQC and STAR output graphs for all samples. (A–B) Phred scores per base and per sequence. (C) Per sequence GC content. (D) STAR output of alignment scores.

Reads were then mapped to the mm10 reference genome using STAR (RRID:SCR_005622) [6], which uniquely mapped between 39 000 000 and 59 000 000 reads, constituting an average 90.15% unique map rate across samples (Fig. 2D). Since the library preparation protocol entails a PCR enrichment step, which can lead to technical duplication and hence an overestimation of observed transcripts, we used Seqmonk (Babraham Bioinformatics: Babraham, England; RRID:SCR_001913) to plot the read density against the duplication levels (i.e., the percentage of duplicate reads) for each transcript. The obtained duplication plots showed a clear positive relation between read density and duplication levels (Fig. 3; Supplementary Fig. S1), suggesting that the origin of read duplication is biological, rather than technical.

Figure 3:

Overlays of duplication plot contours, showing a positive correlation between read density and duplication levels. Depicted contours enclose 90% of the data points.

Based on the above quality control measures, we determined that our RNA-sequencing data was of sufficient quality to be used in downstream analyses; therefore we continued with gene expression analysis.

Analysis of gene expression

Using a 2-read cut-off, we identified 16 900 to 17 600 transcripts per sample (Fig. 4A). Raw gene counts can be found online (see the Supporting Data for [7]). Differential gene expression analyses across groups were performed using EdgeR v. 3.12.1 (RRID:SCR_012802) [8, 9] using only genes with a count per million (CPM) >1 in at least 4 samples (Supplementary Table S1 for details on the commands used). Since laminar identity is an important feature of our experimental setup, we assessed the relative expression of known molecular markers for L2/3 (Cacna1h, Id2, Igfbp4, Igfn1, Mdga1, Plcxd1, Rasgrf2, Rgs8, Tle3) and L4 (Cartpt, Cyp39a1, Kcnh5, Kcnip2, Lmo3, Rorb, Scnn1a) [10-12], which showed selective enrichment of the laminar markers in isolated layers (Fig. 4B).

Figure 4:

Gene expression analyses. (A) Histogram of read counts per transcript per sample. With a cut-off of 2 reads, between 16 900 and 17 600 transcripts could be identified across samples. (B) Relative expression of known molecular markers for cortical laminae. Layer 4 markers are enriched in samples originating from this layer; the same is true for layer 2/3 marker expression in layer 2/3 samples. (C) Cumulative plots of the CV of individual experimental groups. Including only transcripts identified by 50 reads or more, average CVs of <15% are found in ∼85% of transcripts. (D) PCA showing sample clustering by layer, including only transcripts identified by at least 50 reads. PC1 and 2 account for 88% of overall variance.

To assess the variance in transcript counts, we calculated the coefficient of variation (CV) for each transcript with a cut-off of 50 as the minimal read count separately for each group (Fig. 4C). This analysis showed that, on average, 85.93% of transcripts have a CV below 15%, suggesting low variance across transcript counts for individual genes. Principal component analysis (PCA) showed that samples cluster based on layer, and the first 2 components explained ∼88% of the variance in the data (Fig. 4C; Supplementary Fig. S2B).

These quality control routines suggest that we have obtained RNA-sequencing data of high read quality, with individual bases being called confidently throughout the length of reads, which uniquely map to the mm10 reference genome at high rates (>90% average). The laminar origin of our samples could be identified through known molecular markers, confirming our samples are of high anatomical specificity.

Re-use potential

The current RNA-seq dataset might help address the molecular underpinnings of cortical experience-dependent plasticity. For example, it could be used (i) to identify genes whose transcription is modulated in an experience-dependent manner, (ii) to statistically map the transcriptional networks at laminar resolution, (iii) creating synergy with the single neuron RNA-seq datasets [13, 14], to address the molecular diversity of the cortical networks, (iv) combined with the proteomic analysis performed under comparable experimental conditions in the accompanying manuscript (Kole et al., submitted), to systematically study the transcriptional and translational regulation of the genome upon altered sensory experience, and finally (v) to identify and quantify splice isoforms given the sequencing depth of the current dataset. Since splicing and other posttranscriptional mechanisms govern which proteins are ultimately produced, combining the current transcriptomic dataset with a proteomics approach [15] would also be of high importance.

The current dataset focuses on isolated cortical columns and layers, which are necessarily diverse samples containing neuronal and non-neuronal cell classes. In terms of experience-dependent plasticity, although most previous studies focus on excitatory projections, inhibitory cells and even non-neuronal cells have been implicated in plasticity [16-18]. This heterogeneity might be particularly important for L2/3, as also shown by the principal component analysis (Fig. 4D), given the relative diversity of cellular populations in supragranular layers and their heterogeneous connectivity patterns [19].

Researchers reusing our dataset should be aware that comparisons between control column C and spared columns (A/E, B/D) may have to be approached with caution as this would involve 2 different columnar identities (whose transcriptomic dissimilarities are currently unknown), each coming from cortices that have had different sensory experience. However, direct comparisons between the C columns across experimental conditions (i.e., control vs deprived) as well as within-animal across-column comparisons in deprived animals control for these confounding variables.

Taken together, we hope that this data will prove useful in discovering the novel molecular targets responsible for cortical plasticity and will lead to targeted control of plasticity in health and disease.

Availability of the supporting data

All supporting data are available in the GigaScience repository, GigaDB [7].

The raw sequence reads were deposited in the NCBI under GEO accession GSE90929.

Additional files

Supplementary Figure S1. Duplication plots for all samples, produced using SeqMonk (Babraham Bioinformatics: Babraham, England).

Supplementary Figure S2. (A) Cumulative plots of the CVs of each experimental group, including transcripts identified by at least 1 read. Average CVs of <25% are found in ∼85% of transcripts. (B) PCA including transcripts identified by at least 1 read. The majority (88%) of overall variance is explained by principal components (PC) 1 and 2.

Abbreviations

EDP: experience dependent plasticity; L2/3: cortical layer 2/3, also known as supragranular layers; L4: cortical layer 4, i.e., granular layer.

Competing interests

The authors declare that they have no competing interests.

Funding

Funding for the current work was provided by the Faculty of Science of the Radboud University, Nijmegen, the Netherlands (grant number 626830–6200821) as well as the ALW Open Programme of the Netherlands Organization for Scientific Research (grant number 824.14.022).

Author contributions

K.K. performed all experimental manipulations, sample acquisition, biological and bioinformatic quality controls, and prepared the tables and figures. Y.K. and Ja.P. performed bioinformatic analysis. Je.P. performed library prep. V.B. supervised RNA-seq. P.T. contributed bioinformatic analysis and co-supervised the project. T.C. designed and supervised the project. K.K. and T.C. wrote the manuscript. All authors edited otherwise approved the final version of the manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.