Analysis of hippocampal transcriptomic responses to technical and biological perturbations

Cost-effective next-generation sequencing has made unbiased gene expression analysis possible. Single-neuron gene expression studies may be especially important for understanding nervous system structure and function because of the neuron-specific functionality and plasticity that defines functional neural circuits. Cellular dissociation is a prerequisite technical manipulation for single-cell and single cell-population studies, but the extent to which the cellular dissociation process cells affects neural gene expression has not been determined, nor has it been determined how gene expression is altered by the stress that accompanies many of the behavioral manipulations that are required to study learning and memory and other cognitive functions. Here, we determined to which extent cellular dissociation-induced changes in hippocampal gene expression might confound studies on the behavioral and physiological functions of the hippocampus. We processed tissue punch samples from the dentate gyrus (DG), CA3, and CA1 hippocampus subfields using either a tissue homogenization protocol or a cellular dissociation protocol in preparation for RNA sequencing analysis to evaluate the impact of the tissue preparation. Then, we evaluated the effect of stressful experience and cognitive training on hippocampus subfield specific gene expression and determined to which extent these response overlap with the cellular dissociation response. Finally, we assessed the extent to which the subfield-specific gene expression patterns are consistent with those identified in a recently published hippocampus subfield-specific gene expression database. We report substantial differences in baseline subfield-specific gene expression, that 1% of the hippocampal transcriptome is altered by the process of cellular dissociation, that an even weaker alteration is detected 24 h after stressful experience, and that while these alterations are largely distinct from the subfield specific response of the hippocampus transcriptome to cognitive training, there is nonetheless some important confounding overlap. These findings of the concordant and discordant effects of technical and behavioral manipulations should inform the design of future neural transcriptome studies and thus facilitate a more comprehensive understanding of hippocampal function.

learning and memory and other cognitive functions. Here, we determined to which extent 23 cellular dissociation-induced changes in hippocampal gene expression might confound 24 studies on the behavioral and physiological functions of the hippocampus. We processed 25 tissue punch samples from the dentate gyrus (DG), CA3, and CA1 hippocampus subfields 26 using either a tissue homogenization protocol or a cellular dissociation protocol in 27 preparation for RNA sequencing analysis to evaluate the impact of the tissue preparation. 28 Then, we evaluated the effect of stressful experience and cognitive training on hippocampus 29 subfield specific gene expression and determined to which extent these response overlap 30 with the cellular dissociation response. Finally, we assessed the extent to which the 31 subfield-specific gene expression patterns are consistent with those identified in a recently 32 published hippocampus subfield-specific gene expression database. We report substantial 33 differences in baseline subfield-specific gene expression, that 1% of the hippocampal 34 transcriptome is altered by the process of cellular dissociation, that an even weaker 35 alteration is detected 24 h after stressful experience, and that while these alterations are   Hatfield, PA). All punches for RNA sequencing came from the slice corresponding to image 103 74 of the Allen Brain Reference Atlas (RRID:SCR_013286).

104
Animal and tissue preparation for assessing impact of cellular dissociation 105 A 1-year-old female C57BL/6J mouse was used for the cellular dissociation experiment. One 106 tissue punch was designated for the control homogenized processing and the other for the 107 cellular dissociation treatment (Fig 1A). Two adjacent tissue samples were collected from each 108 subfield for each mouse. The 'control sample' was processed using the manufacture instruc-  Male C57BL/6J mice that were 3-4-months old were used. They were obtained from the Jack-119 son Laboratory (Bar Harbor, ME) and housed at the Marine Biological Laboratory. Gene ex-120 pression in tissue from mice taken from the home cage was compared to mice that received 121 mild foot shock, to evaluate how gene expression is affected by stressful experience, which 122 is a common confound of behavioral manipulations such as water maze learning, fear con-    Twenty-four hours after the behavioral manipulations the mice were sacrificed and tissue 148 punches were collected was described above.    We used the VennDiagram R package (Schwenk et al., 1984) for preliminary visualization of 175 differential gene expression, but the final Venn diagrams were drawn with Adobe Illustrator.

