Transcriptional Analysis Implicates Endoplasmic Reticulum Stress in Bovine Spongiform Encephalopathy

Bovine spongiform encephalopathy (BSE) is a fatal, transmissible, neurodegenerative disease of cattle. To date, the disease process is still poorly understood. In this study, brain tissue samples from animals naturally infected with BSE were analysed to identify differentially regulated genes using Affymetrix GeneChip Bovine Genome Arrays. A total of 230 genes were shown to be differentially regulated and many of these genes encode proteins involved in immune response, apoptosis, cell adhesion, stress response and transcription. Seventeen genes are associated with the endoplasmic reticulum (ER) and 10 of these 17 genes are involved in stress related responses including ER chaperones, Grp94 and Grp170. Western blotting analysis showed that another ER chaperone, Grp78, was up-regulated in BSE. Up-regulation of these three chaperones strongly suggests the presence of ER stress and the activation of the unfolded protein response (UPR) in BSE. The occurrence of ER stress was also supported by changes in gene expression for cytosolic proteins, such as the chaperone pair of Hsp70 and DnaJ. Many genes associated with the ubiquitin-proteasome pathway and the autophagy-lysosome system were differentially regulated, indicating that both pathways might be activated in response to ER stress. A model is presented to explain the mechanisms of prion neurotoxicity using these ER stress related responses. Clustering analysis showed that the differently regulated genes found from the naturally infected BSE cases could be used to predict the infectious status of the samples experimentally infected with BSE from the previous study and vice versa. Proof-of-principle gene expression biomarkers were found to represent BSE using 10 genes with 94% sensitivity and 87% specificity.


Introduction
Transmissible spongiform encephalopathies (TSEs), also termed prion diseases, are fatal, neurodegenerative diseases including Creutzfeldt-Jakob disease (CJD) in humans, scrapie in goats and sheep and bovine spongiform encephalopathy (BSE) in cattle [1,2,3,4]. The infectious agent of these diseases is thought to be an abnormally folded isoform (PrP Sc ) of the cellular prion protein (PrP C ) and it is further thought that the accumulation of the misfolded prion protein leads to disease [1]. PrP Sc is characterized by a high b-sheet content and resistance to protease treatment. In addition to the accumulation of PrP Sc , the pathological features of prion diseases in the brain of affected subjects include neuronal cell loss and vacuolation. Prion diseases have long incubation periods prior to the onset of clinical signs.
BSE was first discovered in 1986 [5] and became a major epidemic in the UK, peaking in 1992; to date more than 185,000 cases have been recorded. It is thought to be caused by contaminated meat and bone meal, a dietary supplement for cattle [6]. The BSE strain has also most probably crossed the species barrier to humans and has produced variant CJD [7]. The mean incubation period of BSE in cattle is estimated at about 5 years [3]. The clinical signs are: difficulties in locomotion and behavioural changes. The neuropathology of BSE is characterized by the lesions mainly found in the brain stem where vacuolar changes are found in neurons and the neuropil [5]. However, apoptosis plays a very limited role in neuronal loss in BSE [8].
In recent years atypical bovine spongiform encephalopathy has been identified [9,10]. In these cases the distribution of PrP Sc in the animal differs from that of BSE; there is less PrP Sc accumulation in the brain stem and the biochemical signature of PrP Sc is different.
The pathogenesis of BSE is still poorly understood. In a previous gene expression study using brain tissue samples from cattle experimentally infected with BSE, we have demonstrated that the largest number of differentially regulated genes is detected at 21 months post inoculation, suggesting that there are many pathogenic processes in the animal brain even prior to the detection of infectivity in the CNS of these orally dosed cattle [11]. Moreover, a set of differentially regulated genes could be used to predict the infectious status of preclinical samples.
To further understand the pathogenesis of BSE and to explore the possibility of using gene expression profiles as biomarkers, we analysed brainstem RNA samples from confirmed naturally infected cases of BSE (field cases) in cattle and from healthy controls.

