Brain transcriptomes of zebrafish and mouse Alzheimer's disease knock-in models imply early disrupted energy metabolism

Alzheimer's disease (AD) is a complex and highly heterogenous neurodegenerative disease, defined by the presence of intracellular neurofibrillary tangles (primarily consisting of hyperphosphorylated tau proteins), and extracellular plaques mostly consisting of a small peptide, amyloid β (Aβ) (Jack et al., 2018). The pathological basis of AD has been the subject of research for over 100 years (Alzheimer, 1906). Nevertheless, most treatments tested in clinical trials have shown limited therapeutic benefit.

AD has a strong genetic basis (reviewed by Sims et al., 2020). In some rare cases, early-onset familial forms of AD (EOFAD, occurring before 65 years of age) arise due to dominant mutations in one of four genes: presenilin 1 (PSEN1), presenilin 2 (PSEN2), amyloid β precursor protein (APP) and sortilin-related receptor 1 (SORL1) (reviewed by Barthelson et al., 2020a; Bertram and Tanzi, 2012; Ayodele et al., 2021). However, most AD cases are sporadic, showing symptom onset after the arbitrarily defined threshold of 65 years (late-onset sporadic AD, LOAD). Genetic variants at many loci have been associated with increased risk of LOAD (Jansen et al., 2019; Kunkle et al., 2019; Lambert et al., 2013). The most potent variant is the ε₄ allele of apolipoprotein e (APOE) (Farrer et al., 1997), which has been described as ‘semi-dominant’ (Genin et al., 2011).

An understanding of the early cellular stresses on the brain that eventually lead to AD is necessary to advance the development of preventative treatments. This is difficult to achieve through studying living humans, as EOFAD mutations are rare, and access to young presymptomatic brains is limited. Nevertheless, imaging studies have implicated structural and functional changes to the brain long before diagnosis of AD (Iturria-Medina et al., 2016; Quiroz et al., 2015). Brain imaging cannot provide detailed molecular information about these changes. Transcriptome analysis is, currently, the strategy that can provide the highest resolution molecular description of cells and tissues. However, transcriptome analyses of ante-mortem brains carrying EOFAD mutations can only be performed using brain tissue from animal models.

Our group has exploited the zebrafish to generate a collection of knock-in models of EOFAD-like mutations in order to analyse their young brain transcriptomes (Barthelson et al., 2020b, 2021a, 2021b, 2021c; Dong et al., 2021; Hin et al., 2020a, 2021; Jiang et al., 2020; Newman et al., 2019). Our experimental philosophy has been to replicate, as closely as possible, the single heterozygous mutation state of EOFAD in humans, thereby avoiding possibly misleading assumptions regarding the molecular mechanism(s) underlying the disease. Our overall goal has been to compare a broad range of EOFAD-like mutations in a number of EOFAD genes to define their shared pathological effects in young adult brains where the long progression to AD begins. To assist in this definition (by exclusion), we also created non-EOFAD-like mutations in the same genes as negative controls, i.e. frameshift mutations in the presenilin genes that do not cause EOFAD (reviewed by Jayne et al., 2016, the ‘reading frame preservation rule’). The presentation of EOFAD and LOAD as similar diseases (reviewed by Blennow et al., 2006; Masters et al., 2015) implies similarity, to some degree, at the cellular and molecular levels. Therefore, despite differences in the genetic variants that promote these two diseases, understanding the molecular effects of heterozygosity for EOFAD mutations may give insight into the molecular changes underpinning LOAD.

Here, we summarise our findings of brain transcriptome analyses of EOFAD-like mutations in the zebrafish orthologues of genes implicated in EOFAD: psen1, psen2 and sorl1. EOFAD mutations also exist in APP. However, zebrafish express two APP ‘co-orthologous’ genes, appa and appb, complicating analysis of single heterozygous mutations. Therefore, we re-analysed the best available publicly accessible brain transcriptomic data from a knock-in model of APP mutations: the App^NL-G-F mouse. Finally, we compared whether the brain transcriptome changes occurring due to single heterozygous EOFAD-like mutations in zebrafish are similar to the changes occurring due to the strongest genetic risk factor for LOAD, the ε₄ allele of APOE, using publicly available brain transcriptome data from a humanised APOE targeted-replacement mouse model (APOE-TR) (Sullivan et al., 1997). We identify changes to energy metabolism as the earliest detectable cellular stress due to AD mutations, and demonstrate that knock-in zebrafish models are valuable tools for studying the earliest molecular pathological events in this disease.

We first collated our findings from our zebrafish models of EOFAD-like mutations in psen1 (Barthelson et al., 2021a; Hin et al., 2020, 2021; Newman et al., 2019), psen2 (Barthelson et al., 2021b) and sorl1 (Barthelson et al., 2020b, 2021c). An advantage of using zebrafish for RNA-seq analyses is minimisation of genetic and environmental noise through breeding strategies, such as that shown in Fig. 1A. Large families of synchronous siblings can consist of heterozygous mutant and wild-type genotypes, allowing direct comparisons of the effects of each mutation. So far, we have performed six brain transcriptomic analyses based on various breeding strategies (summarised in Table 1 and Figs S1-S6). The detailed analyses can be found in the publications cited above. However, the outcomes are summarised below and in Fig. 1.