176
The hierarchical clustering analysis was conducted and visualized using the R package  . The PCA analysis was visualized using the ggplot2, cowplot, and using -log(p-value) as a continuous measure of significance to identify GO categories that 185 are significantly enriched with either up-or down-regulated genes. No significance cutoff is 186 required for the analysis, but an arbitrary p-value was set to visualize the top 10 most significant 187 GO terms. Figure 5 was generated in two steps. First, a p-value = 0.1 was set for determining 188 significantly expressed genes in each analysis, and these data were converted into a binary (0 189 or 1) for a typical GO enrichment analysis using Fisher's exact test to determine if GO categories 190 are overrepresented among the significantly expressed genes. The effects of cellular dissociation on hippocampal transcriptomes 202 We identified 162 genes that were differentially expressed between the control and dissoci-203 ated samples, 331 genes that were differentially expressed genes (DEGs) between any of the 204 three hippocampus subfields, and 30 genes were shared between both sets of differentially 205 expressed genes at p-value < 0.05 (Fig. 1A,B).
206 Figure 1. The effect of cellular dissociation on hippocampal transcriptomes. A) From a single female mouse, we collected 2 CA1, CA3, and DG hippocampal tissue samples. One sample was subjected to a cellular dissociation treatment (dissociated) whereas the control samples (control) were standardly homogenized. B) We identified 162 dissociation-induced changes in gene expression, 331 genes with region-specific expression patterns, and 30 genes differentially expressed by both region and treatment (FDR p-value < 0.05). C) Hierarchical clustering separates the hippocampal subfields of the homogenized samples (light gray) but not the dissociated samples (dark gray). D) PC1 accounts for 40% of all gene expression variation and by inspection, separates the DG samples from the CA1 and CA3 samples. PC2 accounts for 22% of the variation in gene expression and varies significantly with treatment. Ellipses are hand-drawn.
A hierarchical clustering analysis of all differentially expressed genes does not give rise 207 to distinct clusters that are separated by subfield or method; however, when examining the 208 control, homogenized samples alone (identified with light grey boxes), the three subfields 209 form distinct clusters, while the dissociated samples do not cluster by subfield (Fig. 1C). 210 Next, we conducted a principal component analysis of all identified genes. PC1 accounts 211 for 40% of the variation and visually separates the DG samples from the CA1 and CA3 sam-212 ples (Fig. 1D). To confirm statistical significance of this visual pattern, we conducted a two- of the three brain regions at p-value < 0.05 ( Fig. 2A, B).
226 Figure 2. The effects of a stressful experience on hippocampal transcriptomes. A) We compared CA1, CA3, and DG tissue samples from control mice taken directly from their home cage to mice that were subjected to a mild foot shock. B) We identified 0 genes that responded to treatment, and 1669 genes that were differentially regulated across regions of the hippocampus (FDR p-value < 0.05). C) Hierarchical clusters groups samples by brain region but distinct treatments clusters are not present. D) PC1 accounts for 31% of the variation and visually separates the DG samples from the CA1 and CA3 samples. PC2 accounts for 1% of the variation and distinguish the three subfields. Ellipses were hand-drawn.
Hierarchical clustering of the differentially expressed genes gives rise to three distinct clus-227 ters corresponding to the three subfields, with CA1 (purple) and CA3 (green) being more similar 228 to one another than to DG (orange), whereas the effects of the stress manipulation were not 229 distinctive (Fig. 2C). 230 Next, we conducted a principal component analysis of all the genes that were measured.