Identification of differentially regulated genes in the BSE field case samples
The expression of genes in the brain of naturally infected BSE samples was compared with negative controls. In order to identify differentially regulated genes in BSE, the following stringent conditions were set with two filters: 2 fold change and one-way ANOVA with the p value being 0.05. 409 probe sets (a technical term that describes a transcript on the microarray) were identified as differentially regulated between the BSE infected (n = 14) and negative (n = 12) samples. After removal of duplications, 230 genes were identified and these genes are listed in Table 1 and the unannotated probe sets are listed in Table S1.
Only 18 genes (8%) were down-regulated and 212 (92%) genes were up-regulated after the repeated or un-annotated probe sets were removed (Table 1). Each step of filtering was re-examined to determine the number of up and down-regulated genes. The 2 fold change filter yielded 2138 probe sets: 792 (37%) of them were down-regulated in BSE field cases and 1346 (63%) probe sets were up-regulated. When the 2 fold change and 1-way ANOVA filters were combined, 409 probe sets were selected: 366 (89%) of them up-regulated and 43 (11%) probe sets down-regulated. Therefore, the up-regulated genes were increased in percentage after the ANOVA filter.
The largest functional group amongst the 230 identified genes was the genes involved in transport (39 genes), followed by the membrane protein group (25 genes), the metabolism group (20 genes) and the DNA and RNA binding group (19 genes; Table 1). The maximal increase was 5.75 fold for Myosin head domain containing 1 and the maximal decrease was 4.92 fold for OCIA domain containing 1.
Clustering analysis using the 409 probe sets showed that the samples were divided into two groups, Group A contained only negative control samples, while Group B contained all the BSE infected samples plus one negative control sample ( Figure 1). This analysis confirmed that the samples of BSE and controls were relatively homogeneous amongst themselves with regard to the genes defined as differentially regulated.
In response to ER stress, the unfolded protein response (UPR) is induced to restore cell function by reduction in newly translated proteins entering into the ER, by an increase in the capacity for protein folding [36]. If ER stress is prolonged, the UPR signaling pathways also initiate apoptosis [36]. In BSE, up-regulation of chaperones Grp94 and Grp170 suggests the induction of the UPR; while up-regulation of CerS2 indicates the inhibition of the UPR. To further explore the involvement of ER stress in the pathogenesis of BSE, Western blotting analysis on two more ER stress markers, Grp78 and Chop, was carried out. Grp78, is an ER chaperone and also known as an ER stress master regulator; while Chop is a transcription factor for induction of apoptosis, often upregulated in response to ER stress [36]. In BSE, the Grp78 protein was up-regulated ( Figure 2). Up-regulation of these ER chaperones: Grp78, 94 and 170 indicates the presence of ER stress and the activation of the UPR. The level of Chop was slightly decreased ( Figure 2) and this is consistent with the evidence that apoptosis plays a very limited role in BSE [8].
Using the gene expression profiles as a biomarker to represent BSE In our previous BSE time course study, 205 differentially regulated probe sets (corresponding to 114 genes) have been used to show that preclinical animals at 45 months post inoculation (mpi) cluster with cases positive for BSE and allowed the prediction that they are indeed preclinical and close to developing BSE [11]. The same 205 probes sets were used here in a clustering analysis to classify the disease status of the samples from the BSE field cases ( Figure 3). These samples fell into two main groups: Group A contained 11 positives and one negative, Group B contained the remaining 11 negatives and three positives. This analysis was therefore able to classify the samples according to infection status with 78.5% (11/14) sensitivity and 92% (11/12) specificity.
In a reverse analysis, the 409 probe sets identified in this study were used for clustering the samples from the time course study [11]. One group included the negatives, the samples from animals 6 mpi and 36 mpi and the other group contained the positives, and the samples from 21, 27 and 39 mpi animals ( Figure 4a). The clustering was similar to the one derived with the 205 probe sets from the time course study [11]. When these 409 probe sets were used to predict the status of the preclinical animals at 45 mpi in the time course study the clustering analysis grouped the individual samples into two groups: one with all the negatives (n = 3) and 6 mpi samples (n = 3) and the other with all the positives (n = 3) and 45 mpi samples (n = 2) (Figure 4b).
The analyses above indicate that either the genes from the time course study or the field case samples could be used to predict the    infection status. However, it would not be practical to apply all 409 or 205 probe sets as biomarkers to represent BSE. A group of 10 genes were sought to represent BSE from these 230 genes listed in Table 1. Initially, the search was carried out using genes associated with prion diseases (10 genes), ER stress (10 genes), the largest fold changes (10 genes) or the smallest p values (10 genes) separately but the sensitivity and specificity of prediction were low. When these 40 genes were combined and10 genes were selected from them by comparing the expression levels of individual samples from both this study (clinical BSE, n = 14; control, n = 12) and the time course study (clinical BSE, n = 3; control, n = 3), only two groups were produced (Figure 5a) These 10 genes above were then used to classify the preclinical samples from the time course study with clustering analysis. The clustering analysis produced two groups: three negatives and three 6 mpi samples being one group and three positives and 45 mpi samples being the other with 100% (5/5) sensitivity and 100% (6/ 6) specificity ( Figure 5b). Therefore, the results of these analyses suggest that these 10 genes might be used to represent the patterns of BSE gene expression at the terminal stages of BSE.