In our previously published analyses, we found that heterozygosity for most of our EOFAD-like mutations does not result in many differentially expressed genes in young adult brains (as would be expected for modelling a disease that becomes overt in middle age) (Barthelson et al., 2020b, 2021a, 2021b, 2021c; Newman et al., 2019). Therefore, we performed gene set enrichment analyses (GSEA) to predict which cellular processes were affected by each of the mutations in each experiment. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) gene sets to determine whether changes to gene expression were observed in any of 186 biological pathways/processes. Additionally, we recently proposed that neuronal iron dyshomeostasis may be an effect-in-common of EOFAD mutations in the context of AD pathogenesis (Lumsden et al., 2018). Therefore, we used our recently defined iron responsive element (IRE) gene sets (Hin et al., 2021) to test for evidence of iron dyshomeostasis. Biological processes found to be affected in at least two different zebrafish mutants are shown in Fig. 1B (the statistical significance of all KEGG and IRE gene sets in each mutant can be found in Table S2). The one gene set consistently altered by all of the EOFAD-like mutations, but not by the non-EOFAD-like mutations examined, is the KEGG_OXIDATIVE_PHOSPHORYLATION gene set (Fig. 1C), supporting that changes to mitochondrial function are an early cellular stress in EOFAD. The KEGG_OXIDATIVE_PHOSPHORYLATION gene set is also affected by heterozygosity for the K97fs mutation of psen1. K97fs is a frameshift mutation and so does not follow the ‘reading frame preservation rule’ (Jayne et al., 2016) of presenilin EOFAD mutations. However, the truncated protein encoded by K97fs resembles a hypoxia-induced isoform of human PSEN2, denoted PS2V, which shows increased expression in LOAD brains (Sato et al., 1999, and see Moussavi Nik et al., 2015, for additional explanation). Therefore, K97fs is still an AD-relevant mutation.

Genes encoding the components of ribosomal subunits, as defined by the gene set KEGG_RIBOSOME, were affected by all the EOFAD-like mutations but also by non-EOFAD-like mutations in psen1 and psen2 (Fig. 1D). Evidence for iron dyshomeostasis was also observed for the relatively severe EOFAD-like mutation psen1^Q96_K97del/+ (under both normoxia and acute hypoxia conditions) and in transheterozygous sorl1 mutants (i.e. with complete loss of wild type sorl1), as shown by significant changes to the expression of genes possessing IRE(s) in the 3′ untranslated regions (UTRs) of their encoded mRNAs (ire_hq and ire_all).

EOFAD is also caused by mutations of the gene APP. Modelling of APP mutations in zebrafish is complicated by duplication of the APP-orthologous gene in this organism. However, brain transcriptome data are available for a knock-in mouse model of EOFAD mutations in APP: the App^NL-G-F mouse model (Castillo et al., 2017). In this model, the murine App sequence is modified to carry humanised DNA sequences in the Aβ region, as well as the Swedish, Beyreuther/Iberian and Arctic EOFAD mutations (Saito et al., 2014). Although these mice do not closely reflect the genetic state of heterozygous human carriers of EOFAD mutations of APP (as the mice possess a total of six mutations within their modified App allele and are usually analysed as homozygotes), they should, at least, not generate artefactual patterns of gene expression change due to overexpression of transgenes (Saito et al., 2016). Castillo et al. (2017) performed brain transcriptomic profiling via microarrays of the brain cortices of male homozygous App^NL-G-F mice relative to wild-type mice at 12 months of age, as well as a transgenic mouse model of AD, 3xTg-AD mice (Oddo et al., 2003) relative to non-Tg mice (Castillo et al., 2017). All mice used in the study were maintained as inbred lines. However, there is no information on whether any of the mice analysed were littermates. It is highly unlikely that the mice used in each comparison between mutant individuals and their wild-type counterparts all arose from the same litter, because obtaining three homozygous and three wild-type male mice in a single litter arising from an incrossing of heterozygous mutant mice (expected to produce a wild type:heterozygote:homozygote Mendelian genotype ratio of 1:2:1) would be a rare event, as litters of mice generally consist of five to ten pups. Therefore, additional variation was introduced into the analysis through the use of mice from different litters and this is likely confounding with genotype. This is important to note, as the results presented here were generated under the assumption that any effects of litter-of-origin are negligible. Our re-analysis of the microarray dataset of Castillo et al. (2017) aimed to address the following questions: (1) are the KEGG gene sets affected in the male homozygous App^NL-G-F mice similar to those affected in EOFAD-like zebrafish?; and (2) is there evidence for iron dyshomeostasis in the brains of these mice? Initially, we attempted to replicate the results of Castillo et al. (2017) using the Affymetrix Transcriptome Analysis Console software. However, we were unable to find sufficient information to achieve this. Therefore, we analysed the microarray dataset in a reproducible manner following a recommended microarray analysis workflow (Carvalho and Irizarry, 2010).

After pre-processing of the raw intensities (Fig. S7), we performed principal component analysis (PCA) to explore the overall similarity between samples (Fig. 2A). Samples separated across PC1 by genotype, suggesting that the homozygous genotypes in this study result in distinct transcriptome states. Notably, the App^{NL-G-F/NL-G-F} samples and their corresponding App^+/+ control samples appeared to separate to a greater extent across PC1 than the 3xTg samples and their corresponding non-Tg wild-type control samples. Additionally, a differential gene expression analysis revealed 158 and 126 genes to be differentially expressed in App^{NL-G-F/NL-G-F} and 3xTg mice, respectively (relative to their corresponding controls, Fig. S8). This suggests that the disturbance to the cortex transcriptome in App^{NL-G-F/NL-G-F} mice is greater than in 3xTg mice.