231
PC1 accounts for 31% of the variation and by inspection, separates the DG samples from the 232 CA1 and CA3 samples (effect of region: F2,15= 42.89; p < 0.001; Fig. 2D). Post hoc Tukey tests 233 confirmed CA1 = CA3 ≠ DG. The effects of stress and the stress x region interaction were not 234 significant. PC2 accounts for 18% of the variation and varies significantly between CA1 and 235 CA3 and CA1 and DG (effect of region: F2,15= 11.41; p < 0.001; Tukey tests: CA1 ≠ DG = CA3). The 236 effects of stress and the stress x region interaction were not significant. PC3 accounts for 15% 237 of the variation and also explains some brain region specific differences (effect of region: F2,15= 238 6.315; p < 0.01; Tukey tests: CA1 = DG ≠ CA3), whereas effects of stress and the stress x region 239 interaction were not significant. PC4 is also influences by region (F2,15= 6.315; p = 0.0102; Tukey 240 tests: CA1 ≠ CA3. PC5 did not account for any significant differences according to region or 241 treatment. PC6 significantly accounted for variance associated with the effect of a stressful 242 experience (F1,16> 4.774; p's < 0.04).

243
The effects of cognitive training on hippocampal transcriptomes 244 We identified that 423 genes were differentially expressed between the yoked control and 245 cognitively trained animals, 3485 genes that were differentially expressed across subfields, and 246 324 showed an interaction at FDR p < 0.05 (Fig. 3A, B). We see a large effect of brain region 247 on gene expression, with 20% of detectable genes being differentially expressed between 248 one or more brain-region comparisons (3485 differentially expressed genes /17320 measured 249 genes). This is an order of magnitude greater than the 2% of the transcriptome that changed 250 in response to learning (423 DEGs /17320 genes measured). Figure 3. Effects of a learned avoidance behavior on hippocampal transcriptomes. A) Mice used in this study were either subjected to random but mild foot shocks (control) or subjected to mild foot shocks conditioned with spatial cues (trained). Tissue samples were collected from CA1, CA3, and DG. B) We identified only 285 genes that were significantly expressed according to the behavioral manipulation and identified 3622 genes that were were differentially expressed between any of the three brain regions. C) Hierarchical clustering of the differentially expressed genes gives rise to three distinct clusters corresponding to the three brain regions, with CA1 and CA3 being more similar to one another than to DG. D) A principle component analysis of all genes in analysis (regardless of level of significance) shows that PC1 accounts for 50% of the variation and distinguishes the DG samples and the CA1 and CA3 samples (Region: F2,19= 199.3; p = 1.78e-13). PC2 accounts for 18% of the variation and distinguishes all three subfields (F2,19= 220.4; p = 7.15e-14). Ellipses were hand-drawn.

251
Hierarchical clustering of the differentially expressed genes separates samples by both 252 subfield and treatment (Fig. 3C). A principal component analysis of all gene expression data 253 revealed that brain region explains 75% of the variance in the data (Fig. 3D) Next, we examined the overlap in genomic response to the technical and biological pertur-263 bations. We identified three specific genes that responded to both cellular dissociation and 264 to cognitive training: Grin2a, Epha6, and Ltbp3 (Fig. 4A). There was no overlap in differentially 265 expressed genes compared to the cellular dissociation treatment (Fig. 4A). 266 Figure 4. Unique and shared responses to technical treatments and biological perturbations. A) The number of genes that responded to chemical dissociation (163 genes), a stressful experience (0 genes), and cognitive training (423 genes). The three genes that respond to both technical and biological perturbation are Epha6, Grin2a, and Itbp3. B, C) The molecular function of gene ontology (GO) categories that are significantly enriches with either up-or down-regulated genes in response to cellular dissociation (B) or cognitive training (C). The top 10 most significant GO terms are visualized, each with a p-value < 0.001. The fraction next to GO category name indicates the fraction of genes in that category that survived a 10% FDR threshold for significance. Zero terms survived a 10% FDR threshold in response to a stressful experience. We next analyzed gene ontology at 5% FDR significant in each of the data sets to identify 267 the molecular function of genes that changed in response to cellular dissociation (Fig. 4B) 268 or cognitive training (Fig. 4C). The process of cellular dissociation results in a significant up-269 regulation of molecular processes related to ribosomal activity, rRNA binding, oxidoreductase 270 activity, and proton transport, while it caused a down regulation of ligase and helicase activ-271 ity (Fig. 4B). The GO analysis detected no Molecular Function GO terms in the significantly downregulation of oxidoreductase and ribosomal activity (Fig. 4C). Notably, learning induces 279 a downregulation of a structural constituent of ribosomes and oxidoreductase, which were 280 both up-regulated in response to cellular dissociation (Fig. 4B,C).

281
Recovering robust subfield-specific gene expression patterns 282 Using the public Cembrowski et al. (2016b) dataset, we identified 10,751 genes that were 283 differentially expressed between hippocampal sub-regions (Fig. 5A). Using meta analyses of 284 the public Cembrowski data with the primary data described herein, we identified 146 genes 285 that showed robust subfield-specific gene expression patterns (Fig. 5A). Those 146 genes are 286 enriched in cellular compartments related to synapses and molecular functions related to 287 calcium signaling, GTP exchange, and proteoglycan binding (Fig. 5B).
288 Figure 5. Meta analysis of primary and public data. A) This Venn diagram shows the overlap in brain-region specific gene expression across all four experiments (cellular dissociation, stressor habitation, cognitive training, and a public dataset examining subfield comparisons). Grey numbers indicate total number of differentially expressed genes between and two-way subfield comparison. Using this approach, we identified 146 genes that were differentially expressed between any two subfields of the hippocampus in all four experiments. B) Those 146 provide robust brain-region specific markers of gene expression belong to molecular function and cellular compartment GO. The top 10 most significant GO terms are visualized, each with a p-value < 0.05. The fraction next to GO category name indicates the fraction of genes in that category that survived a 10% FDR threshold for significance.

