A steady increase in EEG theta power marks RML progression

We have observed that a subset of mice in the commonly used C57BL/6 background are hyperactive at night [6, 17, 18]. In contrast, 129S4 mice are calmer, exhibit more consistent behavior between individuals, and have a very different sleep pattern than C57BL/6 mice [6, 19]. Therefore, to optimize the identification of markers of PrD through gene expression analysis, which could be obscured by bouts of hyperactivity, we used the 129S4 mouse strain for these studies.

The neurodegenerative PrD model was initiated by unilateral, intracranial injection of brain homogenate from either normal (NBH, control) or RML infected mice (Fig 1A). RML induced terminal disease in 129S4 mice following an incubation period of 22.6 weeks, with behavioral features typical for prion diseases (e.g., kyphosis, ataxia, reduced body condition, poor grooming, clasping during tail suspension, etc.) To identify early cellular mechanisms altered by RML, we defined study stages as 1) pre-onset, before any differences and 2) onset, when differences between NBH and RML groups could be detected using a clinically relevant method in living mice (Fig 1B). This was accomplished by measuring EEG changes during sleep-wake cycles, concomitant with core body temperature (Tb) and locomotor activity in freely-moving mice implanted with wireless, telemetric devices (Data Sciences Inc). Twenty-four-hour recordings were taken every two to four weeks until terminal disease was reached. Body weight was measured weekly.

EEG measures the summarized electrical activity of cortical and subcortical neurons, providing a very sensitive tool to detect brain dysfunction [20, 21], including theta power increases in several human PrDs [16]. Analyses of EEG spectral frequency distributions revealed a pronounced and progressive increase in the theta band (5–10 Hz) relative to other frequencies (0–50 Hz) that reached statistical significance at 18 weeks post-injection (WPI) (Fig 1C). Theta power increase was apparent in wake, non-rapid eye movement (NREM) sleep, and rapid eye movement (REM) sleep states (S1 Fig). To examine if EEG theta reflected disease progression, we calculated relative theta power (5–10 Hz/0-50 Hz) against WPI. We found that relative theta power remained consistent in controls yet incrementally increased in RML with disease progression, across all vigilance states (Fig 1C). Mixed model ANOVAs revealed interaction effects of treatment group (RML and NBH) and timepoint (WPI): wake (F_6,48 = 8.3, p = 0.003), NREM (F_6,48 = 6.9, p = 0.011), and REM (F_6,48 = 8.9, p = 0.002). At 18 WPI, a peak effect was seen, reflected by statistically significant higher theta in RML-infected mice compared to NBH mice (Fig 1C). By focusing on the evolution of EEG theta across WPI, we determined that theta power was comparable between RML and NBH mice up to 11 WPI, after which the groups slowly separated (Fig 1C).

Next, we examined brains of RML and NBH injected mice for neuropathological changes. PrDs lead to PrP aggregation (PrP^res) and morphological changes of astrocytes and microglia readily detectable with specific antibodies (PrP, GFAP and Iba1, respectively). Each of these neuropathological hallmarks of PrD were detected in RML-infected mice at 18 but not 10 WPI (Fig 1D). Accordingly, and combined with the EEG theta marker, 18 WPI for RML-infected mice appeared to be a clinically relevant timepoint indicating disease onset as reflected by evident change in brain activity and histology. Nevertheless, mice demonstrated an overtly normal phenotypic behavior, with no differences found in overall sleep parameters, core body temperature, locomotor activity, or response to 6 h of sleep deprivation by the gentle handling method (S1 Fig), and still lacked the typical signs of RML (e.g., kyphosis, ataxia, etc.). Lastly, by combining histological and EEG analyses, we identified 10 WPI as an early pre-onset disease time point. At this timepoint, there was no divergence in theta power or body weight (Fig 1C), and no overt histological changes (Fig 1D) between RML and NBH mice.

Cell types demonstrate specific responses to RML

We directed the expression of epitope-tagged ribosomes [12] to astrocytes and four neuronal types (glutamatergic, GABAergic, PV and SST neurons) using knock-in mouse lines expressing Cre from endogenous Gja1 (Cx43), Slc17a6 (vGluT2), Gad2, Pvalb (PV), and Sst genes, respectively (Fig 2A and 2B). These cell types were chosen primarily as GABAergic neurons, particularly PV-expressing cells, are reported to be the most vulnerable in human and rodent PrDs [4, 22–24]. SST neurons, while generally non-overlapping with PV neurons [25], constitute another important GABAergic subtype. In contrast, glutamatergic neurons are excitatory, driving neural activity, and the overall network excitability is balanced by inhibitory neurons, which are often GABAergic, although some GABAergic neurons are not strictly inhibitory. Lastly, we chose to study astrocytes as they play a key role in the removal of toxic substances from extracellular space [26] and undergo reactive gliosis in PrDs and other NDs [27–29]. An inducible Cre line was employed for astrocytes since constitutive Cre expressing lines are active in neural precursor cells and thus also label neurons [30].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. RiboTag activation and validation.

(A) The RiboTag system is activated in cells with Cre activity through site-directed recombination in Rpl22 to make Rpl22-HA (RiboTag). Frozen mouse brain tissue is homogenized, creating a mix of mRNAs attached to either normal ribosomes (white oval bubbles) or RiboTagged ribosomes (white oval bubbles with pink tags) made specifically by Cre+ cells. Cell-type-specific, ribosome associated mRNA is isolated by RiboTag IP with HA antibody and protein G magnetic beads. mRNA from Cre+ cells is then sequenced. (B) Cell types used in these studies. Shaded areas reflect the expression overlap between Gad2, PV and SST Cre lines. (C) Confocal images of cortex following double immunofluorescence labeling of HA (RiboTag) and representative histological cell-type markers. Cell type markers used (from 1 to 5 respectively): Astrocytes (Sox9), GABAergic cells (Satb2; glutamatergic cell marker to highlight the lack of colocalization with GABAergic neurons), PV neurons (PV), SST neurons (Sst), glutamatergic neurons (Satb2). Examples for double-positive cells are indicated by arrowheads, HA-/marker+ cells by upward pointing arrows, HA+/marker- cells by downward pointing and smaller arrows. Scale bar indicates 50 μm (identical for all panels). C1: specific HA-labeled astrocyte-like shaped cells colocalize with Sox9 in Cnx43-Cre/RiboTag mice. C2: Complete separation of HA and the glutamatergic marker, Satb2, in Gad2-Cre/RiboTag mice was observed. C3: PV and HA show substantial colocalization in PV-Cre/RiboTag mice. C4: All Sst+ cells were co-labeled for HA, although there were several HA+/SST- cells in SST-Cre/RiboTag mice. C5: Nearly all cells were co-labeled for HA and the nuclear marker for glutamatergic neurons, Satb2, in the vGlut2-cre/RiboTag mice. (D) Violin plots showing RiboTag IP specificity; each violin summarizes absolute expression (TPM) ratios (X-axis) of known cell type markers in IP samples and corresponding input (total RNA) samples (Y-axis).