We did not observe alteration of any similar gene sets between the App^NL-G-F mice and our EOFAD-like zebrafish (Figs S9-S12). However, any similarities may well have been masked by the overwhelming effects of greater age, variable environment (mouse litter-of-origin) and the effects of their six App mutations on brain cortex cell-type proportions and inflammatory processes. The most statistically significantly affected cellular process in 12-month-old App^NL-G-F mice was lysosomal function, as represented by the KEGG_LYSOSOME gene set (Fig. 2D). Additionally, a plethora of inflammatory gene sets were also affected, with changes in the relative proportions of glial cells, particularly microglia, contributing to the appearance of increased levels of these gene transcripts in the bulk cortex RNA analysed (Fig. 2B,C; Fig. S13). Note that changes to cell-type proportions are not observed in our zebrafish models of EOFAD (Barthelson et al., 2020b, 2021b, 2021c) (see Fig. S14 for two examples).

A puzzling observation from the genome-wide association studies (GWAS) of LOAD, is that none of the risk variants identified fall within the EOFAD genes PSEN1, PSEN2 or APP (Jansen et al., 2019; Kunkle et al., 2019; Lambert et al., 2013). This has led to speculation that EOFAD and LOAD may be distinct diseases despite their histopathological and cognitive similarities (reviewed by DeTure and Dickson, 2019; Tellechea et al., 2018). Only one gene identified by GWAS of LOAD, SORL1, is suspected to harbour mutations causative of EOFAD. Mutations in SORL1 cause AD with ages of onset typically later than many mutations in PSEN1 or APP (or maybe incompletely penetrant) (Pottier et al., 2012; Thonberg et al., 2017). Nevertheless, as shown in Fig. 1B, we identified changes in the KEGG gene set for oxidative phosphorylation in young adult zebrafish heterozygous for EOFAD-like mutations in sorl1, as well as in zebrafish modelling overexpression of the PS2V isoform that is upregulated in LOAD (K97 fs).

The strongest and most common genetic risk factor for LOAD is the ε₄ allele of the gene APOE (Corder et al., 1993; Genin et al., 2011; Jansen et al., 2019; Kunkle et al., 2019; Lambert et al., 2013; Saunders et al., 1993). Like APP, the APOE orthologous gene in zebrafish is refractory to analysis due to duplication. Therefore, to compare our zebrafish mutant data to early brain transcriptome changes caused by the ε₄ allele of APOE, we analysed data from a set of human gene-targeted replacement mouse models, APOE-TR (Sullivan et al., 1997). These mouse models transcribe human APOE alleles from the endogenous murine Apoe promotor: the predominant human allele ε₃; the rare AD-protective ε₂ allele; and the AD-risk allele ε₄. Zhao et al. (2020) performed a comprehensive brain transcriptome profiling experiment across aging in both male and female mice to assess the effect of homozygosity for the ε₂ or ε₄ alleles relative to the risk-neutral ε₃ allele. In that analysis, pairwise comparisons between the ε₂ or ε₄ alleles relative to the ε₃ allele (e.g. Fig. 3A) were not conducted at each age and sex. Only genes/pathways that were influenced overall by APOE genotype, age, sex and interactions between these factors were reported. As our aim is to identify the early changes occurring due to AD-related mutations, we re-analysed only the 3-month brain samples from the Zhao et al. (2020) dataset (i.e. omitting the samples from 12- and 24-month-old mice) to investigate which processes are affected by homozygosity for the ε₂ or ε₄ alleles relative to the ε₃ allele. Hereafter, we refer to these homozygous mice as ‘APOE2’, ‘APOE3’ and ‘APOE4’.

After pre-processing of the APOE-TR RNA-seq data (Figs S15, S16), we performed PCA to visualise the overall similarity between APOE-TR brain transcriptomes. The plot of PC1 against PC2 revealed that samples separated into two distinct clusters of sex across PC2 (Fig. 3B). This suggests that the effect of sex on the murine brain transcriptome is substantial and cannot be ignored in the differential gene expression analysis. Among the male samples, APOE4 samples formed a cluster distinct from the APOE2 and APOE3 samples, suggesting that the APOE4 genotype has a distinct effect on the transcriptome compared to APOE2 relative to APOE3 in males. This was not observed to the same extent in the female samples. However, the male APOE4 and APOE3 samples appeared to have been taken from distinct litters, as implied from the date of birth of each sample (Fig. 3B). This confounds the effect of genotype and complicates interpretation of whether any effects observed in a pairwise comparison between male APOE3 and APOE4 mice are due to APOE genotype or litter-of-origin (or, most likely, both). Indeed, χ² tests for independence revealed that there is a highly significant dependence of APOE genotype and litter across the entire 3-month-old dataset (χ²=82.7, d.f.=20, P=1.4×10⁻⁹), as well as within only male samples (χ²=43.1, d.f.=14, P=8.2×10⁻⁵) and within only female samples (χ²=39.3, d.f.=14, P=3.0×10⁻⁴). Some litters did not contain sufficient mice to remove the effect (i.e. some coefficients could not be estimated during the generalised linear model fitting procedure due to the design matrix not having full rank). Therefore, we continued the analysis under the assumption of a negligible effect of litter.

To determine which genes were dysregulated in APOE4 mice and APOE2 mice relative to APOE3 mice, we performed a differential gene expression analysis using edgeR (McCarthy et al., 2012; Robinson et al., 2009). Many genes were found to be differentially expressed in each comparison, particularly in male APOE4 mice. Additionally, the biases noted by Zhao et al. (2020) in the original analysis for increased GC content and longer transcript length among differentially expressed genes was also apparent in our analysis (Fig. S17). Therefore, we corrected for these observed biases using conditional quantile normalisation (CQN) (Hansen et al., 2012). After CQN, many genes were identified as differentially expressed in each comparison, and the percentage of GC and gene length biases were decreased (Fig. S18).