289
The main purposes of this study were 1) to test whether analysis of gene expression in hip-290 pocampus subfields is changed by tissue preparation procedures (cellular dissociation versus 291 homogenization) and 2) to evaluate the effects of a stressful experience relative to cognitive 292 training on analysis of gene expression. The work was designed to evaluate the extent to 293 which technical (i.e. cellular dissociation) and biological confounds (i.e. stressful experience) 294 can impact efforts to assess the transcriptomic response to cognitive processes. This is po-  (Fig. 1D, Fig. 2D, Fig. 3D). The samples that were subjected to 301 cellular dissociation show the least amount of region-specific variation, suggesting that this 302 process might alter the genes that normally distinguish the hippocampal subfields from one of genes with subfield specificity; this is likely due to the cell sorting process that generates 305 a relatively homogenous rather than a heterogeneous population of neurons. These results 306 indicate that cell-type specific differences may be recovered by sorting cells into homogenous 307 populations.

308
Across the data set, many genes consistently or robustly demonstrate hippocampal sub-309 field specificity in their expression (Fig. 5B). The meta-analysis of the primary data and the  (Fig. 1B). This is smaller than the 2.4% (423 /17,320) change we detected in response to 320 cognitive training (Fig. 3B). The stressful experience produces a negligible response (i.e. no 321 significant genes or GO terms were detected), indicating that the mild stress that likely ac-322 companies most behavioral tasks does not have a lasting influence on hippocampal gene 323 expression (Fig. 1B). 324 The extent to which cellular dissociation and unintended stress impacts the expression of 325 particular genes and signaling pathways, limits the feasibility of investigating how genes con-326 tribute to behavior and other responses to organismal manipulations. We found that Grin2a 327 responded to both cellular dissociation and cognitive training (Fig. 4A). Grin2a encodes sub- Alzheimer's Disease (Neuner et al., 2017). 336 We can look beyond the specific genes and examine which pathway responses are con-337 cordant or discordant to multiple treatments. In this case, we saw upregulation of ribosomal 338 activity and rRNA biding in response to cellular dissociation, but we saw an opposing downreg-339 ulation in ribosomal activity and mRNA binding in response to cognitive training (Fig. 5B,C).  We found no detectable transcriptional response in the CA1, CA3, or DG following the stress-347 ful experience (Fig. 2B). The shock experience we used causes a large increase in plasma 348 corticosterone levels, comparable to exposure to predator threat, that is observed after the 349 initial shock exposure session but is absent 24-h later after the second training session (Les- cage controls that do not experience shock, their stress response is indistinguishable from the 353 trained mice but their sensory experience is very different. In contrast, the shock-yoked mice 354 have identical sensory experience as the experimental mice, but they experience stress that 355 the experimental animals do not (Lesburguères et al., 2016). It may be that untrained control 356 mice are optimal, because they would have the identical experience of the environment as 357 experimental mice, except at the brief 500 ms moments of shock that account for roughly 358 3% of the task assuming 20 shocks in 600 s. Depending on the question one control may be 359 preferable over the others, but as demonstrated here, when assessed 24 h after the training 360 experience, they appear to be equivalent in terms of their gene expression profiles (Fig. 2).

362
Many factors contribute to variation in gene expression. We set out to identify the extent to 363 which the process of cellular dissociation -which allows for single cell analysis of neurons 364 -has an appreciable effect on our ability to detect biologically meaningful variation in hip-365 pocampal gene expression. We conclude that there are specific dissociation-induced and 366 cognitive training-induced changes in gene expression that are largely non-overlapping. It 367 is encouraging that the overlap between cellular dissociation and cognitive training is small, 368 indicating that these technical and biological processes affect different transcriptional pro-369 cesses. It is also encouraging to know that the stressful experience had no substantial effect 370 on hippocampal gene expression, which if generalizable to other tasks will allow for using 371 behavioral control groups and behavioral manipulations that also induce modest, potentially