Abstract
Background:
5XFAD humanized mutant mice and Trem2 knockout (T2KO) mice are two mouse models relevant to the study of Alzheimer’s disease (AD)-related pathology.
Objective:
To determine hippocampal transcriptomic and polyadenylation site usage alterations caused by genetic mutations engineered in 5XFAD and T2KO mice.
Methods:
Employing a publicly available single-nucleus RNA sequencing dataset, we used Seurat and Sierra analytic programs to identify differentially expressed genes (DEGs) and differential transcript usage (DTU), respectively, in hippocampal cell types from each of the two mouse models. We analyzed cell type-specific DEGs further using Ingenuity Pathway Analysis (IPA).
Results:
We identified several DEGs in both neuronal and glial cell subtypes in comparisons of wild type (WT) versus 5XFAD and WT versus T2KO mice, including Ttr, Fth1, Pcsk1n, Malat1, Rpl37, Rtn1, Sepw1, Uba52, Mbp, Arl6ip5, Gm26917, Vwa1, and Pgrmc1. We also observed DTU in common between the two comparisons in neuronal and glial subtypes, specifically in the genes Prnp, Rbm4b, Pnisr, Opcml, Cpne7, Adgrb1, Gabarapl2, Ubb, Ndfip1, Car11, and Stmn4. IPA identified three statistically significant canonical pathways that appeared in multiple cell types and that overlapped between 5XFAD and T2KO comparisons to WT, including ‘FXR/RXR Activation’, ‘LXR/RXR Activation’, and ‘Acute Phase Response Signaling’.
Conclusion:
DEG, DTU, and IPA findings, derived from two different mouse models of AD, highlight the importance of energy imbalance and inflammatory processes in specific hippocampal cell types, including subtypes of neurons and glial cells, in the development of AD-related pathology. Additional studies are needed to further characterize these findings.
Keywords
INTRODUCTION
Alzheimer’s disease (AD) is the most common form of dementia [1], and is one of society’s largest healthcare burdens. Given the limitations associated with studies in humans, mouse models are valuable tools for studying AD. One such model is the humanized 5XFAD mouse, which harbors five mutations known to predispose people to familial AD [2]. 5XFAD mice exhibit deficits in spatial memory, spatial learning, fear conditioning, and impaired social conditioning/social behavior [2–8]. Additionally, they exhibit early amyloid pathology, gliosis, neuronal loss, and myelin deficits [2–4, 8].
Another mouse model relevant to AD is the Trem2 knockout (T2KO). The human TREM2 gene encodes a protein known as triggering receptor expressed on myeloid cells 2. It is found in a variety of immune cells in the periphery and central nervous system. In the brain, TREM2 is expressed primarily in microglia, a subtype of brain cell important for normal processing of amyloid deposits [9]. Certain loss of function variants of TREM2 increase risk for AD, suggesting a role in pathophysiology (for review, see [10]). Furthermore, loss of Trem2 function in mouse models of AD results in worsening amyloid pathology [11–14]. Together, T2KO and 5XFAD mice are recognized as valuable models to study the etiology of AD, and in fact, the two models have been integrated recently to investigate the interaction among the relevant genetic factors [14].
Single-cell (or single-nucleus) RNA sequencing (scRNAseq or snRNAseq) allows transcriptomic analysis of single cells and provides insight into the function of specific cell types. It also allows comparison between normal and disease states and identification of differentially expressed genes (DEGs). Relevant to AD, a recent snRNAseq study analyzed transcriptomes from the prefrontal cortices of 48 individuals with AD, and identified DEGs in non-neuronal cells involved in myelination and inflammation, suggesting a role for neuronal-glial interactions in AD pathology and highlighting the importance of distinguishing gene expression characteristics in specific cell subtypes [15, 16]. SnRNAseq studies conducted in AD mouse models can support and extend human studies aimed at identifying the cells and genes with altered transcription that are associated with development of AD-like phenotypes.
Common methods of analyzing single cell/nucleus data for DEGs do not distinguish between alternative mRNA isoforms. However, new bioinformatic methods demonstrate that data from scRNAseq or snRNAseq experiments can be mined to detect differential transcript usage (DTU) [17]. DTU specifically refers to differential polyadenylation site usage. We recently detected DTU in hippocampal cells from Apoe knockout mice [18], another mouse model of relevance to AD [19]. DTU analysis of scRNAseq and snRNAseq data from other AD models could provide further insight into the molecular mechanisms underlying disease pathophysiology.
We based the present study on the hypothesis that the genetic mutations engineered in 5XFAD and T2KO mice cause transcriptomic and polyadenylation site usage alterations of genes implicated in AD pathology, and that these alterations would have both unique and overlapping features between the models. To test this hypothesis, we analyzed publicly available snRNAseq data from pooled hippocampi of 15-month-old Trem2 knockout (T2KO) mice, 5XFAD (named WT-5XFAD in the original dataset [14]) mice, and wild type (WT) mice. While many areas of the brain are relevant to understanding the complex phenotype that is AD, we chose to look at mouse hippocampus data as the hippocampus is one of the main brain areas affected in AD, is involved in memory and cognitive function, and is impacted in the early stages of AD [20]. For these reasons, it is a logical region to investigate when looking for transcriptomic effects of mutations that cause AD in humans. As hypothesized, results reveal partially overlapping transcriptomic profiles in the two models. Specifically, we confirm DEGs in microglia and oligodendrocytes reported in prior studies, identify functional pathways in astrocytes related to inflammatory processes, and characterize novel and potentially important transcriptomic alterations in a subtype of excitatory neurons. Additionally, we obtained DTU findings which implicate novel genes of potential relevance to AD pathogenesis.
MATERIALS AND METHODS
Datasets
All analyses were performed on GEO dataset GSE140399 [dataset] [14] comprised of snRNAseq data generated by droplet-based Chromium Single Cell 3’ Reagent Kit (10x Genomics) on pooled RNA samples. We analyzed three pooled RNA samples (one pool of RNA for each genotype), and each pool contained RNA from three mice. These RNA pools came from hippocampi samples from 15-month-old C57BL/6J WT male mice (Sample GSM4160647), 15-month-old T2KO mice created on a C57BL/6J background [12, 13] (Sample GSM4160648), and 15-month-old 5XFAD mice produced on a C57BL/6J background [3, 4] (Sample GSM4160649).
Seurat analysis
We loaded the filtered matrix of gene-by-cell expression, barcodes file, and genes file for the three samples into Seurat (version 4.0.3 [21]) and used this to generate three Seurat objects. The merge method was then used to combine these three Seurat objects into a single Seurat object for analysis. Quality control (QC) metrics included filtering out nuclei with over 2500 or less than 200 genes, and/or containing greater than 5% of total transcripts from mitochondrial genes. We normalized the data (using the default “LogNormalize”) and identified the 2000 most highly variable genes via the FindVariableFeatures command with the “vst” selection method. After scaling by linear transformation, we performed linear dimension reduction using principal component analysis (PCA) on the scaled data. We clustered cells using the first 50 PCs at resolution = 0.5, yielding 18 distinct clusters. Two analyses were done using the FindMarkers function. First, clusters were compared to identify gene markers for each cluster, using default statistical settings of Wilcoxon Rank Sum test and Bonferroni correction, but with logfc.threshold = 0. Secondly, the FindMarkers function was used to identify DEGs between the three genotypes which were used in our analysis, again using default statistical settings of Wilcoxon Rank Sum test and Bonferroni correction, but with logfc.threshold = 0 and not including mitochondrial genes. Log2(fold-change) values of the average expression between the two groups in each comparison were obtained (positive values indicate that the gene is more highly expressed in the first group within a particular comparison). Adjusted p-values, based on Bonferroni correction using all genes in the dataset, were obtained for reported DEGs. We also used these data in combination with known markers for hippocampal cell types [22–24] to determine the cell type of each cluster.
We used Pearson’s chi-squared contingency table analysis to compare the number of cells in each cluster between WT and 5XFAD mice and between WT and T2KO mice (Fig. 1A; Fig. 1B; Supplementary File 1). Each chi-squared analysis was run comparing the number of cells analyzed in the indicated cluster to the rest of the cells in that snRNAseq sample.