Discussion
In this study, 230 genes were found to be differentially regulated between BSE field cases and controls (Table 1). These genes belong to many functional groups from apoptosis to transport. Seventeen genes were associated with the ER and 10 of them may be involved in stress related situations, especially up-regulation of ER chaperones Grp94 and Grp170 as they are ER stress markers. Since ER stress triggers the UPR [37,38,39], the level of protein expression of Grp78, another ER stress marker, was increased in BSE. Up-regulation of Grp78, Grp94 and Grp170 is induced by ER stress response transcription factors XBP1 and ATF6 as all three of them have an ER stress response element (ERSE) in their regulatory regions [36]. These analyses suggest the presence of ER stress and the activation of the UPR in the disease process of BSE. This is in agreement with increasing evidence of the involvement of ER stress in prion diseases [40,41,42]. In this study only changes in gene and protein expression of these chaperones were measured to indicate activation of the UPR. There are other methods to measure the induction of the UPR as many proteins are activated or inactivated through phosphorylation cascade in the UPR  signalling pathways. For example, the release of Grp78 bound to PERK triggers autophosphorylation of PERK which in turn phosphorylates elf2a to attenuate protein translation [36].
To cope with accumulation of misfolded proteins, ER stress induces ER associated protein degradation I (ERAD I, ubiquitin/ proteasome) [43] and ERAD II (autophagy/lysosome) [44], possibly through the UPR. ERAD I is closely linked to the ER quality control system [45] as unfolded or misfolded proteins are targeted for degradation after the failed attempt of folding by ER chaperones. In BSE, ubiquitin-activating enzyme E1 (UBE1) and three E3 ligases: WW domain containing E3 ubiquitin protein ligase 2 [46], ariadne homolog 2 [47] and ubiquitin carboxylterminal esterase L1 [48] were found to be up-regulated (Table 1). Recently, the E3 ligase HECTD2 has been identified as genetically associated with vCJD and kuru [49].
ERAD II is also known as autophagy. It is a pathway of selfdegradation of cellular components in which autophagosomes sequester organelles or protein aggregates and fuse with lysosomes for degradation. When the production of misfolded proteins exceeds the capacity of ER chaperones and ERAD I, misfolded and aggregated proteins are targeted by the aggresome-autophagy pathway [50]. In BSE, up-regulation of several genes (Table 1) suggests that this pathway might be induced. In the lysosome, both cathepsin B and D (lysosomal hydroases) were up-regulated [44]. On the membrane of the lysosome, the increased levels of lysosomal-associated membrane protein 2 (LAMP2) suggest autophagy initiation [51]. In the cytosol, there were also several up-regulated genes related to ERAD II, such as ubiquitin carboxyl-terminal esterase L1 (aggresome initiation in proteasome inhibition) [48], sphingomyelin phosphodiesterase 1, acid lysosomal (SMPD1; autophagy promotion) [52] and vimentin (cytoskeleton) [14] (Table 1). This association between ERAD II and BSE has been shown in both mice and cattle [53,54].
In this study, the analyses suggest that ER stress might be involved in BSE pathogenesis and that the UPR, ERAD I and II might all be activated in a concerted effort to rid the cell of harmful PrP Sc . The question, therefore, is how much these ER related pathogenic events contribute to fatal prion diseases in general. When the GPI anchor of the PrP protein is removed, the transgenic mice infected with scrapie, also a prion disease, can survive up to 400-600 days post infection (dpi) without clinical scrapie, while the wild type controls develop clinical signs within 140-160 dpi [55]. Some animals with this anchorless PrP have up to 40% more PrP Sc than clinically sick controls. The results indicate that infectivity (PrP Sc accumulation) and toxicity can be uncoupled. One model to explain it is intra neuronal generation of a toxic intermediate [32]. Here we offer another explanation of prion neurotoxicity using ER stress. The reason for PrP Sc accumulation in the ER is because the ER quality control system senses the misfolded forms of PrP and ER chaperones retain them in the ER for folding or degradation by ERAD I. PrP Sc is protease resistant so that the rate of removing the misfolded protein is slow; while more and more PrP C converts to PrP Sc . Eventually, PrP Sc Figure 3. Clustering analysis of sample status in the BSE field case study. The analysis was performed by GeneSpring using 205 differential regulated probe sets generated from the time course study [11]. The similarity was measured using the spearman correlation. N: BSE negative controls; P: clinical BSE samples. doi:10.1371/journal.pone.0014207.g003 Figure 4. Clustering analysis of sample status in the BSE time course study. The analysis was performed by GeneSpring using 409 differential regulated probe sets generated from this study and the samples were from the BSE time course study [11]. The similarity was measured using the spearman correlation. Neg: accumulation causes ER stress and the subsequent activation of the UPR and ERAD I and II. Prolonged ER stress leads to cell death [45,56]. Hence, ER stress related responses might be the major source of prion toxicity. What happens when misfolded PrP Sc bypasses the ER quality control? There are lines of evidence that the anchorless prion protein is not detected by the ER quality control system [57,58]. As the anchorless PrP Sc can pass the ER efficiently, there is no toxicity to cause clinical scrapie. Since the cell is not under ER stress, the ERAD pathways are not activated. As a result, more PrP Sc accumulates in the brain of transgenic mice with anchorless PrP than in the brain of the wild type controls.
Both the current field case study and the previous time course study were carried out with brainstem tissues infected with BSE [11]. Although these two sets of samples differed in age, in infectious dose and in stages of disease development, many differentially regulated genes were expected to be shared between these two studies. Nonetheless, when the two gene lists were compared, there were only two genes overlapping. However, the profiles generated from one study could be used to predict the sample status of the other study as biomarkers, suggesting that there were some underlying links between these two gene lists (Figures 3 and 4). One possible explanation is that there are more differentially regulated genes than those identified by the analytical method. In order to define a gene list that is relevant with a condition or a disease within a study, the p value is often set at 0.05 or less. However, by doing so, much of the coverage is lost and many differentially regulated genes are not considered. In order to make the list more manageable, an additional 2 fold change filter was introduced to reduce the number of probe sets to 409. If the fold change filter had not been introduced and the p value had been set at 0.1, the number of probe sets would have been 1604. By definition, only 160 of them were selected randomly and the rest of 1446 probe sets should be truly differentially regulated. The remaining 1037 (1446-409) probe sets were not analyzed. Figure 6 provides a simple graphical model for this situation. The small inner circles (stringent settings) overlap only marginally. If all differentially regulated genes had been considered (large circles), there would have been many genes shared by these two studies and that is the most likely reason why the profiles from one study could be used to predict sample status from the other study. In recent years, there have been many publications on gene expression analyses of prion diseases. It is a surprise that relatively few differentially regulated genes are shared between these studies [11,12,13,16,59,60]. However, the explanation above for the BSE studies may also apply to gene expression studies of prion diseases in general.
Considerable efforts have been made to find biomarkers for the prion diseases, especially in the early stage of the incubation period. To date, the detection of PrP res is still the only reliable method. There are some reported potential biomarkers for the disease such as 14-3-3 protein [61], galectin-3 [62], SCRG1 [63], clusterin [64] and cystatin C [65]. However, none of them has been developed for routine diagnosis. One of the reasons may be the natural variation for the single marker within a population. Clustering analysis suggested that a prediction could be made by comparing the gene expression profiles of a sample with those of known BSE positive and negative samples. The analysis also showed a proof of principle that a prediction for a given sample could be made with high sensitivity (94%) and specificity (87%) using just 10 genes as biomarkers although the tissues used in this study were from the brainstem which may not be suitable for diagnose. These ten gene markers might represent the diseased state better than any single markers as they might allow some variations in expression. In Huntington's disease, gene expression profiling of blood reveals a subset of 12 up-regulated mRNAs Figure 5. Clustering analysis of possible biomarkers for BSE. The analysis was performed by GeneSpring using 10 differential regulated genes generated from this study. The similarity was measured using the change correlation with value 1 for separation ratio and value 0.001 for minimum distance in merge similar branches. (a), All samples. N: BSE negative controls and P: clinical BSE samples in this study. Pos: clinical BSE and Neg: the negatives from the BSE time course study [11]. which have been shown to be able to distinguish controls, presymptomatic Huntington's disease gene carriers and symptomatic Huntington's disease patients [66].
In conclusion, gene expression analysis suggests that BSE infection caused ER stress and the UPR, ERAD I and II might be induced in response to ER stress. Clustering analysis showed that the differentially regulated genes could be used to predict infection status. Ten genes were selected to represent gene expression state in BSE, which might eventually be used as biomarkers.