https://doi.org/10.1371/journal.ppat.1010747.g002

We verified the fidelity of RiboTag expression through co-localization of the RiboTag epitope with known cell type markers in histological sections (Figs 2C and S2). RiboTag IP samples and corresponding total mRNA (input control) analyzed by RNA-sequencing (RNA-seq) confirmed cell type-specificity with enrichment or depletion of known cell type marker genes (Fig 2D). Principal component analysis revealed that samples from the same cell types clustered tightly together (S3 Fig).

Since we expected changes to be very mild at 10 WPI, we defined differentially expressed genes (DEGs) to have a Benjamini-Hochberg adjusted p-value (hereafter false discovery rate or FDR) < 0.1 without applying a fold change cut-off. A summary of the DEGs that are unique for each cell type or shared between multiple cells, at each time point, is presented in Fig 3A, and an interactive browser of all DEGs for each cell type is available online (https://shiny.it.liu.se/shiny/Scrapie_RML_RiboTag).Althoughat 18 WPI RML-infected mice appeared normal via passive observation, there were severe gene expression changes in all cell types, with 8 to 40% of all mRNAs altered (FDR ≤ 0.1), and thousands of DEGs having an absolute log₂ fold change|L2FC|≥ 1 (Fig 3B), consistent with gross EEG and immunohistochemistry (IHC) (Fig 1). In contrast, at 10 WPI, SST and PV neurons in RML-infected mice were barely altered (1 and 3 DEGs, respectively, considered unaffected cell types hereafter), contrarily to our predictions based on literature [4, 22–24], but consistent with a recent report [31]. Both glutamatergic and GABAergic neurons showed more changes, with 38 and 83 DEGs respectively, but astrocytes demonstrated the greatest number of DEGs with 139 (Fig 3A and 3B). Glutamatergic DEGs were primarily upregulated, whereas GABAergic and astrocytic DEGs were primarily downregulated (Fig 3A and 3B). Importantly, in contrast to 18 WPI when over half of DEGs were changed in multiple cell types, only 10% of 10 WPI DEGs were changed in multiple cells, indicating that early in disease, cells activate unique responses (Fig 3A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Overview of translatome changes.

(A) Numbers of differentially expressed genes (DEGs). Values represent numbers of genes found to be differentially expressed (FDR ≤ 0.1) in a single or combination of cell types studied. Combinations are marked by dumbbell-shaped symbols below the plots. The shade of grey corresponds to the directionality of changes (see legend above the panel). For 18 WPI all clusters containing < 100 genes were omitted. Note that for 10 WPI, the changes in the 3 affected cell types: astrocytes, glutamatergic cells and GABA cells are mostly non-overlapping (respectively, 120, 66, and 34) (B) Visualization of DEGs (FDR ≤ 0.1; dark blue) and non-DEGs (FDR > 0.1; light grey) in all cell types and timepoints. Y-axis represents log₂ fold change, X-axis represents cell type. (C) Overrepresentation analysis of DEGs in the affected cell types at 10 WPI. Significantly enriched categories were shown for molecular function (MF) and biological process (BP) gene ontology subsets. Categories related to translation, cytoskeletal regulation and circadian rhythm were the most enriched in astrocytes, glutamatergic neurons and GABAergic neurons, respectively. (D) Gene set enrichment analysis (GSEA) in all cell types analyzed at 18 WPI. Significantly enriched categories were shown for the MF gene ontology subset. The most prominent categories were related to translation and terminal oxidation (mitochondrial respiration). Cell type legend for all panels: Cx43: astrocytes, Gad2: GABAergic cells, vGluT2: glutamatergic neurons, PV, parvalbumin neurons, SST: Somatostatin neurons. NES is Normalized Enrichment Score [91].

https://doi.org/10.1371/journal.ppat.1010747.g003

We expected that if cell responses were coordinated, as opposed to being a result of gene expression dysregulation, multiple DEGs for a given cell type would be associated with a specific cellular pathway or process. Since the data from the two timepoints have vastly different numbers and magnitudes of DEGs, we optimized our analysis by using the most appropriate method for the properties of the data. Since the 18 WPI data were characterized by a large percentage of genes being differentially expressed, and consequently multiple functional modalities in the cells were affected, we analyzed those data with gene set enrichment analysis (GSEA). GSEA considers the ranking of all genes (DEGs and non-DEGs), thereby maximizing the potential to determine the most affected modalities (higher gene ranking contributes more to the overall GSEA scoring). In contrast, overrepresentation analysis (ORA) uses unranked DEG subsets which is more appropriate for the 10 WPI data since relatively small changes may not be appropriately reflected in the ranking due to a lower signal-to-noise ratio. Genes were ranked based on the degree of differential expression and log₂fold change (L₂FC) using the formula: rnk = LFC*(−log (FDR)). Noteworthy, when applied to 10 WPI sample data, GSEA failed to yield significant enrichment scores. Corroborating the differences between the cell types with respect to differential gene expression (Fig 3A), this analysis further confirmed that there were no common pathways at 10 WPI (Fig 3C), although the cells responded similarly at 18 WPI, with ribosomal and electron transfer molecular functions in common (Figs 3D and S4). Intrigued by the changes at 10 WPI, we focused the remainder of our analysis on those data. To explore specific cellular changes further, we included interactome information in our analysis.

DEGs occupy neighboring positions in the interactome

By identifying genes in regions of the interactome connected to DEGs, we could deduce which cellular changes are relevant to disease progression. The proteins encoded by genes sharing disease associations often colocalize within protein-protein interactions (PPI) networks, forming communities, referred to as disease modules. Identification of disease modules typically relies on exploring the interactions of known disease genes, i.e., the seed genes. Herein, we defined the seeds as the subset of gene network nodes overlapping with DEGs and examined whether they provide sufficient clues on disease module localization. To this end we collected high confidence PPI interactions (cumulative confidence score > 0.7) from the STRING database and constructed a separate interactome for each cell type, excluding the genes that were not detected in either diseased or control mice. We then determined the sizes of the largest connected components (LCCs) within subnetworks composed of seed genes and compared each to a reference LCC size distribution of randomly selected nodes in the interactome. In astrocytes, GABAergic, and glutamatergic cells the LCC sizes were greater than for random subnetworks of equivalent size (Z-scores 2.53, 5.42, 2.57 respectively, empirical p-values 0.005, 0.001, 0.043) (S5A Fig). The general lack of DEGs in PV and SST neurons was corroborated with a lack of significant connectivity between the top ranked genes (0.01 Wald test p-value cutoff was used) in those cell types (Z-scores 0.30 and 1.69 respectively, empirical p-values 0.342 and 0.077, S5B Fig).

We then sought to identify putative disease modules specific to each affected cell type in early pre-onset disease stage (10 WPI) by integrating the interactome and RiboTag data sets using a modified DIAMOnD method [32]. To account for bias towards well studied genes in the STRING-based network, our modified version of the DIAMOnD algorithm used a two-stage ranking system to iteratively add genes to a disease module (see the supplementary methods section ‘Modification of the DIAMOnD method’). Identified disease modules for all cell types were enriched with genes which share functional annotations with seed genes and therefore represent network communities likely affected by disease. The results are described in the following sections, organized by cell type.

Astrocytes downregulate ribosomal and mitochondrial mRNAs

The most prominent feature of DEGs in astrocytes was a strong down regulation of mRNAs encoding 26 ribosomal proteins, (average L₂FC = -0.64; range -0.41 to -0.94). Additional proteins related to ribosome synthesis (Fbl) or translation in general (Ccdc124, Denr, Dohh, Eif5b, Upf2) were also downregulated, supporting the view that fewer ribosomes were being made. Mitochondrial proteins (Glrx5, Lars2, mRpl18, mRps12, Romo1, Tomm6) and components of the electron transport chain including cytochrome c, and proteins in complexes I (n = 3), III (n = 4), IV (n = 3) and V (n = 3), were all downregulated in astrocytes, except for Lars2 which was upregulated. The most strongly downregulated gene, Pcsk1n (L₂FC -1.25, FDR < 3E-5), encodes proSAAS, a secreted chaperone that has anti-aggregation activity and is found associated with protein aggregates in brains of Alzheimer’s disease (AD) patients and an AD mouse model [33]. Interestingly, RML induced a decrease in macrophage migration inhibitory factor (Mif), within astrocytes. This factor is induced in neurons in AD patients and mouse models [34].

There were also three up-regulated DEGs of interest: C4b, Serpina3n, and Tnrc6a. C4b is a marker for aging astrocytes and signals for inflammation [35], Serpina3n is also a marker for aging astrocytes [35] and is strongly up-regulated in all forms of human PrDs [36], and Tnrc6a associates with a complex involving PrP and Argonaute2 for a role in miRNA-mediated silencing [37].

These observations of 10 WPI DEGs being closely related to ribosome biogenesis and oxidative phosphorylation were confirmed by our ORA analysis (Fig 4A). IHC experiments verified Rps21 was reduced by RML at 10 WPI in cerebellar astrocytes (Fig 4B) which preferentially express Rps21 [6]. Notably, numerous components of ribosomes and mitochondria were reduced in postmortem human PrD brains [38].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Changes in astrocytes (RiboTag::Cx43-Cre).

(A) Overrepresentation analysis of DEGs in astrocytes (Cx43+ cells) at 10 WPI. Significantly enriched GO categories were shown for molecular function (MF), biological process (BP) and cellular component (CC) GO subsets. Most significantly enriched categories were related to translation (e.g., ribosome biogenesis) and terminal oxidation (e.g., electron transfer). X-axis is the reciprocal of gene ratio. (B) RML causes reduced staining of Rps21, one of the ribosomal proteins downregulated at 10 WPI, especially in Bergmann glia between molecular (M) and granular (G) cell layers. (C) Visualization of 10 WPI disease DIAMOnD modules constructed around 64 of 128 DEGs (large circles) in astrocytes. Small circles denote genes added by the DIAMOnD algorithm (see also S6 Fig for the same visualization with gene labels). The module consists of 6 functionally related clusters (identified with the fast greedy modularity optimization algorithm). In each cluster, the genes annotated to the representative GO term are colored. See also S6 Fig for additional GO results.

https://doi.org/10.1371/journal.ppat.1010747.g004

128 DEGs from astrocytes were included in our astrocyte cell interactome and used as seed genes for inference of disease-relevant genes using the modified DIAMOnD, producing a module with 7 major communities. Within the two largest communities, the most overrepresented functional categories were ‘oxidative phosphorylation’ and ‘nuclear transcribed catabolic process nonsense mediated decay’, which consisted primarily of mitochondrial and ribosomal protein gens, respectively. In small communities, the overrepresented categories were ‘proteasomal protein catabolic process’, ‘RNA splicing via transesterification reactions’, ‘heart process’, ‘myeloid leukocyte mediated immunity’ and ‘mitochondrial translational termination’ (Fig 4C). Thus, at 10 WPI, astrocytes primarily downregulated genes for mitochondrial and ribosomal biogenesis in response to RML.

Glutamatergic neurons upregulate cytoskeleton genes

Of the 38 DEGs altered by RML in glutamatergic neurons, many were linked to upregulation of cytoskeleton-shaping proteins (Gsn (gelsolin), Limch1, Mprip, Ppp1r9a) and cytoskeletal components including GFAP and four spectrins (Sptan1, Sptb, Sptbn1, Sptbn2) (Fig 5A). In human PrD, mass spectrometry analysis of synaptasomes revealed Gsn increased while Sptan1 decreased [39]. Gesolin and spectrin breakdown products are common in human NDs [40, 41] and the mRNA upregulation seen here may be compensatory. The GABA receptor protein Gabrb2 was decreased which may lead to increased neuronal activity. Consistent with this notion, Arc (activity-regulated cytoskeleton) commonly induced for synaptic plasticity [42], was increased, which in turn may have led to the increased synthesis of cytoskeleton related proteins. The ORA analysis confirmed this notion with biological process terms including ‘regulation of long-term synaptic depression’, ‘membrane raft organization’, and ‘endoplasmic reticulum to Golgi vesicle mediated transport’ (Fig 5A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Changes in glutamatergic neurons (RiboTag::vGluT2-Cre).

(A) Overrepresentation analysis of DEGs in glutamatergic neurons (vGluT2+ cells) at 10 WPI. Significantly enriched GO categories were shown for molecular function (MF), biological process (BP) and cellular component (CC) GO subsets. Most significantly enriched categories were related to cytoskeleton (e.g., regulation of actin filament organization) and synaptic function (e.g., regulation of long-term synaptic depression). X-axis is the reciprocal of gene ratio. (B) Visualization of disease DIAMOnD modules constructed around 15 of 33 DEGs (circles) in glutamatergic neurons. Pentagons denote genes predicted by the DIAMOnD algorithm. The module consists of 7 functionally related clusters (identified with the fast greedy modularity optimization algorithm). In each cluster, the genes annotated to the representative GO term are colored. See also S7 Fig for additional GO results.

https://doi.org/10.1371/journal.ppat.1010747.g005

We next used the modified DIAMOnD algorithm for inference of disease-relevant genes which incorporated 33 of 38 DEGs in the glutamatergic cell interactome as seed genes. The resulting module contained 10 communities, within which the most overrepresented functional categories were ‘ER to Golgi vesicle mediated transport’, ‘RNA splicing via transesterification reactions’, ‘myeloid leukocyte mediated immunity, ‘endocytosis’, ‘response to peptide hormone’, ‘actin filament-based movement’, ‘humoral immune response mediated by circulating immunoglobulin’ ‘cytoplasmic translation’, ‘post-translational protein modification’ and cell-cell signaling by Wnt’ (Fig 5B). Therefore, at 10 WPI the changes observed in vGluT2+ cells were coordinated and quite different from those in Cx43+ cells.

GABAergic neurons modulate circadian rhythm genes

While RML induced astrocytes to downregulate ribosomal and mitochondrial biogenesis and glutamatergic neurons increased cytoskeleton genes at 10 WPI, GABAergic neurons modulated 11 genes implicated in circadian rhythmicity: Bhlhe41 (Dec2), Ciart (chrono), Dbp, Nr1d1, Nr1d2, Per1, Per2, and Usp2, were down-regulated while, Arntl (Bmal1), Tardbp, and Vip were up-regulated. These differences are related to RML and not time of death, as the latter was controlled for in our experimental design. Remarkably, circadian rhythm genes are commonly dysregulated in neurodegenerative diseases [43]. Many of the downregulated genes normally suppress Arntl expression, potentially explaining its upregulation. Arntl forms dimeric complexes with Clock or Npas2 to drive expression of circadian clock and redox genes, and its disruption in mice leads to neurodegeneration with prominent reactive astrocytosis [44, 45]. Studies have highlighted a role for the circadian rhythm network in antioxidation and neuroprotection [44] and we see a potential connection in our data with down regulation of proteins involved in mitochondrial morphology (ROMO1) and electron transport complex III (Uqcc2 and Uqcrq). It is worth noting that timing may be modified by RML independently to sleep if driven by GABAergic cells.

The most significant pathways from ORA included ‘RNA catabolic process’, and ‘circadian regulation of gene expression’ for BPs, ‘postsynaptic density membrane’ for CC, and ‘structural constituent of ribosome’ for MF (Fig 6A). Analysis of 10 WPI data with the modified DIAMOnD resulted in a module built on 69 of 83 DEGs with 6 communities with functional categories of ‘ATP synthesis coupled electron transport’, ‘cotranslational protein targeting to membrane’, ‘RNA splicing’, ‘protein polyubiquitination’, ‘rhythmic process’, and ‘myeloid leukocyte mediated immunity’ (Fig 6B). Thus, although Gad2+ cells also changed expression of some ribosomal and mitochondrial proteins, they predominately changed additional pathways and demonstrated an overall response quite different from Cx43+ and vGluT2+ cells.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Changes in GABA neurons (RiboTag::Gad2-Cre).

(A) Overrepresentation analysis of DEGs in GABAergic neurons (Gad2+ cells) at 10 WPI. Significantly enriched GO categories were shown for molecular function (MF), biological process (BP) and cellular component (CC) GO subsets. Most significantly enriched categories were related to circadian rhythm (e.g., circadian regulation of gene expression) and translation (e.g., cytosolic large ribosomal subunit). X-axis is the reciprocal of gene ratio. (B) Visualization of disease DIAMOnD modules constructed around 27 of 69 DEGs (large circles) in GABAergic neurons. Small circles denote genes added by the DIAMOnD algorithm. The module consists of 6 functionally related clusters (identified with the fast greedy modularity optimization algorithm), the largest of which involved ‘rhythmic process’. In each cluster, the genes annotated to the representative GO term are colored. See also S8 Fig for complete GO results.

https://doi.org/10.1371/journal.ppat.1010747.g006

Identification of EEG biomarker that precedes morphological disease changes

A better insight into mechanisms of selective vulnerability in neurodegenerative diseases is needed to develop therapies. For these therapies to stop or slow disease progression, it is important to understand disease mechanisms at early stages, before clinical signs. We chose mouse RML prion disease as a model as the timeline from induction until clinical illness and neuropathological changes is highly precise [46]. To define disease onset, we measured theta waves longitudinally with EEG. As in human PrDs [16] we found theta increased. Despite spectral changes in this band, and other spectral composition alterations at 18 WPI, we did not find changes in sleep-wake states, which suggests for this RML model, that degeneration remains too mild to affect sleep or that compensation mechanisms are in place that protect sleep-wake behavior up to 18 WPI. Nonetheless, the observed increase in theta as disease progressed might be suitable for detecting effects of disease-modifying experimental treatments in mice.

RiboTag elucidates cell type-specific responses to early ND

RML in RiboTag mice permitted the detection of early molecular responses elicited by different cell types. To verify if the RiboTag IPs contained mRNAs from the cells of interest we assessed the enrichment or depletion of known cell type marker genes compared to the total RNA input. Cell type markers are generally not absolutely specific for one cell type, but high abundance can predict the population of interest. GFAP is a notable example; as it is widely recognized as an astrocyte marker, but we found expressed in our neuronal samples, and RML induced a significant upregulation of GFAP in glutamatergic neurons at 10 WPI. Neuronal expression of GFAP has been reported by chromatin marks in healthy young adult mouse brains [47], human AD by IHC [48] and single nucleus RNAseq [49], in motor neurons in a mouse model of amyotrophic lateral sclerosis [50], and in synaptasomes of human PrD [39]. Therefore, increased GFAP expression in neurons seems common to many neurodegenerative diseases. Furthermore, since several other astrocyte markers were strongly depleted/absent in neuron-derived samples, we conclude that this GFAP observation was not due to cross-contamination.

Unexpectedly, at 18 WPI there were extensive changes in gene expression even though by passive observation RML-infected mice had no overt phenotype. This implies that even when the first clinical symptoms appear in patients of NDs, there may already be significant cellular perturbations, across several cell types making successful treatment interventions difficult. Furthermore, at 10 WPI, prior to EEG and histological abnormalities, we identified that astrocytes, glutamatergic, and GABAergic cells (presumably PV- and SST-) were responding. While at 18 WPI each cell type changed similar molecular pathways and sets of genes, at 10 WPI each population demonstrated unique DEGs and pathways and was, therefore, altered specifically. These results support our hypothesis that the phenomenon of selective vulnerability is caused by specific brain cells adopting unique responses to cope with an emerging neurodegenerative disease.

Distinct translatome responses in glutamatergic and GABAergic neurons

Glutamatergic neurons at 10 WPI altered their translatomes to maintain or enhance actin and cytoskeleton networks, which are important for cellular transport to the synapse and signal transduction. Moreover, Arc responds to synaptic activity by driving expression of cytoskeletal proteins at synapses [51], and its upregulation suggests glutamatergic neurons are abnormally active, which may relate to excitotoxicity [52], a molecular mechanism long postulated to be involved in neurodegenerative diseases.

GABAergic neurons prominently modulated their translatomes for the circadian rhythm pathway, potentially corroborating findings for other NDs in humans [43, 53, 54]. Interestingly, a method analagous to RiboTag also found circadian rhythm genes prominately changed in a series of mouse models of Huntington’s disease [55]. Future experiments are needed to reveal why the circadian rhythm pathway is modulated during neurodegeneration.

Astrocytes reduce synthesis of ribosomes and mitochondria

Human PrDs are characterised by reduced synthesis of mitochondria and ribosomes [38, 56], and our cell-type-specific data suggest that during early disease, these changes occur in astrocytes, even before other signs of disease. There are several possible explanations for why astrocytes deploy this strategy.

The strong reduction of ribosome protein mRNAs could be related to ribosome specialization, where the stoichiometric ratio of ribosomal proteins differs between tissues and brain regions, possibly to steer the preferential translation of specific mRNAs [2, 6, 57–59]. However, the reduction in many ribosome proteins in combination with reduced mitochondrial proteins suggests it is not simply ribosome specialization. Since ribosomes are highly abundant and energy intensive to make, their synthesis is tightly coupled to that of mitochondria, and the reduction in both implies that complete ribosome synthesis is decreased. In the mouse brain, the half-life for ribosomal and mitochondrial proteins is approximately 10 and 20 days, respectively [60]. Thus, while acute intracellular signaling (e.g., phosphorylation), occurs without mRNA expression changes, drivers of long-lasting reduction in mitochondria and ribosomes will likely be directed by gene expression changes, which might be reflected in our data.

Glycogen metabolism, a key function of astrocyte energetics, was not affected as adenylate cyclases, phosphodiesterases and other genes associated with protein kinase A signaling were unchanged. The mTOR (mechanistic target of Rapamycin) pathway also senses energy levels and controls growth, in part through coordination of ribosome biogenesis. However, the only DEG related to mTOR, Fkbp8, is an inhibitor and its decrease should therefore stimulate ribosome biogenesis. Thus, a lack of glycogen and mTOR signaling changes suggest there is not a change in metabolism.

The unfolded protein response acts to relieve protein misfolding stress by reducing overall protein synthesis, especially ribosomes, and increasing production of protein quality control proteins (e.g., chaperones) [61]. Such a response was noted in an over-expression transgenic model of acquired prion disease [62, 63]. However, data from another study [31] and ours did not demonstrate increased synthesis of protein quality control genes, suggesting that an artifactual unfolded protein response may have been triggered in the overexpression transgene study [63].

A final explanation for the reduced ribosome synthesis in our study relates to molecular crowding. Cells tune ribosome abundance to sculpt their physical properties, including diffusion rates of intracellular protein aggregates or cellular stiffness [64, 65]. It is conceivable that astrocytes would benefit from such a manipulation in early stage PrD.

Comparison with previous studies

Previous reports have also contributed to our understanding of selective vulnerability by applying alternative methods. The laser capture microdissection approach was used to cleanly remove sections of the CA1 region of the hippocampus or granule cells of the cerebellum from mice at various stages of PrD [66, 67]. Purified RNA was then measured with microarrays. These regions primarily, or exclusively, contain glutamatergic neurons and are thus most comparable to our glutamatergic RiboTag samples. Interestingly, they found RML upregulated genes related to actin-cytoskeleton in both regions early in disease, with a stronger response in the cerebellum. Our data also hint at the cerebellum being the location of affected glutamatergic neurons as some of the DEGs (e.g., Zic1 and Zic2) are preferentially expressed there, according to the Allen Brain Atlas.

Another study employed fluorescent activated cell sorting to study microglia from prion diseased mice [68], a cell type providing an important neuroprotective function to brains during prion diseases [69]. Interestingly, microglia sharply increase expression of mitochondrial and ribosomal protein encoding genes [68], in direct contrast to what we report here for astrocytes. This difference may relate to the concept that ribosome abundance is employed to regulate intracellular crowding [64, 65] and/or due to differences in phagocytic capacity by micro- and astroglia during PrDs [70].

During the completion of our work a report of experiments like ours was published [31]. They employed a Cre dependent ribosome tagging system analogous to RiboTag to study the responses of excitatory and PV neurons, astrocytes, and microglia to intraperitoneal RML injection. Despite the overall similarity of the experiments, the results in the two studies were surprisingly different. For example, in their study 24 WPI corresponded to 75% of complete disease progression, with less than 300 DEGs, whereas the 18 WPI timepoint of our study corresponded to 80% of disease progression, with thousands of DEGs, many with |L2FC| >> 1. Furthermore, we detected over 200 DEGs at 10 WPI representing 44% disease progression, while they detected no DEGs at 16 WPI, representing 50% disease progression. Last, they detected no changes in astrocytes and microglia until late in disease, and none in neurons until mice were near death. One explanation for the discrepancy is the application of a |L2FC| > 1 cutoff in their study, which we did not apply since we expected only small changes in mild disease stages. We consider it likely that mRNAs changing less than 2-fold can have substantial biological impact. For example, given how abundant they are in the first place, it is easy to envision how the coordinated reduction of ribosomes and mitochondria in astrocytes at 10 WPI may impact on cellular metabolism, or as discussed earlier, intracellular crowding. Additional explanations for the different outcomes between the studies are technical in nature, including different mouse genetic backgrounds, different treatments of RNA (e.g., the application of RNase in their procedures), different transgenic mouse lines, etc. All explanations ought to be considered when comparing these studies.

Limitations of this study

As with most experiments of this size and complexity there are some limitations of our study. One limitation is that death of neurons could give the impression that surviving cells have changed their translatomes even if they have not. Although the RML prion strain has been reported to not cause significant neuronal loss even at disease stages further progressed than those we studied [31], it is conceivable that a subpopulation of neurons that died went undetected. Future experiments employing stereology methods, in combination with markers from our data, could be used to find these missing cells. A second limitation is that the capture of ribosomes informs on the mRNAs undergoing translation, but misses other mRNAs, such as those in processing bodies or sequestered in granules. To satisfy the need to capture total RNAs and ribosome-associated mRNAs from the same cells of the same mouse, the Tagger mouse line can be employed [71]. Another limitation is that although PrP^res was not detected at 10 WPI using conventional IHC methods, it is likely that more sensitive methods would demonstrate a small amount of PrP^res was present. Furthermore, the GFAP labeling at 10 WPI could have missed important morphological changes to astrocytes which our data suggest may be present. Despite these limitations, the conclusion that different cell types manipulate their translatomes in a coordinated fashion at a very early stage of prion disease is strong.

Outlook

With these results in hand, it will now be interesting to study additional, early time points, especially considering the drastic changes in gene expression at 18 WPI, as well as selected brain regions to further study localized subgroups of astrocytes, GABAergic and glutamatergic neurons. Furthermore, in addition to actively translated mRNAs, other levels of gene expression including miRNAs and lncRNAs as well as epigenetic features could be analyzed from specific cell types using the Tagger mouse line, analogous to the RiboTag system, but expressing four tagging modalities [71]. Studying multiple levels of gene expression in a cell-type-specific manner will greatly enhance our understanding of disease mechanisms and may point towards therapeutic targets. The work described here, validating such a system in our mouse model in the context of an EEG-based early biomarker of disease progression, is an important first step in this direction.

Ethics statement

Ethical permissions for animal work were granted by the Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen, permission #84–02.04.2013.A128 and #84–02.04.2013.A169. All experimental procedures were performed in accordance with the internal regulations of the German Center for Neurodegenerative Diseases and Linköping University. Surgical anaesthesia was accomplished with isoflurane and mice were sacrificed by aspyxiation from carbon dioxide.

Mouse lines and husbandry

Mice were housed in individual ventilated cages (22°C; 60% humidity) on a standard 12h light / 12h dark cycle and had access to food and water ad libitum. RiboTag mice (Jax #11029; [12]) were crossed to cell-type-specific Cre driver lines: Vglut2-Cre (Jax #016963; [72], Gad2-Cre (Jax #010802; [73], SST-Cre; Jax #013044, PV-Cre (Jax #008069; [74]) and Cx43-CreER^T a gift from Martin Theis, [75]. To induce Cre activity in Rtag/Cnx43-CreER^T mice, starting one week before sacrifice, Tamoxifen (T5648, Sigma-Aldrich) prepared as 10 mg/ml in 10% ethanol / 90% sterile sunflower oil, was intraperitoneally injected, 10 μl/g body weight, once per day for three days. Experimental mice were heterozygous for both RiboTag and Cre recombinase. For each condition (5 RiboTag/Cre combinations, injected with 20 μl 0.1% RML or NBH, sacrificed at 10 or 18 WPI) four mice were used (total n = 80). Mice were genotyped for the presence of Cre and the RiboTag allele using an established PCR protocol [12]. Prior to intercrossing, each Cre and RiboTag mouse line was bred to S4 for at least 8 generations and were verified to be at least 99% 129S4 [6] by a commercial high throughput single nucleotide polymorphism service (Envigo/Harlan).

NBH/RML prion injections

The RML prion homogenate used in these experiments was originally acquired from Gregory Raymond (Rocky Mountain Labs). It was injected into 129S4 mice at the location where this line originated, the Whitehead Institute for Biomedical Research, Cambridge, Massachusetts, resulting in an incubation period to terminal disease of 22.8 weeks. Brain homogenates prepared from these RML-infected 129S4 mice resulted in an incubation period to terminal disease of 22.6 weeks and were the source of RML for the EEG and RiboTag experiments reported here. All three experiments included a common injection procedure: 20 μl 0.1% brain homogenate from normal brain (NBH) or RML-infected mice was injected into the right brain hemisphere at the bregmatic suture using hand-held syringes with needle guides to control injection depth of approximately 3.5 mm.

EEG and sleep recordings

For recordings of EEG, EMG, body temperature, and locomotor activity, ten male wild type 129S4/J (Jax #009104) mice were bred in-house and injected with NBH (n = 5) or RML (n = 5). Only males were used to reduce the variability of sleep parameters influenced by sex [76, 77]. They were injected at 1 year of age since the study was done in parallel with a similar study of inherited prion diseases that required older mice (manuscript in preparation). Mice injected at this age die from RML only 5% faster than young mice [78]. Surgical procedure: After at least twelve and a maximum of 22 d incubation time, mice were implanted intraperitoneally with F20-EET transmitters (Channel bandwidth 1–50 Hz, Data Sciences International (DSI)) under isoflurane anesthesia. EEG leads were routed subcutaneously to the skull, placed epidurally above the left frontal cortex (AP: 1.5, ML: 1.5 mm from bregma, negative lead) and the right parietal cortex (AP: -2.5, ML: 2.0, positive lead) and fixed in place with dental acrylic. EMG leads were anchored in the neck muscles. Mice were allowed at least two weeks of recovery before recordings. Data acquisition and analysis: Undisturbed 24 h baseline recordings were performed in the home cages. Mice were kept in individual cages but could hear, smell and see at least one other mouse of the same experiment. The cages were placed in ventilated cabinets with a 12 h light / 12 h dark cycle. Sleep scoring and analysis was done as reported before [79–81]. EEG and EMG were recorded via telemetry using DQ ART software (DSI). Sampling frequencies were 500 Hz. EEG low-pass filter cut off was 100 Hz (in addition to the 1 Hz high pass and 50 Hz low pass anti-aliasing filtering built into the transmitter). EEG and EMG recordings were scored in 10 s epochs as wake, rapid eye movement sleep or non-rapid eye movement sleep by an expert scorer who examined the recordings visually using NeuroScore 3.0 software (DSI). EEG spectra were analyzed with a fast Fourier transform algorithm using a Hanning Window without overlap (NeuroScore) on all epochs without stage transition or artifact. For direct comparisons of EEG power spectra, power was expressed as relative power, i.e., each frequency bin (0.122 Hz) was divided by the sum of the values between 0 and 50 Hz. Spectra were compared using permutation ANOVAs with 5000 iterations (Manly). Relative theta power was calculated as the power between 5 and 10 Hz (summed values of the respective frequency bins) divided by the sum of the values between 0 and 50 Hz.

Tissue preparation

For gene expression studies 2-4-month-old RiboTag mice were injected, 4 mice each for NBH and RML. Only females were used for Gad2, PV, and Cx43 lines; both sexes were used in nearly equal proportion for vGlut2 and SST lines. RiboTag mice were sacrificed at 10 or 18 WPI by CO2 asphyxiation. Brains were removed from the skull, the olfactory bulbs were removed and discarded, the brain stem was cut at the level of the most caudal tip of the cerebellum, the two hemispheres were separated; one was snap frozen for mRNA isolation and the other formalin fixed for IHC. Prion infectivity in formalin fixed hemispheres was inactivated with formic acid prior to embedding in paraffin cassettes. Histology samples from mice injected with NBH were treated the same way in parallel. Additional technical details are in the supplementary materials section “Immunohistochemistry”.

RiboTag translational profiling

Each frozen brain hemisphere was weighed and put into a tissue grinder with a corresponding volume of polysome buffer (50 mM Tri pH = 7.5, 100 mM KCl, 12 mM MgCl2, 1% Nonidet P-40, 1 mM DTT, 1 x Protease inhibitor (SigmaFast Protease Inhibitor Tablets (S8820), Sigma-Aldrich), 100 U/ml RNAse inhibitor (RNAase inhibitor (N8080119), Thermo Fisher Scientific), 100 μg/mL cycloheximide (Cycloheximide (C7698), Sigma-Aldrich) to prepare a 10% brain homogenate. Homogenate was prepared with help of a motorized tissue grinder (Heidolph, 600rpm, ~30s) on ice and centrifuged (10000 g, 10 min, 4°C). Supernatant was transferred to a new reaction tube and used as IP input and for isolation of total RNA. For IP with magnetic beads (Dynabeads Protein G 10004D, Thermo Fisher Scientific) pre-cleared supernatant (25 μl beads and 200 μl supernatant, 30 min, 4°C) was first incubated with anti-HA antibody 12CA5 (Roche Life Science; 200 μl pre-cleared supernatant, 10 μl antibody, 45min, 4°C) and this mixture then added to the magnetic beads (50 μl beads, 1–2 h, 4°C). Magnetic beads were washed three times with PBS (~500 μl) before use and incubation steps of IP were done on a rotator. IP samples were put on a magnetic rack and the magnetic bead pellets were washed three times with high salt buffer (50mM Tris pH = 7.5, 300 mM KCl, 12 mM MgCl2, 1% Nonidet P-40, 1mM DTT, 100 μg/ml cycloheximide; ~500 μl). Cell-type-specific mRNA was eluted from the magnetic beads by RLT buffer supplemented with 2-mercaptoethanol from the RNeasy Mini Kit (Qiagen; 200 μl; Thermomixer: 700 rpm, 5–10 min, RT) and afterwards isolated with this kit. For each supernatant two technical replicates of IP were done, and total RNA from the input supernatant (200 μl) was isolated in parallel. Quality and quantity of immunoprecipitated mRNA and total RNA were verified by Qubit Fluorometer (Thermo Fisher Scientific) and Agilent 2100 Bioanalyzer (Agilent Technologies).

RNA-seq library preparation and RNA-sequencing

300 ng of total RNA or 150 ng of immunoprecipitated mRNA was used for RNA-sequencing. Each sample was checked for amount and quality using a Nanodrop 2000 (Thermo Fisher Scientific) and Agilent 2100 Bioanalyzer. For each condition (5 Cre/Rtag lines, injected with RML or NBH, sacrificed at 10 or 18 WPI) we used four individual samples (exceptions: SST RML 10 WPI: n = 3; PV NBH 10 WPI: n = 3; Astro RML 18 WPI: n = 2). RNA was converted to cDNA using the Transcriptor High Fidelity cDNA Synthesis Kit (Roche Aplied Science). RNA-sequencing libraries were prepared using the TruSeq RNA Sample Preparation 2 Kit (Illumina). The library quality was checked using an Agilent 2100 Bioanalyzer and concentration was measured by a Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific) and adjusted to 2nM before sequencing (single end, 50 bp) on a HiSeq 2000 Sequencer (Illumina) using TruSeq SR Cluster Kit 3-cBot-HS (Illumina) and TruSeq SBS Kit 3-HS (Illumina) according to the manufacturer’s instructions.

RNA-seq data analysis

Quality assessment was based on the raw reads using the FastQC (0.10.1, Babraham Bioinformatics) quality control tool. The sequence reads (single-end, 50 bp) were aligned to the mouse reference genome (mm10) with Bowtie2 (2.0.2) [82] using RSEM (1.2.29) [83] with default parameters. First, the mouse reference genome was indexed using the Ensembl annotations (84.38) with rsem-prepare-reference from RSEM software. Next, rsem-calculate-expression was used to align the reads and quantify the gene abundance. Differential expression analysis was carried out using total gene read counts with DESeq2 package (1.12.4) [84]. Genes with less than five reads (baseMean) were filtered out and false discovery rate (padj/FDR) was recalculated with Benjamini-Hochberg procedure for the remaining genes.

Overrepresentation analysis (ORA) and gene set enrichment analysis (GSEA)

ORA was performed with the enricher() function from clusterProfiler R package [85]. GSEA was done using WebGestallt [86]. Gene sets were sourced from Molecular Signature Database 3.0 [87] and sets containing less than 10 or more than 500 genes were excluded from the analysis. For ORA, all annotated genes expressed in the analyzed disease time point and cell-type were used as background gene set. Redundant categories were filtered out using simplifyEnrichment R package (v.1.3.1) (doi: https://doi.org/10.1101/2020.10.27.312116) and visualized using ggplot2 R package (v.3.3.5) (https://ggplot2.tidyverse.org).

Gene network construction for disease module analysis

Protein-protein interaction network for Mus musculus was sourced from the STRING database (v. 11.0) [88]. ENSEMBL protein IDs representing node names were replaced with corresponding ENTREZ gene IDs using the getBM() function (biomaRt package; version 2.64.3). Data were filtered for interactions with combined confidence score ≥ 0.7 and a separate graph object was created for all cell types (igraph package; version 1.2.6), including only the vertices with corresponding Wald test p-values (genes that were not expressed in either RML or NBH mice were removed). Topological characteristics of all generated networks are presented in S1 Table. Interactome information was sourced from the STRING database (v. 11.0) Szklarczyk, Gable [88]. The interactome igraph object used in this study is available from GitHub (https://tinyurl.com/tcyryath).

Seed clusters validation for disease module analysis

We have considered the genes with FDR < 0.1 as members of the seed cluster and used a network-based approach to infer cell-type-specific disease modules. For a given set of seed genes LCC size was measured with components() function (igraph R package) and compared to a reference distribution obtained by measuring the LCC size of 10,000 random node subsets of size equal to that of the original seed gene set. Randomization was done so that node selection was limited to genes with degrees comparable to those of the seed cluster. To avoid repeated selection of high degree nodes (in biological networks such nodes are rare, thus many have unique degree values), nodes within a specified degree interval were organized into bins such that there were at least 100 nodes in each bin. Random sampling of nodes was then performed from bins containing seed genes (degree-preserving randomization). Statistical significance of the observed LCC size was assessed by calculating the empirical p-value (the fraction of random sets of nodes with LCC size equal to or greater than the size of LCC formed by seed genes), and Z-scores (Z-score = (LCC_O—mean(LCC_R))/sd(LCC_R), where LCC_O is observed LCC, LCC_R is random LCC and sd is standard deviation).

Construction and validation of cell type-specific disease modules

Candidate disease module genes were identified for selected cell types using the modified DIAMOnD algorithm, consisting of the following steps: 1) for all genes connected to any of the seed nodes, Wald test p-values (the minimum value was selected in cases where multiple Wald test p–values corresponded to the same ENTREZ gene ID) and connectivity significance were determined, 2) candidate proteins were ranked based on their respective Wald test p-values and connectivity significance (ties within individual rankings were resolved using fractional ranking), 3) Individual rankings were combined into a single score given as the sum of the inverted ranks. 4) the gene with the highest combined score was included into the disease module and became a member of the new seed cluster, 6) steps 1–5 were repeated 500 times. The threshold for the total number of module genes was set as the iteration at which the number of genes with Wald test p-value < 0.01 in the module reached a plateau (S5B Fig). The relevance of individual thresholds was verified by determining the number of module genes sharing functional annotations with seed nodes. Annotated gene sets from Gene Ontology (GO), KEGG, Reactome and Biocarta were retrieved using the msigdbr R package (v7.4.1) and all GO terms/pathways significantly enriched within a given set of seed genes were identified using Fischer’s exact test (Benjamini–Hochberg adjusted p-value < 0.05) (Gene sets containing less than 10 and more than 500 genes were excluded from the analysis). Module genes annotated to any of the identified GO terms/pathways were considered true positives. The statistical significance of the number of true positive genes at a given iteration step was determined using Fischer’s exact test). Gene modules were validated using annotated gene sets from Gene Ontology (GO), KEGG, Reactome and Biocarta were retrieved using msigdbr R package (v. 7.4.1) and used in ORA analysis using clusterProfiler R package. Genes annotated to any of the identified GO terms/pathways were included in the validation gene set. The statistical significance of the number of functionally related module genes at a given iteration step was determined using one-sided Fisher’s exact test.

Availability of supporting data

Data generated or analysed during this study are included in this published article and its supplementary information files. The raw and preprocessed data have been deposited in NCBI’s Gene Expression Omnibus [89] and are accessible through GEO (GSE189527; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE189527). The interactome igraph object used for network analysis can be downloaded from GitHub: https://tinyurl.com/tcyryath. Scripts used for data analysis are available from GitHub: https://tinyurl.com/5sbexp89. An Interactive graphical browser of all RiboTag data for each cell type is available online: https://shiny.it.liu.se/shiny/Scrapie_RML_RiboTag.