We next performed enrichment analysis of the KEGG (Kanehisa and Goto, 2000) and IRE (Hin et al., 2021) gene sets to determine whether changes are observed in APOE4 mice similar to those in our EOFAD-like model zebrafish (Figs S19-S23 and Table S1). We found statistical evidence for significant changes in the expression of oxidative phosphorylation (Fig. 3C) and ribosome (Fig. 3D) gene sets in mice homozygous for the humanised ε₄APOE allele, consistent with our zebrafish models of EOFAD (although different genes appear to drive the statistical enrichment of the KEGG_OXIDATIVE_PHOSPHORYLATION gene set in the two organisms, Fig. S25). Interestingly, we observed highly statistically significant changes in the gene set ire3_all only in APOE4 male mice, reminiscent of similar signals in some of the young adult EOFAD-related zebrafish (Barthelson et al., 2020b; Hin et al., 2021; Fig. 1B; Fig. S22C) and supporting the existence of iron dyshomeostasis. These effects were not observed for the AD-protective ε₂ allele (Fig. 3; Tables S2, S3). However, the effects of APOE genotype were highly dependent on the litter-of-origin of the samples, and changes to cell-type proportions were observed in the male APOE4 mice (Fig. S24). Therefore, future replication of this analysis with better-controlled transcriptome data is desirable to confirm that the effects observed are due to the APOE genotype.

Energy production is the most fundamental of cellular activities. Life cannot be sustained without energy, and all other cellular activities depend upon it. The human brain, in particular, has very high energy demands and consumes the majority of the glucose of the body when at rest (reviewed by Zierler, 1999). Within the brain, the majority of energy use is to maintain the Na⁺-K⁺ membrane potential of neurons (Attwell and Laughlin, 2001), and neurons are assisted in meeting these energy demands by support primarily from astrocytes [e.g. via the astrocyte-neuron lactate shuttle (Pellerin and Magistretti, 1994)]. All cells allocate considerable portions of their energy budgets to protein synthesis to maintain their structure and activity (Buttgereit and Brand, 1995). Energy is also required to maintain the low pH and high Ca²⁺ concentration of the lysosome (Christensen et al., 2002), the organelle which mediates the uptake and recycling (autophagy) of cellular structural constituents (e.g. the amino acids for protein synthesis) (reviewed by Yim and Mizushima, 2020). Lysosomes are important for the uptake and recycling of ferrous iron (Yambire et al., 2019), which is essential for oxidative phosphorylation by mitochondria (Oexle et al., 1999). On the lysosomal membrane, mammalian target of rapamycin kinase (mTOR) complexes sense nutrient and energy status to regulate protein synthesis, autophagy and mitochondrial activity (reviewed by Lim and Zoncu, 2016).

The EOFAD genes PSEN1, PSEN2, APP and SORL1 all encode proteins expressed within the endolysosomal pathway of cells (Andersen et al., 2005; Kawai et al., 1992; Pasternak et al., 2003; Sannerud et al., 2016) and within the mitochondrial associated membranes (MAMs) of the endoplasmic reticulum (Area-Gomez et al., 2009; Lim, 2015). MAMs are responsible for the regulation of ATP production [through Ca²⁺ signalling (Duchen, 1992)], oxidative protein folding (reviewed by Simmen et al., 2010) and the initiation of autophagy (Hamasaki et al., 2013). Interestingly, like EOFAD mutant forms of PSEN1 (Lee et al., 2010) and the C99 fragment of APP (Jiang et al., 2019), the ε4 allele of APOE has been shown to affect both lysosomal pH (Prasad and Rao, 2018) and the MAM (Tambini et al., 2016). Our analyses of young adult brain transcriptomes in zebrafish have found that five EOFAD-like mutations in a total of three EOFAD gene orthologues (psen1, psen2 and sorl1) all cause statistically significant effects on the expression of genes involved in oxidative phosphorylation, whereas non AD-related mutations in psen1 and psen2 do not. Therefore, effects on oxidative phosphorylation are a common, early ‘signature’ of EOFAD. Intriguingly, we previously observed downregulation of the oxidative phosphorylation genes due to heterozygosity for the psen1^Q96_K97del mutation in whole-zebrafish larvae at 7 days post fertilisation (dpf) (Dong et al., 2021), suggesting changes to mitochondrial function are a very early cellular stress in EOFAD pathogenesis. Additionally, we observed that the ‘semi-dominant’ ε4 LOAD risk allele (Genin et al., 2011), like EOFAD mutations, also affects the expression of genes involved in oxidative phosphorylation and ribosome function in young adult brains. Thus, changes in oxidative phosphorylation (and so energy production) appear to be a common early disturbance associated with both early- and late-onset forms of AD.