Tissue samples
Brainstem tissues from 100 confirmed cases of BSE in cattle were supplied by the TSE archive at the Veterinary Laboratories Agency, UK. The animals were females, between 4 and 10 years old that had been diagnosed clinically and killed on farm. The major breed was Holstein/Friesian and other breeds were: Limousin Cross, Guernsey, Hereford Cross and Brown Swiss. The negative controls (100 brainstem tissue samples) were from LGC Forensics (Queens Road, Teddington, Middlesex, TW11 0LY, UK) and were comparable in breed, sex and age with the naturally infected BSE samples. Since all samples were from the Archives, approval from the Ethics Committee was not necessary.

Microarrays analysis
The preparations of samples and reagents were carried out according to the Affymetrix GeneChip Expression Analysis manual and as described in the previous study [11]. The RNA samples were resolved by 1% agarose gels and selected according to the integrity of ribosomal RNA bands. Since the tissues used in this study were from cattle naturally infected with BSE (field cases), the quality of RNA was generally poor. From 100 cases each, the best quality RNA samples, 12 controls and 14 BSE infected, were selected for microarray analysis with Affymetrix GeneChip Bovine Genome Arrays. The raw data were first imported into the Affymetrix GeneChip operating software version 1.4. All array data were MIAME compliant and the raw data were deposited in ArrayExpress with the accession number: E-MTAB-302. After initial analysis, the pivot formatted data were further analysed with the GeneSpring version 7 software (Silicon Genetics). The data were normalized in three steps: 1. Data transformation set measurements less than 0.01 to 0.01; 2. Each measurement was divided by the 50.0th percentile of all measurements in that sample; 3. Each measurement for each gene in test samples was divided by the median of that gene's measurements in the corresponding control samples. The value for each gene was divided by the median of its measurements in all samples. If the median of the raw values was below 10 then each measurement for that gene was divided by 10. If the numerator was above 10, the measurement was discarded. These steps were the default settings for the GeneSpring package.
Two filters were used to find differently regulated genes: 2 fold change and the one way ANOVA statistical analysis with the parameters of 0.05 for p-value cutoff, multiple testing correction and Student-Newman-Keuls for the post hoc tests, without assume variances equal for the parametric test.