Cell counts per cluster (cell type) in (A) WT versus 5XFAD comparison and (B) WT versus T2KO comparison. Pearson’s chi-squared contingency table analysis was run, comparing the number of cells analyzed for each comparison. After multiple testing correction, there were statistically-significant differences for some of the cluster comparisons. These are indicated by asterisks. Cluster 0 = C1ql2+ Excitatory neurons (Excitatory neuron 1); Cluster 1 = Fibcd1+ Excitatory neurons (Excitatory neuron 2); Cluster 2 = Plxdc1+ Excitatory neurons (Excitatory neuron 3); Cluster 3 = Cpne4+ Excitatory neurons (Excitatory neuron 4); Cluster 4 = Slc1a3+ Resting astrocytes; Cluster 5 = Pde4d+ Excitatory neurons (Excitatory neuron 5); Cluster 6 = Sst+ Inhibitory neurons (Inhibitory neuron 1); Cluster 7 = Tmem91+ Inhibitory neurons (Inhibitory neuron 2); Cluster 8 = Erbb4+ Inhibitory neurons (Inhibitory neuron 3); Cluster 9 = Map2+ Excitatory neurons (Excitatory neuron 6); Cluster 10 = Ndst4+ Excitatory neurons (Excitatory neuron 7); Cluster 11 = Trf+ Myelin-forming mature oligodendrocytes (MFOLs); Cluster 12 = Pdgfra+ Oligodendrocyte precursor cells (OPCs); Cluster 13 = Slc1a3+, Gfap+ Reactive astrocytes; Cluster 14 = Cldn5+, Dcn+ Endothelial cells/vascular leptomeningeal cells (VLMCs) (Endothelial cells/VLMCs); Cluster 15 = Kcnj13+ Choroid plexus epithelial cells; Cluster 16 = C1qa+ Microglia; Cluster 17 = Calb2+ Excitatory neurons (Excitatory neuron 8).
Sierra analysis
We used Regtools [25] from the Github repository to extract the transcript junction information in BED (browser extensible data) format from the WT, T2KO, and 5XFAD BAM alignment files. We generated the reference file from the GRCm38/mm10 genome build. The whitelist.bc files were from the WT, T2KO, and 5XFAD barcodes.tsv files obtained from GSE140399. We called peaks using Sierra (version 0.99.26 [17]) separately for each of the WT, T2KO, and 5XFAD bamfiles. We then merged the three files of called peaks in Sierra. We counted peaks for WT, T2KO, and 5XFAD samples separately and aggregated them into one file. We annotated peaks with BSgenome.Mmusculus.UCSC.mm10 genome file and transferred to a Seurat object. We examined DTU in the 18 clusters using the DUTest command with the parameter “exp.thres = 0.1”. We filtered DUTest results using an FDR adjusted p-value < 0.05 and log2(fold-change) threshold of 0.25 (approximately 20% difference). We then used the DUTest command to generate a list of gene peaks differentiating the WT, T2KO, and 5XFAD samples with thresholds as noted above. We generated lists of relative peak expression by cluster for the WT versus T2KO samples and the WT versus 5XFAD samples for those peaks with a fold change ≥2.0 in at least one cluster.
We compared FindMarkers gene hits— DEGs in each cell cluster with an FDR-adjusted p-value < 0.05 to genes containing DTU peaks with a fold change > 2.0 in the corresponding cluster in order to find lists of genes which overlap in thisanalysis.
Ingenuity pathway analysis: Canonical pathways
After removing mitochondrial DNA gene transcripts, and transcripts without a known Ensembl ID, we used DEGs in each cluster with an FDR-adjusted p-value < 0.05 in IPA (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis) [26]. We performed canonical pathway analysis and generated a list of terms with a q-value < 0.05 after FDR correction for multiple testing.
RESULTS
Seurat analysis: Transcriptomics
Analysis of the combined WT, T2KO, and 5XFAD data identified 18 separate clusters corresponding to different cell type identities (Fig. 2; separate uniform manifold and projection (UMAP) plots of WT, 5XFAD, and T2KO samples are shown in Supplementary Figure 1). The choroid plexus appears to have been captured alongside the hippocampus. Chi-squared analysis showed statistically significant differences for many of the cluster comparisons between the WT and 5XFAD mice in numbers of neurons (excitatory and inhibitory), astrocytes (resting and reactive), oligodendrocyte precursor cells (OPCs), and microglia (Fig. 1A; Supplementary File 1). When comparing the WT and T2KO mice, we found differences in numbers of subtypes of excitatory neurons (Fig. 1B; Supplementary File 1).