The majority of the heterozygous EOFAD-like mutations we have studied in zebrafish cause an overall downregulation of the oxidative phosphorylation gene set in young adult brains relative to wild-type brains (Fig. 1C). Only heterozygosity for the T141_L142delinsMISLISV (reading frame-preserving) mutation of psen2 has been seen to give overall upregulation of these genes. Another complex yet probable EOFAD-like mutation in psen2 we have studied in zebrafish, psen2^S4ter (which likely produces Psen2 proteins lacking N-terminal sequences), also showed strong overall upregulation of the oxidative phosphorylation gene set (Jiang et al., 2020); however, that dataset contains technical artefacts that complicate interpretation, so it was not included in the current analysis. Transheterozygosity for mutations in sorl1 also results in overall upregulation of oxidative phosphorylation genes. We are uncertain as to why this variability in effects on the oxidative phosphorylation gene set occurs. However, for the Presenilins, the single most consistent characteristic of the hundreds of known EOFAD mutations is that they maintain the ability of the genes to produce at least one transcript isoform with the original reading frame [the ‘reading frame preservation rule’ (Jayne et al., 2016)]. This strongly supports that all these mutations act via a dominant gain-of-function molecular mechanism (to interfere with a normal cellular function). Alternatively, it may be that the disruption of this gene set that is consistently observed is a product of both genotype and environmental factors, i.e. the mutant fish may be more or less responsive to environmental variation, such as changes in water quality, microbiome, handling, etc. Additionally, we note that it can be misleading to infer the direction of change in a particular cell activity, such as oxidative phosphorylation, based on the majority behaviour of a (somewhat arbitrarily) defined set of genes. Obviously, actual measurement of, for example, respiratory rates in the zebrafish mutant brains would be needed to establish, with certainty, how the mutations are affecting oxidative phosphorylation. Note, however, that the subtlety of the gene regulatory effects we have observed in the fish models means that discernment of physiological oxygen consumption differences between mutant fish and their wild-type siblings may be challenging. (Simultaneous measurement of differences in the expression levels of the ∼100 genes in the oxidative phosphorylation gene set gives great statistical sensitivity for detection of subtle differences.)

Importantly, male mice homozygous for LOAD risk allele APOE ε₄ showed altered expression of the oxidative phosphorylation gene set, and female APOE4 mice showed a similar trend that did not reach the threshold for statistical significance. In fact, male APOE4 mice appeared to have more disrupted brain transcriptomes than female APOE4 mice (Fig. 3B), including alteration in the abundance of transcripts with IREs in their 3′ UTRs, supporting the possibility of early brain iron dyshomeostasis (as was previously observed in psen^Q96_K97del/+ fish and those lacking wild-type sorl1 function). This was unexpected, as human females are more susceptible to the effects of APOE ε₄ than males (Farrer et al., 1997; Wang et al., 2020). However, other sources of variation (i.e. litter-of-origin, changes to cell-type proportions and possibly transcript length) may be masking the true effects of APOE genotype in these mice. Reassuringly, a recent single-cell RNA-seq analysis (in which discrepancies due to differences in cell-type proportions are overcome) of female APOE4 mice showed that expression of oxidative phosphorylation genes is decreased at 12 months of age, particularly in astrocytes. Further dissection of the metabolic phenotype of female APOE4 mice revealed a shift away from oxidative phosphorylation and towards glycolysis for ATP production (Farmer et al., 2021). Remarkably, this result was consistent with observations in human ε₄-carrying females, who show lower energy expenditure, decreased oxygen consumption and alterations to their plasma metabolomes indicative of increased glycolysis (Farmer et al., 2021).

The changes to gene expression in the oxidative phosphorylation pathway in both EOFAD-like zebrafish and APOE4 mice demonstrates the similarity, at the molecular level, between the cellular effects of genetic variants causing EOFAD and the most significant variant promoting LOAD. The differences in disease onset age between EOFAD and LOAD may be due to the severity of the effects on energy metabolism of the different genetic variants that promote each disease in concert with environmental variables. This is consistent with the fact that classification of AD into these two subtypes appears arbitrary as there is no discernible early age-dependent peak in the population prevalence of dementia (of which AD contributes the majority of cases). The observed molecular pathway similarity between knock-in models of EOFAD and LOAD genetic variation supports the validity and utility of analysing early molecular events in AD pathogenesis using knock-in models in both zebrafish and mice, and that analysis of endogenous EOFAD mutations can contribute information that is valuable for understanding LOAD pathogenesis. The mystery of why GWAS has failed to detect variation in PSEN1, PSEN2 or APP in LOAD remains, although some mutations in PSEN2 and APP genes do cause later-onset familial forms of the disease (Cruchaga et al., 2012) and/or show incomplete penetrance (Finckh et al., 2000; Rossor et al., 1996; Sherrington et al., 1996; Thordardottir et al., 2018) or recessive inheritance (Di Fede et al., 2009; Tomiyama et al., 2008).

Brain hypometabolism is a diagnostic criterion for AD and can be visualised by 2-deoxy-2-(¹⁸F)fluoro-D-glucose positron emission tomography (¹⁸FDG-PET) (reviewed by Marcus et al., 2014). This technique has been used previously to investigate adults at risk of developing AD (e.g. Ou et al., 2019), whereas blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (MRI) has identified brain regional changes in activity (and hence energy consumption) in child carriers of the PSEN1 E280A (‘Paisa’) EOFAD mutation (Quiroz et al., 2015). This suggests that ¹⁸FDG-PET, or techniques such as dynamic glucose-enhanced MRI (Tolomeo et al., 2018), might enable the screening of individuals to determine the risk of later development of AD and be useful tools for investigating early energy changes in animal models of AD. Indeed, such techniques have been exploited to monitor energy consumption in the brains of transgenic mouse models (Huang et al., 2020; Luo et al., 2012; Poisnel et al., 2012; Tolomeo et al., 2018). However, they have not yet been applied to mice with knock-in EOFAD mutations. The observations of change in in vivo brain energy metabolism discussed above are consistent with the observed changes in the expression of genes of the oxidative phosphorylation pathway in the postmortem brains of early and late AD subjects relative to age-matched controls (Manczak et al., 2004). Neuronal cells derived from human induced pluripotent stem cells of LOAD patients also show increased expression of oxidative phosphorylation proteins and oxidative stress (Birnbaum et al., 2018), and neurons derived from a patient carrying the PSEN1^S170F EOFAD mutation showed mitochondrial abnormalities (Li et al., 2020).