Western blotting
Cell-free extracts (60 mg protein) were loaded on 12% 1-D SDS PAGE (Invitrogen) and resolved proteins from several mini-gels were transferred to the same PVDF membrane (Millipore) so that one set of samples was used to monitor protein loading using b-Actin. The blots were immuno-stained with mouse monoclonal anti-b-Actin IgG (Santa Cruz Biotech), rabbit polyclonal anti-Grp78 (US Biological) and rabbit polyclonal anti-Chop (BioLegend). The protein bands were visualized by using secondary antibodies, alkaline phosphatase conjugated IgGs (anti-mouse, Santa Cruz Biotech; anti-rabbit, Sigma) and the ECL developer kit (Amersham). The images were captured by Fluor-S Multi-Imager (Bio-Rad) and the protein bands were quantified by the Quantity One software (Bio-Rad).

Quantitative PCR
The RNA samples were treated with the DNA free TM kit (Ambion) for 1 h at 37uC to remove any trace of DNA. The treated RNA was then used as a template for cDNA synthesis with the TaqMan reverse transcription kit (Applied Biosystems). The real time PCR was carried out by denaturing at 95uC for 15 s, annealing at 50uC for 2 min and extension at 60uC for 1 min for 40 cycles using an ABI Prism 7700 Sequencing Detector. The GAPDH gene was used as an internal control to normalize the expression levels of target mRNA. The primer sets were chosen by the Primer Express 1.5 for TaqMan software. The sequences of the primer sets were as following: for CD47, 59