Uniform manifold and projection (UMAP) plot of wild type (WT), Trem2 knockout (T2KO), and 5XFAD (5XFAD) hippocampus samples combined. A) The transcriptome profiles of 10230 (3325 WT, 2906 T2KO, and 3999 5XFAD) single nuclei were utilized for unbiased clustering in Seurat and are presented as a UMAP dimension reduction plot of all single nuclei. Nuclei are color-coded by cluster, which corresponds to cell type. B) A dot plot was generated in Seurat and clusters were annotated using genes known to be markers of different cell types. The size and color of the dots is proportional to the percentage of nuclei expressing the gene of interest (Pct. Exp.) and the average expression level of the gene (Avg. Exp.), respectively. Cluster numbers and colors match those of the UMAP in Fig. 2A. Cluster 0 = C1ql2+ Excitatory neurons (Excitatory neuron 1); Cluster 1 = Fibcd1+ Excitatory neurons (Excitatory neuron 2); Cluster 2 = Plxdc1+ Excitatory neurons (Excitatory neuron 3); Cluster 3 = Cpne4+ Excitatory neurons (Excitatory neuron 4); Cluster 4 = Slc1a3+ Resting astrocytes; Cluster 5 = Pde4d+ Excitatory neurons (Excitatory neuron 5); Cluster 6 = Sst+ Inhibitory neurons (Inhibitory neuron 1); Cluster 7 = Tmem91+ Inhibitory neurons (Inhibitory neuron 2); Cluster 8 = Erbb4+ Inhibitory neurons (Inhibitory neuron 3); Cluster 9 = Map2+ Excitatory neurons (Excitatory neuron 6); Cluster 10 = Ndst4+ Excitatory neurons (Excitatory neuron 7); Cluster 11 = Trf+ Myelin-forming mature oligodendrocytes (MFOLs); Cluster 12 = Pdgfra+ Oligodendrocyte precursor cells (OPCs); Cluster 13 = Slc1a3+, Gfap+ Reactive astrocytes; Cluster 14 = Cldn5+, Dcn+ Endothelial cells/vascular leptomeningeal cells (VLMCs) (Endothelial cells/VLMCs); Cluster 15 = Kcnj13+ Choroid plexus epithelial cells; Cluster 16 = C1qa+ Microglia; Cluster 17 = Calb2+ Excitatory neurons (Excitatory neuron 8).
We identified multiple DEGs between 5XFAD and WT mice and between T2KO and WT mice (Fig. 3; Supplementary File 2). The two mouse models share 13 DEGs in specific neuronal and glial cell types including Ttr, Fth1, Pcsk1n, Malat1, Rpl37, Rtn1, Sepw1, Uba52, Mbp, Arl6ip5, Gm26917, Vwa1, and Pgrmc1 (Table 1). Among the hippocampal cell types analyzed, we found the most DEGs in C1ql2+ excitatory neurons. Correspondingly, we found the most overlapping DEGs (i.e., found in both mouse models) in this cell type as well. We note that all but two of these DEGs are found in this cluster. DEGs found in the most clusters are Ttr (6 of 18 clusters), Fth1 (5 of 18 clusters), and Rtn1 (3 out of 18 clusters) (Table 1; Supplementary File 2).