Our findings regarding EOFAD mutations in zebrafish were not consistent with findings from our analysis of transcriptome data from homozygous App^NL-G-F mice. However, these transcriptome data were generated from ‘middle-aged’ (12 months old) mice rather than young adults, and the endogenous App gene of the mouse was altered with a total of six mutations (three that humanise the sequence of the Aβ region and three EOFAD mutations), motivated by the idea that the more aggregation-prone human Aβ sequence plays a critical role in the pathogenic mechanism of AD. Therefore, it is not directly comparable with our zebrafish EOFAD models, which contain single EOFAD-like mutations within single alleles of endogenous genes. To our knowledge, a transcriptome analysis has not been performed with App^NL-G-F mice at a younger age. However, the expression of genes involved in lysosomal function (KEGG_LYSOSOME) was observed to be highly significantly upregulated in homozygous App^NL-G-F brains. This is not unexpected, as acidification of the endolysosomal system is impaired by increased levels of the β-CTF fragment of APP [also known as C99 and generated by β-secretase cleavage of APP (Jiang et al., 2019)]. Increased β-CTF has been observed in the brains of App^NL-G-F mice (Saito et al., 2014). In a mouse model of a lysosomal storage disorder [glycogen storage disease type 2, which most seriously affects muscle (Yambire et al., 2019)], lysosomes failed to become sufficiently acidic, and this resulted in an intracellular ferrous iron deficiency and a pseudohypoxic response, mitochondrial dysfunction and inflammation (Yambire et al., 2019). [Degradation of HIF1-α, the master transcriptional regulator of the cellular response to hypoxia, is dependent on both oxygen and ferrous iron (Ivan et al., 2001).] Additionally, the App^NL-G-F mouse model shows increased levels of Aβ from a young age (Saito et al., 2014), and the deposition of Aβ into plaques was shown to be associated with the increased expression of genes in the complement system in a comprehensive spatial transcriptomics study of aging in the App^NL-G-F mouse model (Chen et al., 2020), providing another avenue for these mutations to trigger inflammation. Mitochondrial dysfunction has not been observed directly in App^NL-G-F mice. However, increased levels of oxidative stress have been observed at 12 months of age (Izumi et al., 2020), suggestive of increased reactive oxygen species that can be generated by dysfunction of mitochondrial respiration. Therefore, we suspect that processes similar to those in our zebrafish models and the APOE4 knock-in mice are being affected in App^NL-G-F mice. However, subtle signs of mitochondrial dysfunction in the transcriptome may be obscured by noise from the strong inflammatory signals in the bulk brain transcriptomic data (as well as confounding influences on the transcriptome analysis, such as litter-of-origin effects).

All of the mutations studied have also resulted in changes in the levels of transcripts required for ribosome formation. Protein translation is one of the most energy-costly processes within a cell (Buttgereit and Brand, 1995), and so expression of ribosomal proteins is modulated by the mTOR system, which surveys cellular nutrient status to adjust cellular metabolism (reviewed by Mayer and Grummt, 2006; Zhou et al., 2015). mTOR signalling, which appears to be increased in late-stage AD brains compared to controls (Griffin et al., 2005; Li et al., 2005; Sun et al., 2014), is regulated by growth factors, nutrients, energy levels and stress. In addition to ribosome biogenesis, mTOR signalling plays a role in various other cellular processes also implicated in AD pathogenesis, such as autophagy and lipid metabolism (reviewed by Saxton and Sabatini, 2017).

The mTOR proteins are localised at lysosomes within the mTORC1 and mTORC2 protein complexes (Sancak et al., 2010). Intriguingly, the v-ATPase complex that acidifies the endolysosomal pathway is required for mTORC1 activation (Zoncu et al., 2011). Proper assembly of the v-ATPase at the lysosome requires the PSEN1 protein (and this process is impaired in EOFAD patient fibroblasts) (Lee et al., 2010). Stimulation of mTOR signalling has been observed in response to accumulation of Aβ (Caccamo et al., 2010), and hyperactivation of mTOR is observed in Down's syndrome (in which the dosage of the APP gene is increased because it resides on chromosome 21 and early-onset AD is common) (Bordi et al., 2019; Iyer et al., 2014). Intriguingly, Bordi et al. observed that inhibition of mTOR signalling (specifically mTORC1) rescues autophagy and mitophagy defects in the fibroblasts of Down's syndrome individuals (Bordi et al., 2019). Among our transcriptome analyses of AD models, we only observed statistically significant changes to the expression of genes in the KEGG gene set for mTOR signalling in transheterozygous sorl1 mutants, in psen1^Q96_K97del/+ mutant zebrafish after acute hypoxia exposure, and in both male and female APOE4 mice. However, the majority of regulation of mTOR signalling occurs at the protein level, so that it is perhaps unsurprising that, in the EOFAD mutants, we could only detect significant changes in the transcriptional response to altered mTOR signalling rather than in the mTOR gene set itself.

