Distinct translatome changes in specific neural populations precede electroencephalographic changes in prion-infected mice

Selective vulnerability is an enigmatic feature of neurodegenerative diseases (NDs), whereby a widely expressed protein causes lesions in specific cell types and brain regions. Using the RiboTag method in mice, translational responses of five neural subtypes to acquired prion disease (PrD) were measured. Pre-onset and disease onset timepoints were chosen based on longitudinal electroencephalography (EEG) that revealed a gradual increase in theta power between 10- and 18-weeks after prion injection, resembling a clinical feature of human PrD. At disease onset, marked by significantly increased theta power and histopathological lesions, mice had pronounced translatome changes in all five cell types despite appearing normal. Remarkably, at a pre-onset stage, prior to EEG and neuropathological changes, we found that 1) translatomes of astrocytes indicated reduced synthesis of ribosomal and mitochondrial components, 2) glutamatergic neurons showed increased expression of cytoskeletal genes, and 3) GABAergic neurons revealed reduced expression of circadian rhythm genes. These data demonstrate that early translatome responses to neurodegeneration emerge prior to conventional markers of disease and are cell type-specific. Therapeutic strategies may need to target multiple pathways in specific populations of cells, early in disease.


Introduction
NDs typically involve the misfolding and aggregation of proteins which lead to toxicity and cellular stress. Most NDs are characterized by late onset of clinical signs and, conceivably, protective mechanisms must be proactive to counterbalance the toxic processes that eventually manifest as disease symptoms. Neurodegeneration progresses when the capacity of those protective proteostatic mechanisms is exceeded [1]. NDs initially affect defined regions in the brain, a phenomenon known as selective vulnerability [1][2][3][4][5]. Expression of disease-causing proteins is typically ubiquitous, and often higher in regions and cell types resistant to disease [2,6], indicating the observed regional and cellular selectivity is determined by other factors. We hypothesize that the differences in cellular microenvironments (e.g., specific compositions of the protein trafficking and quality control machineries) influence the vulnerability of individual cells, which adopt distinct responses to cope with the disease-related protein conformers [1,2,7]. As brain regions are largely composed of neural networks consisting of multiple functionally interacting cell types, failure of one cell type likely triggers a domino effect that results in region-wide changes observed in human NDs [2].
In this study we dissected responses of five neural cell types to proteostatic stress using a mouse model of neurodegenerative PrD based on the mouse adapted prion strain passaged at the Rocky Mountain Labs (the model is hereafter referred to as RML). In contrast to other NDs, PrDs naturally affect several mammalian species, and mouse models present neuropathological, biochemical and gene expression changes very similar to humans. Moreover, knowledge gained by studying prion diseases is important for other NDs since the prion protein (PrP) that causes PrDs, has a mechanistic role in them [8] and many are thought to spread through the brain via prion-like mechanisms [9][10][11]. RML is an ideal model since it is extraordinarily precise, causing lethality to an entire study group within a week following a fivemonth disease course. Such precision in disease progression is well suited for gene expression analyses. We thus studied cell type-specific responses to RML using RiboTag profiling [12]. RiboTag is a Cre-directed transgenic tool that enables the isolation of ribosome-associated mRNA from cells expressing an epitope-tagged ribosomal protein (Fig 1A) [12]. By capturing actively translated mRNAs rather than total RNA, RiboTag data more closely reflect the changes in the cellular proteome of specific cell types [13]. Moreover, tissue samples can be immediately frozen upon removal, constraining batch effects and preserving mRNA significantly better than methods based on physical separation of cells or cell bodies [14,15]. To determine the most informative analytical timepoints, we tracked theta frequency waves, which increase in human PrD, using longitudinal EEG analysis [16]. Changes to gene expression were analyzed for five brain cell types at the disease onset stage, when theta was significantly increased, and at an earlier, pre-onset stage when theta for both groups was still identical and neuropathological changes were absent.
Even at the pre-onset stage, three cell types demonstrated coordinated, but distinct, responses although the cell type we most expected to respond did not. However, all five cell types were strongly altered at the later, disease onset stage, which coincided with increased EEG theta power. Some of the differentially expressed genes (DEGs) from the early time point have been related to PrDs and other NDs; here we present such alterations in the context of specific cell types. More significantly, this work indicates that disparate cell types unleash unique strategies in response to neurodegeneration and that parallel engagement of multiple therapeutic targets may provide the most efficacious treatments.

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). Twentyfour-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][23][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][28][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].
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 2 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 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/ and 3 DEGs, respectively, considered unaffected cell types hereafter), contrarily to our predictions based on literature [4,[22][23][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).
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 2 fold change (L 2 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 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 twostage 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 2 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 2 FC -1.25, FDR < 3E-5), encodes pro-SAAS, 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]. 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 2 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

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice 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).
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

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice 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 DIA-MOnD 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.

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice

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

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice 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 crosscontamination.
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][58][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 microand 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

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice 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.

Conclusion
In humans and several mouse models, PV neurons show selective vulnerability in late or terminal stages of PrDs. This study shows that, in response to early stage acquired PrD, PV cells

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice do not modify their translatomes, even though other cell types do by modulating distinct pathways. These changes occur even when the EEG profiles of control and experimental mice are indistinguishable. Once the EEG profile becomes abnormal, the translation of many thousands of genes is already affected in each cell type, even though mice still appear normal. If the same occurs in humans, once a patient's brain has degenerated to the point that EEG signs are present, gene expression will be so drastically changed that multiple therapeutic targets will need to be engaged. Pathways we found to be engaged early by resistant cells are composed of a narrow and testable set that might be therapeutically beneficial to PV neurons.

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.

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][80][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 Neuro-Score 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".

PLOS PATHOGENS
RiboTag and EEG analysis of prion-infected mice 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 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

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 DIA-MOnD 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. Examples of frequency power spectra during wake, acquired serially as disease progressed. The power in each frequency bin is expressed as percentage of the cumulative power of all frequencies (0-50 Hz). Curves depict group averages, shaded areas depict SEM. The degrees of freedom are 409 and 3272. The p-values for post hoc uncorrected bin-bybin t-tests are indicated below the spectra. The average wake power spectra up to 11 WPI (C) shows no difference between the groups but at 13 WPI (D) differences begin to emerge and at 18 WPI (F) RML-infected mice show a significant increase in relative power in the theta range (5-10 Hz) and a decrease in relative power in lower and higher frequencies. A summary of data in panels A-G is depicted in H, also shown in Fig 1C. A similar increase in theta at 18 WPI was also seen for NREM (J) and REM (L), but not at 5 WPI (I, K). Principal Component Analysis (PCA) of RiboTag data from all cell types and time points analyzed. PCA was done on variance-stabilized normalized counts (vst() transformation using DEseq2 R package). Clustering was dominated by differences between cell types rather than disease states. (B) Heatmap showing RiboTag specificity. Each row corresponds to a cell type marker for one of the cell types listed on the left side of the heatmap. Each column represents one biological replicate (one mouse), and the specification of the sample is encoded by the colored legend below the heatmap. Input samples are total RNA collected specifically for this analysis for comparative purposes. Z-score for each row was calculated so that the mean input level for each row was set to 0 (Z = (x-row_mean input )/row_SD, where row_SD is row standard deviation and x is the TPM value). In case of PV neurons (D), no terms enriched at FDR < 0.1 were found and therefore all reported terms were shown. Categories with GSEA FDR < 0.1 were shown after removing similar terms based on semantic similarity using GOSemSim R package. Categories related to terminal oxidation (highlighted in red) and translation (highlighted in blue) were amongst the most enriched, and were common across multiple analyzed cell types, suggesting that the 18 WPI timepoint is marked by canonical changes to mitochondria and ribosomes, typical for neurodegenerative diseases. nodes, s-number of seed nodes in the final module, d-edge density, k-number of edges, ks-number of edges including predicted and seed nodes (ks). The modules produced by the modified algorithm had more seed genes and a higher percentage of connections between initial seed genes and predicted genes than those obtained from the original algorithm. (TIFF) S1 Table. Topological characteristics of all generated networks are presented. In contrast to a previous report [90], we observed that node removal does not have a major impact on network connectivity. As indicated by the number of connected components, filtered networks are more fragmented in comparison to the generic network. However, most of the genes constitute a single connected component (96.9% on average). The diameter and mean shortest distance did not increase with node removal, suggesting that the "small world" property of the original network was preserved. (PDF) S2 Table. Seed gene information and validation results. Only a fraction of seed genes have corresponding vertices in the PPI network. Most of the seeds occupy neighboring interactome positions, as indicated by the significance of their LCC size (empirical p-value < 0.05). Note that sizes of LCCs formed by seeds were compared to a reference distribution which was obtained by measuring the LCC size of random genes with degrees similar to the ones of the initial seed cluster. Thus, the observed local clustering cannot be attributed solely to the network's structural characteristics, but points to the network localization of the seed genes themselves, highlighting their relevance in marking the disease-related interactome community. (PDF)