Number of differentially expressed genes (DEGs) per cluster (cell type) in both comparisons: WT versus 5XFAD and WT versus T2KO. Analysis using the FindMarkers function in Seurat, which generated DEGs between WT versus 5XFAD mice and WT versus T2KO mice, respectively, for each of 18 separate clusters corresponding to different cell identities. The numbers listed here correspond to the number of DEGs with a Bonferroni corrected p-value < 0.05. Cluster 0 = C1ql2+ Excitatory neurons (Excitatory neuron 1); Cluster 1 = Fibcd1+ Excitatory neurons (Excitatory neuron 2); Cluster 2 = Plxdc1+ Excitatory neurons (Excitatory neuron 3); Cluster 3 = Cpne4+ Excitatory neurons (Excitatory neuron 4); Cluster 4 = Slc1a3+ Resting astrocytes; Cluster 5 = Pde4d+ Excitatory neurons (Excitatory neuron 5); Cluster 6 = Sst+ Inhibitory neurons (Inhibitory neuron 1); Cluster 7 = Tmem91+ Inhibitory neurons (Inhibitory neuron 2); Cluster 8 = Erbb4+ Inhibitory neurons (Inhibitory neuron 3); Cluster 9 = Map2+ Excitatory neurons (Excitatory neuron 6); Cluster 10 = Ndst4+ Excitatory neurons (Excitatory neuron 7); Cluster 11 = Trf+ Myelin-forming mature oligodendrocytes (MFOLs); Cluster 12 = Pdgfra+ Oligodendrocyte precursor cells (OPCs); Cluster 13 = Slc1a3+, Gfap+ Reactive astrocytes; Cluster 14 = Cldn5+, Dcn+ Endothelial cells/vascular leptomeningeal cells (VLMCs) (Endothelial cells/VLMCs); Cluster 15 = Kcnj13+ Choroid plexus epithelial cells; Cluster 16 = C1qa+ Microglia; Cluster 17 = Calb2+ Excitatory neurons (Excitatory neuron 8).
Differentially expressed genes (DEGs) detected in both comparisons. WT versus 5XFAD\\ and WT versus T2KO
Apart from the overlap between multiple DEGs in several subtypes of excitatory neurons, we found DEGs in common between the mouse models in only two other cell types: Sst+ inhibitory neurons, in which Ttr and Fth1 were differentially expressed, and resting astrocytes, in which Ttr was differentially expressed. Thus, although we observe major transcriptomic differences in several excitatory neuronal hippocampal cell types in mouse models of AD, we also observe major differences in inhibitory neurons and non-neuronal cell types.
Consistent with prior work [14], we find that, compared to WT, complement component C4b, the serine (or cysteine) peptidase inhibitor, clade A, member 3 N (Serpina3n) and histocompatibility 2, D region locus 1 (H2-d1) are upregulated in myelin-forming mature oligodendrocytes (MFOLs) from 5XFAD mice. In addition, we find significant upregulation of Cd9, Klk6, and B2m in MFOLs from 5XFAD mice. Although Thy1 also appears upregulated, this is likely an artifact due to the transgene structure used in the creation of 5XFAD mice [2–4]. We did not detect changes in any of these genes when comparing T2KO to WT, suggesting they are not relevant to pathologies associated specifically with loss of Trem2 function.
Changes in expression of C4b occur in other cell types apart from oligodendrocytes. Prior work indicated that astrocytes from 5XFAD mice upregulated Gfap and C4b, but to a lesser extent than did oligodendrocytes [14]. We also find that Gfap and C4b are upregulated in resting astrocytes when comparing 5XFAD mice with WT. In our analysis, the increase in C4b expression is greater in MFOLs than in resting astrocytes, but the increase in Gfap expression is greater in resting astrocytes. Differences in C4b and Gfap expression are not observed when we compare T2KO with WT mice. Prior studies also report strong upregulation of C4b in OPCs, with a fold-change and p-value similar to genes activated in oligodendrocytes [14]. We see C4b upregulation in both OPCs and MFOLs, but again, only in the comparison between 5XFAD mice and WT. Elevations in C4B in cortex are associated with increased synaptic pruning and schizophrenia [27, 28]. We do not observe differences in expression of Gfap or C4b in reactive astrocytes in either 5XFAD or T2KO mice, suggesting that these genes do not play a major role in this cell type with respect to AD pathogenesis.
Sierra analysis: Differential transcript usage
Comparison of differential transcript usage (DTU) between WT and 5XFAD and between WT and T2KO generated lists of transcripts for each gene in each cluster for which there was a fold change of ∼20% or more (log2FC≥0.25) (Table 2; Fig. 4; Supplementary File 3). We identified DTU that overlapped between the two comparisons in neuronal and glial cell types. Specifically, these occurred in the genes Prnp, Rbm4b, Pnisr, Opcml, Cpne7, Adgrb1, Gabarapl2, Ubb, Ndfip1, Car11, and Stmn4 (Table 2).
Gene peaks with Differential Transcript Usage (DTU) overlapping in both comparisons. WT versus\\ 5XFAD and WT versus T2KO