In a highly sensitive analysis method such as RNA-seq, minimising external sources of variation increases resolving power. Our analysis has revealed that zebrafish can be highly advantageous for transcriptome profiling in the context of RNA-seq, as large numbers of progeny can be produced from a single pair mating, and these can subsequently be raised together in a single aquarium system, thus reducing genetic and environmental variation. This has allowed us to observe subtle effects due to the EOFAD-like mutations we have analysed. In contrast, a female mouse can only birth relatively small litters of 5-10 pups, making it particularly difficult to obtain sufficient numbers of synchronous sibling samples (particularly when genotypes of interest are produced by crossing of heterozygotes). Our re-analysis of the APOE-TR mouse brain transcriptomes was unable to distinguish with great certainty whether the effects we observed were due to Apoe genotype or litter-of-origin. Also, information on whether App^NL-G-F mice were littermates was not available, and this required us to assume that the effects of litter were negligible in order to perform the analysis. Another contrast between brain transcriptome analysis in zebrafish compared to mice is the influence of sex. Mouse brain transcriptomes show very significant differences due to sex, whereas sex has a negligible effect on bulk brain transcriptomes from zebrafish (Barthelson et al., 2020b, 2021a,b,c; Drew et al., 2012), and can generally be ignored in a differential expression analysis. We also found evidence for changes to cell-type proportions in both APOE-TR and App^NL-G-F mice, a phenomenon that can create the artefactual appearance of gene expression change. We have not observed cell-type proportion differences in 6-month-old or 24-month-old zebrafish (Barthelson et al., 2020b, 2021a,b,c; Hin et al., 2020, 2021; Fig. S14), possibly because of the resistance to damage of the highly regenerative zebrafish brain (Kroehne et al., 2011). Although this regenerative ability may hinder the use of zebrafish for studying overt neurodegeneration, it can facilitate analysis of young bulk brain transcriptomes before overt pathological processes would be expected.

The advantages of zebrafish for analysing the early effects of EOFAD mutations are countered, occasionally, by disadvantages. The teleost lineage in which zebrafish arose underwent an early whole-genome duplication event (reviewed by Meyer and Van de Peer, 2005), such that many human genes are represented by duplicate ‘co-orthologues’ in zebrafish (e.g. the co-orthologues of APP and APOE in zebrafish are appa/appb and apoea/apoeb, respectively). This complicates interpretation of the effects of mutations in these genes. Additionally, zebrafish have never been shown, definitively, to be capable of producing Aβ, a pathological hallmark of AD. The β-secretase (BACE) site of human APP does not appear to be conserved in zebrafish Appa and Appb (Moore et al., 2014). Whether Aβ accumulation is a cause or consequence of AD pathological processes continues as a matter of debate within the AD research community (reviewed by Morris et al., 2018). If zebrafish cannot produce Aβ, then the changes we have observed in the brains of our zebrafish models may illuminate Aβ-independent effects of EOFAD mutations.

Knock-in mouse models of single EOFAD mutations were generated 15-20 years ago (Guo et al., 1999; Kawasumi et al., 2004) but their brain transcriptomes have never been analysed in detail. This is probably because these mice showed only very mild cognitive phenotypes and lacked the AD histopathology currently used to define the disease [Aβ deposition and neurofibrillary tangles of tau protein (Jack et al., 2018)]. By expressing multiple mutant forms of EOFAD genes in transgenic mice, Aβ plaques can be detected and cognitive changes observed (reviewed by Esquerda-Canals et al., 2017; Myers and McGonigle, 2019). However, experience of using of many such ‘mouse models of AD’ has shown a lack of correlation of cognitive changes with Aβ levels (Foley et al., 2015) [Aβ levels do not closely correlate with cognitive changes in humans either (Giannakopoulos et al., 2003)], and transcriptome analysis of their brains has shown little to no concordance with transcriptomes from postmortem AD brains, or between the models themselves (Hargis and Blalock, 2017). In two papers, Saito et al. (2014, 2016) described phenotypic disparities between transgenic and APP EOFAD mutation knock-in mouse models. In the 2016 paper, they went so far as to declare that, “We recently estimated using single App knock-in mice that accumulate amyloid β peptide without transgene overexpression that 60% of the phenotypes observed in Alzheimer's model mice overexpressing mutant amyloid precursor protein (APP) or APP and presenilin are artifacts (Saito et al., 2014). The current study further supports this estimate by invalidating key results from papers that were published in Cell. These findings suggest that more than 3000 publications based on APP and APP/PS overexpression must be reevaluated.”

Nevertheless, since 2016, thousands more papers have been published using transgenic mouse models of AD. In this light, we were surprised to find that App^NL-G-F homozygous mice display a young adult brain transcriptome that is more severely disturbed than in the triple transgenic 3xTg-AD model – although that apparent disturbance is likely somewhat artefactual and due to changes in the relative proportions of different cell types in the model. As mentioned above, changes to cell-type proportions are not observed in our zebrafish models (Barthelson et al., 2020b, 2021b,c; Fig. S14).

Frustration with the difficulties of exploiting both transgenic and knock-in models of EOFAD mutations in mice has contributed to the drive for examining knock-in mouse models of LOAD risk variants, such as now conducted by the Model Organism Development and Evaluation for Late-Onset Alzheimer's Disease Consortium (Oblak et al., 2020). The brain transcriptome similarities seen between our single mutation heterozygous EOFAD mutation-like knock-in zebrafish models and the knock-in APOE ε4 mice strongly support the informative value of these models, and imply that heterozygous EOFAD mutation knock-in mouse models offer a path forward, particularly in understanding the earliest molecular events that lead to AD.

For analysis of 6-month-old and 24-month-old genome-edited Tübingen zebrafish (sample sizes and sexes are indicated in Table 1 and Figs S1-S6), we obtained the differential gene expression analysis outputs and harmonic mean P-values (statistical significance of gene sets) from each individual analysis (see Barthelson et al., 2020b, 2021a,b,c). For these analyses, differential gene expression analysis was performed using edgeR (Robinson et al., 2009), and enrichment analysis was performed by calculation of the harmonic mean P-value (Wilson, 2019) of the raw P-values of three methods of ranked-list-based enrichment analyses: fry (Wu et al., 2010), camera (Wu and Smyth, 2012) and GSEA (Subramanian et al., 2005) [as implemented in the fgsea R package (Sergushichev, 2016 preprint)]. We used the harmonic mean P-value to determine the overall significance of changes to gene expression within gene sets because this method does not assume that component P-values are independent (Wilson, 2019). We have previously validated the use of the harmonic mean P-value on simulated RNA-seq datasets (Barthelson et al., 2020b). We considered a gene set to be significantly altered if the false discovery rate (FDR)-adjusted harmonic mean P-value remained below 0.05. The gene sets used for enrichment analysis were the KEGG (Kanehisa and Goto, 2000) gene sets, so that any changes to gene expression in any of 186 biological pathways/processes could be determined. Additionally, we used our recently defined IRE gene sets (Hin et al., 2021) to test for evidence of iron dyshomeostasis. For the K97fs and Q96_K97del analyses, enrichment analysis was not performed on the KEGG gene sets in the original analyses. Therefore, we performed the enrichment analysis as described above for these datasets. For the K97fs analysis, we obtained the gene-level counts and the results of the differential gene expression analysis described by Hin et al. (2020a) from www.github.com/UofABioinformaticsHub/k97fsZebrafishAnalysis. For the Q96_K97del analysis, we obtained the gene-level counts and the results of the differential gene expression analysis described by Hin et al. (2020b) from the first author of the cited paper. Note that, for each zebrafish analysis, the sample size was usually n=6 zebrafish per genotype, based on our previous calculation that this sample size should give ∼70% power to detect the majority of expressed transcripts in a zebrafish brain transcriptome at a fold change of >2 and at an FDR of 0.05 (Barthelson et al., 2020b).

The raw .CEL files were obtained from GEO and analysed with R (https://www.r-project.org/). Pre-processing was performed using the RMA (Irizarry et al., 2003) method as implemented in the oligo package (Carvalho and Irizarry, 2010). We omitted any probesets that contained a median log2 intensity value of <3.5 (lowly expressed) and also any probesets assigned to multiple genes. Differential gene expression analysis was performed using limma (Ritchie et al., 2015), specifying pairwise contrasts between the App^NL-G-F homozygous mice or the 3xTg homozygous mice with their respective controls by using a contrasts matrix. We considered a probeset to be differentially expressed in each contrast if the FDR-adjusted P-value was <0.05. For over-representation of the KEGG and IRE gene sets within the differentially expressed genes, we used kegga (Young et al., 2010). We also performed ranked-list-based enrichment analysis using the harmonic mean P-value as described for the zebrafish analyses.

We obtained the raw fastq files for the entire APOE-TR RNA-seq experiment from the AD Knowledge Portal (accession number syn20808171, https://adknowledgeportal.synapse.org/). The raw reads were first processed using AdapterRemoval (version 2.2.1) (Schubert et al., 2016), with following options selected: --trimns, --trimqualities, --minquality 30 and --minlength 35. Then, the trimmed reads were aligned to the Mus musculus genome (Ensembl GRCm38, release 98) using STAR (version 2.7.0) (Dobin et al., 2013), using default parameters to generate .bam files. These bam files were then sorted and indexed using SAMtools (version 1.10) (Li et al., 2009). The gene expression counts matrix was generated from the bam files using featureCounts (version 1.5.2) (Liao et al., 2014). We only counted the number of reads that uniquely aligned to, strictly, exons with a mapping quality of at least ten to predict expression levels of genes in each sample.

We then imported the output from featureCounts (Liao et al., 2014) for analysis with R. We first omitted genes that were lowly expressed (and are uninformative for differential expression analysis). We considered a gene to be lowly expressed if it contained, at most, 2 counts per million in 8 or more of the 24 samples we analysed. We also assessed whether the sex of each sample was correctly classified by examining the expression of genes that are located on the Y chromosome. Three samples appeared to be classified incorrectly and were subsequently corrected (Fig. S16).

To determine which genes were dysregulated in APOE4 and APOE2 mice relative to APOE3, we performed a differential gene expression analysis using a generalised linear model and likelihood ratio tests using edgeR (McCarthy et al., 2012; Robinson et al., 2009). We chose a design matrix that specified the APOE genotype and sex of each sample. The contrasts matrix was specified to compare the effect of APOE2 or APOE4 relative to APOE3 in males and in females. In this analysis, we considered a gene to be differentially expressed if the FDR-adjusted P-value was <0.05. A bias for longer transcript length and higher GC content percentage was observed in this dataset. Therefore, we corrected for this bias using CQN (Hansen et al., 2012). We calculated the average transcript length per gene, and a weighted (by transcript length) average GC content percentage per gene, as input to CQN to produce the offset to correct for the bias. This offset was then included in an additional generalised linear model and likelihood ratio tests in edgeR with the same design and contrast matrices. For over-representation of the KEGG and IRE gene sets within the differentially expressed genes, we used goseq (Young et al., 2010), specifying average transcript length to generate the probability weighting function, which corrects for the probability that a gene is classified as differentially expressed based on its transcript length (average transcript length per gene) alone. We also performed ranked-list-based enrichment analysis as described for the zebrafish analysis.

Visualisation of gene expression data throughout this analysis was performed using ggplot2 (Wickham, 2016), pheatmap (available at https://CRAN.R-project.org/package=pheatmap) and pathview (Luo et al., 2017). The code used to perform the analysis in this study can be found at www.github.com/karissa-b/AD-signature.

We thank Dr Nhi Hin for providing the Q96_K97del and K97fs zebrafish datasets, and original analyses. We also thank Dr Stephen Pederson for his advice and guidance throughout the project. This work was supported by supercomputing resources provided by the Phoenix High Performance Computing service at the University of Adelaide. The results for the young APOE-TR mice were based on data obtained from the AD Knowledge Portal (https://adknowledgeportal.synapse.org/). We thank Drs Patrick Sullivan and Nobuyo Maeda for generating human APOE-TR mice and providing access through Taconic.