Differential transcript usage per cluster (cell type) in both comparisons: WT versus 5XFAD and WT versus T2KO. Analysis using the DUTest function in Sierra generated lists of peaks for each gene in each of 18 separate clusters corresponding to different cell types. The numbers listed here are the number of gene peaks that had an absolute fold change (AbsFC)≥2.0 for each comparison. Cluster 0 = C1ql2+ Excitatory neurons (Excitatory neuron 1); Cluster 1 = Fibcd1+ Excitatory neurons (Excitatory neuron 2); Cluster 2 = Plxdc1+ Excitatory neurons (Excitatory neuron 3); Cluster 3 = Cpne4+ Excitatory neurons (Excitatory neuron 4); Cluster 4 = Slc1a3+ Resting astrocytes; Cluster 5 = Pde4d+ Excitatory neurons (Excitatory neuron 5); Cluster 6 = Sst+ Inhibitory neurons (Inhibitory neuron 1); Cluster 7 = Tmem91+ Inhibitory neurons (Inhibitory neuron 2); Cluster 8 = Erbb4+ Inhibitory neurons (Inhibitory neuron 3); Cluster 9 = Map2+ Excitatory neurons (Excitatory neuron 6); Cluster 10 = Ndst4+ Excitatory neurons (Excitatory neuron 7); Cluster 11 = Trf+ Myelin-forming mature oligodendrocytes (MFOLs); Cluster 12 = Pdgfra+ Oligodendrocyte precursor cells (OPCs); Cluster 13 = Slc1a3+, Gfap+ Reactive astrocytes; Cluster 14 = Cldn5+, Dcn+ Endothelial cells/vascular leptomeningeal cells (VLMCs) (Endothelial cells/VLMCs); Cluster 15 = Kcnj13+ Choroid plexus epithelial cells; Cluster 16 = C1qa+ Microglia; Cluster 17 = Calb2+ Excitatory neurons (Excitatory neuron 8).
Our analysis revealed a number of clusters having genes with DTU in both mouse models. Unlike our findings with DEGs, T2KO mice had more genes with DTU versus WT than did 5XFAD mice; however, there are a number of genes with cluster-specific DTU and several genes with cluster-overlapping DTU in the two mouse models. Figure 4 and Supplementary File 3 list the number of genes with DTU in each cluster for each mouse model. In general, we detected the greatest number of genes with DTU in the largest clusters (i.e., subtypes of excitatory neurons), and this was independent of mouse model.
When comparing FindMarkers gene hits to genes containing DTU peaks with a fold change > 2.0, there are two genes that overlap in this analysis in the WT versus 5XFAD comparison. Thy1 overlaps in Fibcd1+ excitatory neurons, Cpne4+ excitatory neurons, and Ndst4+ excitatory neurons. Again, this is likely an artifact due to the transgene structure used in the creation of 5XFAD mice [2–4]. Ntm overlaps in resting astrocytes. No overlaps between DEGs and DTU genes were found in the WT versus T2KO comparison (Supplementary File 4).
Ingenuity pathway analyses: Canonical pathways
Input of the cluster-specific differentially expressed gene (DEG) lists into Ingenuity Pathway Analysis (IPA) software resulted in the identification of sixteen canonical pathway terms that appeared in more than one cluster in the WT versus 5XFAD comparison, and six canonical pathway terms that appeared in more than one cluster in the WT versus T2KO comparison (Supplementary File 5). Thirty-three canonical pathway terms appeared in one cluster in the WT versus 5XFAD comparison, and one appeared in one cluster in the WT versus T2KO comparison (Supplementary File 5). The canonical pathways ‘FXR/RXR Activation’, ‘LXR/RXR Activation’, and ‘Acute Phase Response Signaling’ appeared in more than one cell type and overlapped between 5XFAD and T2KO.
Input of the cluster-specific gene list which contain DTU into IPA software resulted in the identification of twenty-five canonical pathway terms that appeared in more than one cluster in the WT versus 5XFAD comparison, and fourteen canonical pathway terms that appeared in more than one cluster in the WT versus T2KO comparison (Supplementary File 6). Fifty-seven canonical pathway terms appeared in one cluster in the WT versus 5XFAD comparison, and twenty-six appeared in one cluster in the WT versus T2KO comparison (Supplementary File 6). The canonical pathways ‘Synaptogenesis Signaling Pathway’ and ‘SNARE Signaling Pathway’ appeared in more than one cell type and overlapped between 5XFAD and T2KO.
DISCUSSION
This secondary data analysis provides deeper characterization of transcriptome alterations associated with AD-related pathophysiology in 5XFAD and T2KO mouse models by performing DTU analysis on existing datasets. We performed DTU analysis as a downstream extension of the pipeline used to generate data in Seurat, and used workflow parameters similar to those used in our previous analysis [18]. Notably, these parameters differ from those used for analysis of the dataset originally [14], and although we confirmed some results, there are certain differences between DEGs reported in that study and those we report in this reanalysis.
An important overlap in results between the original analysis [14] and the present work is in the expansion of the microglial cell pool in 5XFAD mice (Fig. 1A; Fig. 2; Supplementary File 1), which likely reflects Aβ plaque-induced or neurofibrillary tangle-induced activation of the immune defense system. Since Trem2 expression partially drives this response, the size of the microglial pool in T2KO mice is not different from WT. In addition, we detected 3 DEGs in microglia when we compared 5XFAD mice to WT: Prox1, a circadian transcription factor, Mgat1, an acetylglucosaminyltransferase, and Med30, a peroxisome proliferator activator-related transcription factor. Expression of these genes is not altered in T2KO mice, suggesting they are not relevant to Trem2 deficiency per se. Disease-associated microglia (DAM) DEGs described in the original study [14] are not detected in our reanalysis, likely due to differences in bioinformatic workflow.
Overall, we detected more cell type-specific DEGs when analyzing data from 5XFAD mice compared to T2KO mice, suggesting that expression of human AD mutations has a greater impact on the transcriptome of brain cells than loss of Trem2 function. Interestingly, a small subset of genes is impacted in both models, and we hypothesize that DEGs found in the same cell type in both 5XFAD and T2KO mice are particularly relevant to AD-like pathology. Similarly, we suggest that the overlapping DEGs in these common cell types are of primary importance with regards to pathophysiological mechanisms. Pearson’s chi-squared contingency table analysis shows statistically significant differences for many of the cluster comparisons between the WT and 5XFAD mice, and between the WT and T2KO mice (Fig. 1A; Fig. 1B; Supplementary File 1). Confirming this by replication could be the focus of future studies.
The cell types with shared DEGs in both mouse models are mainly subtypes of excitatory neurons, a finding consistent with recent work indicating that alterations in gene expression in excitatory neurons are relevant to phenotypes associated with 5XFAD mice [29]. Interestingly, the genes represented by the excitatory neuron DEGs detected in our analysis are involved in various proposed mechanisms of AD pathophysiology. In particular, Rtn1 is involved in amyloid precursor protein processing [30], probably by blocking BACE1 activity [31], and its protein expression level is reduced in AD brain as compared to controls [32]. Our data show that Rtn1 is downregulated in 5XFAD mice, which is consistent with the human AD data. However, in T2KO mice, Rtn1 is upregulated compared to controls. Ttr is involved in transport of endogenous substances including thyroid hormone and retinol [33]. TTR protein levels are decreased in CSF [33–35] and plasma [33, 37] of AD patients. It is strongly linked to AD through both genetic [38] and epigenetic studies [39]. It is also suggested to have a role in neuroprotection in AD [40]. Stable TTR protein has been shown to be important in clearing Aβ, and stability of TTR is reduced in AD [41]. This is consistent with our data showing that Ttr is downregulated in both 5XFAD and T2KO mice as compared to controls. Iododiflunisal, an agent which stabilizes TTR, is being investigated in the development of AD therapeutics [42]. Fth1 is involved in iron utilization and in control of cellular redox state [43, 44]. It is a paralog of Ftmt, a gene implicated in inflammatory brain pathology of relevance to AD [45]. An increase in Ftmt mRNA and associated protein levels are shown to be increased in the brains of AD patients [45, 46]. This is consistent with our data in which Fth1 levels are increased in both 5XFAD and T2KO mice. Thus, our data suggest that alteration in expression level of these three genes in specific subtypes of excitatory neurons in the hippocampus of two different AD-related mouse lines may provide clues to final common pathways of pathology caused by dysregulation of amyloid processing.
The original analysis of these mouse model data [14] also reported several of these genes as DEGs, but with some differences noted. Findings were not reported with respect to Rtn1, Ttr, or Fth1 DEGs in 15-month-old 5XFAD or T2KO mice; however, that study did find that Fth1 was a DEG in 7-month-old 5XFAD mice in neuronal and oligodendrocyte clusters [14]. Similarly, Rtn1 is a reported DEG in oligodendrocytes and astrocytes from 5XFAD mice, and also in neurons and microglia from both 5XFAD and T2KO mice [14]. Interestingly, FTH1 is reported as being downregulated in microglia in human AD samples [14]. Overall, the present reanalysis of the previously published data [14] is consistent with data from human studies of AD and provides additional support for the use of mouse models to gain further insight into these genes as factors in AD pathophysiology.
Several of the canonical pathways derived from our differential gene expression analysis have been implicated in AD. With respect to 5XFAD mice, one such pathway, the ‘EIF2 Signaling’ pathway refers to genes that interact with or regulate eukaryotic initiation factor 2 (EIF2), an initiation factor involved in the regulation of protein synthesis. This ‘EIF2 signaling’ pathway is implicated in AD [47]. Similarly, the ‘mTOR (mammalian target of rapamycin) Signaling’ and ‘Regulation of eIF4 and p70S6K Signaling’ pathways are also involved in regulation of protein synthesis and degradation [48], and eukaryotic translation factor 4E (eIF4E) and p70S6K are downstream targets of mTOR [48]. Given the mechanistic interrelationships of these protein translation-related pathways and their statistically significant alteration in excitatory neuron cell types, it is possible that these pathways confer AD risk through their actions in these specific cell types. Many of the genes implicated in the ‘mTOR Signaling’ pathway are upregulated in 5XFAD mice, suggesting that increased activity of these genes could contribute to AD pathology. Interestingly, rapamycin, an mTOR inhibitor, is being investigated for use in AD treatment [49].
The canonical pathways, ‘Mitochondrial Dysfunction’ and ‘Oxidative Phosphorylation,’ both relate to nuclear DNA-encoded genes that function in mitochondria to generate energy for cellular functions, including synapse formation. Identification of these pathways as relevant in AD mouse models is consistent with previous studies showing that mitochondrial pathophysiology contributes risk to the development of AD [50]. The fact that both pathways are statistically significant in excitatory neuron subtypes implicates energy imbalance in these cells specifically.
Alteration of the pathway ‘Reelin Signaling in Neurons’ is an interesting finding because reelin plays a role in early development in the cerebral and cerebellar cortices [51] and changes in reelin expression have been observed in AD [51–53]. One study suggests that reelin expression is upregulated in the brains of AD patients due to impaired reelin activity that results from Aβ protein accumulation [52]. Another study reports that altered reelin signaling or processing might play a role in neuronal abnormalities in AD [53].
The canonical pathways, ‘LXR/RXR Activation’ and ‘FXR/RXR Activation,’ are related to receptors involved in cholesterol regulation. Liver X Receptors (LXRs) form heterodimers with Retinoid X receptors (RXRs), which are nuclear receptors that bind retinoids, compounds related to vitamin A. Retinoid signaling is active in the hippocampus and prefrontal cortex, where it is involved in neuronal plasticity [54]. Retinoids are thought to play a role in AD through their effects on oxidative stress (which may relate to mitochondrial dysfunction), Aβ accumulation, inflammation, and neurotransmission [54]. The pathways, ‘LXR/RXR Activation’ and ‘FXR/RXR Activation,’ were both implicated in excitatory neurons, inhibitory neurons, and resting astrocytes. Alteration of these pathways is statistically significant in both 5XFAD and T2KO mice, suggesting that perturbation of mechanisms related to LXR and RXR pathways pervade multiple cell types and may be of central importance in the pathophysiology of AD.
We identified statistically significant alteration of only six canonical pathways after multiple testing correction in T2KO mice, in contrast to sixteen found in 5XFAD mice. However, pathways altered in common in the two models relate to inflammatory processes, highlighting the importance of inflammation and immune response to the development of AD pathology. This view is supported by results of other studies of AD mouse models using RNAseq data [55]. We find these immune-related and inflammatory pathways to be particularly enriched in resting astrocytes, which is also consistent with findings from previous studies [29, 56].
A novel aspect of our study involved investigating cell type-specific DTU. Several gene peaks with DTU also overlapped between mouse models. All but one of the model-independent genes with DTU are found in subtypes of excitatory neurons. The exception is Stmn4 for which we detected DTU in MFOLs. Stmn4 encodes Stathmin 4, which is involved in destabilization of microtubules and is linked to AD [57]. This gene is expressed in various cell types including neurons, but the significance of an Stmn4 DTU specifically in oligodendrocytes is unclear. Srrm2 encodes Serine/Arginine Repetitive Matrix 2, which is involved in RNA splicing. It was detected in clusters of excitatory and inhibitory neurons in the WT versus T2KO comparison. Srrm2 has been implicated in AD, though additional studies are needed to determine the role Srrm2 plays in AD pathology [58–60].The model-independent gene with DTU detected in the greatest number of clusters is Rbm4b (found in excitatory neuron clusters). Rbm4b encodes RNA Binding Motif Protein 4B which plays a role in vitamin A signaling, and can also bind to TTR [61]. This result may be linked to the Ttr DEG detected in some of the same excitatory neuron clusters as discussed above. Retinol binding protein is a carrier protein that binds retinol and transports it in the bloodstream [62]. Retinol binding protein level is increased in brain [63] and decreased in cerebrospinal fluid [64, 65] of individuals diagnosed with AD.
The other gene with DTU that we find in subtypes of excitatory neurons in both mouse models is Opcml. This gene encodes opioid binding protein/cell adhesion molecule like, a member of the IgLON subfamily in the immunoglobulin protein superfamily of proteins that is localized in the plasma membrane. It may play a supporting role in opioid receptor function [66, 67]. Interestingly, there is evidence supporting a role for OPCML in cognitive function. Initial work showed that OPCML SNP rs563954 (A allele) associated with cognitive function measures in a GWAS of educational attainment [68], and subsequent work identified the C allele of OPCML SNP, rs16606197, as a top SNP associated with cognitive decline in older adults free of dementia [69]. More recent work shows that SNP-SNP interactions in this gene associate with levels of neurofibrillary tangles and total paired helical filament tau (PHF-tau), measures of neuropathology, in brain of patients with AD [70]. Studies also show it is a key gene in transcriptional pathways that are altered in brain of individuals with AD [71, 72], and it has potential value as an AD biomarker [73].
An important paralog of Opcml is Ntm, which we find as a DEG (in both mouse models) and also a gene with DTU (in 5XFAD mice), specifically in resting astrocytes. Ntm encodes neurotrimin, which like Opcml, is a member of the IgLON subfamily of immunoglobulin (Ig) domain-containing glycosylphosphatidylinositol (GPI)-anchored cell adhesion molecules. These genes are located on chromosome 11q25, a region notable for statistically significant genetic linkage in studies of AD and phenotypes related to cognitive function [74–76]. Genetic association studies of NTM also support its relevance to AD. A SNP in NTM (rs73035809) reached a genome-wide “suggestive” level of statistical significance for neurofibrillary tangle pathologic burden [77]. A different NTM SNP (rs1040103) also associated significantly with AD pathology [78]. In a third study, interactions between multiple pairs of SNPs in NTM showed statistically significant association with total PHF-tau measurements [70]. Our finding that Ntm, and its paralog Opcml, are both altered specifically in resting astrocytes in two mouse models of AD supports and further refines the potential role of these genes in brain pathology related to cognitive dysfunction in humans.
Carbajosa et al. (2018) [79] generated RNA-Seq data from cortex and hippocampus of T2KO and WT mice at 4 and 8 months of age. Gene network construction was conducted and modules of co-expressed genes were expressed within gene networks. When cross-referenced to our list of DTU genes, the module associated with “RNA splicing and cell-cycle related functions” contained the gene Zranb2, a significant DTU finding in the WT versus T2KO comparison in subtypes of excitatory neurons (Supplementary File 3). This is notable, as the A allele in ZRANB2-AS2 SNP, rs2163503, has been found to associate with a measure of cognitive function in a recent GWAS study [68]. However, we note that the age of the mice used, and the brain regions studied differ in Carbajosa et al. (2018) [79] from our study.
A recently published study that conducted a secondary analysis of RNAseq datasets of human AD brain identified 435 genes with DTU that overlapped in the two datasets in the temporal lobe samples which they analyzed [80]. When compared to the DTU list from our analysis, App was found in clusters of various subtypes of excitatory neurons in the WT versus 5XFAD comparison. Osbpl6 and Ttc14 were found in clusters of excitatory neurons in the WT versus T2KO comparison. However, there are certain caveats to this analysis as the methodology used to generate the data in the referenced study [80] differed from what we used in our analysis as it is unclear from what cell types these DTU originated in that study.
After IPA of gene lists containing DTU was performed, we identified statistically significant alteration of only fourteen canonical pathways after multiple testing correction in T2KO mice, compared to twenty-five in the 5XFAD mice. However, pathways altered in common in these models highlight synaptic transmission in AD pathogenesis. Notably, in 5XFAD mice, ‘Amyloid Processing’ and ‘Neuroprotective Role of THOP1 in Alzheimer’s Disease’ are altered in subtypes of excitatory neurons. The consequences of pathway changes related to DTU could be the focus of future studies.
Limitations
The results of this study provide new information relevant to the biological mechanisms underlying AD; however, there are several limitations. First, the slightly different genetic backgrounds of the WT, T2KO, and 5XFAD mice is an uncontrolled factor that could influence effects on DEGs or genes with DTU. Furthermore, we have n = 1 for the number of libraries analyzed for each genotype, which derives from the study that generated the original dataset [14]. An additional limitation is the inability to detect splice variants occurring more 3’ than the final 100 bases of the transcripts. Moreover, our data are correlative, as we did not conduct experiments to show causative effects. Another limitation is that we did not confirm the DEGs or genes with DTU, as this is beyond the scope of an exploratory study. A final limitation to note is the caveats associated with using mouse models for a human brain disease as species differences could mask important findings or identify effects in mice that are not relevant to humans.
Conclusion
Single-cell based analyses are useful for learning about the biological factors underlying AD as they promote a higher level of understanding than previously possible. IPA findings from our DEG lists highlight the importance of energy imbalance and inflammation in the development of AD. When those processes malfunction, perhaps in response to amyloid and/or tau deposition, this could promote the onset of AD. Our DTU analysis of mouse snRNAseq data implicate several genes which may play a role in AD. Single-cell (nuclei) findings, such as those from this study, have the potential to inform new treatments that target specific cell types found in the brain. DTU might be important in gene expression changes at the protein level. If the presence or absence of miRNA binding sites, mRNA stability elements, or RNA localization codes within the 3’UTR of transcripts is considered, one might expect differential protein levels in general and also at specific locations within neurons. Additional studies examining differential gene expression and transcript usage at the single-cell level are needed to identify specific cellular roles and processes in the pathogenesis of AD for the eventual development of novel ADtherapeutics.
Footnotes
ACKNOWLEDGMENTS
The preparation of this manuscript was supported by a National Institutes of Health grant R01 MH109260 (WHB is PI). AEW was supported by T32 MH014654 (WHB is PI) and by a Tobacco Settlement Act grant from the Pennsylvania Department of Health (WHB is PI). The authors wish to thank those providing support for this work. BCR was supported by a 2017 NARSAD Young Investigator Grant (#26634) from the Brain and Behavior Research Foundation as the Patrick A. Coffer Investigator, funding for which was generously provided by Ronald and Kathy Chandonais. The authors thank the researchers Yingyue Zhou, Marco Colonna MD, and their colleagues at Washington University in St. Louis for generating and graciously making publicly available the dataset GSE140399 which was used in our analysis.
