Single-Cell RNA-seq reveals transcriptomic modulation of Alzheimer’s disease by activated protein C

Integration and annotation

Generally following the Satija Lab integration vignette [19], integrated datasets were made utilizing four sample conditions from wild type C57BL/6 mice w/o the administrator of APC (WT, WT+APC), and 5xFAD C57BL/6 mice with the administrator of APC (5xFAD, 5xFAD+APC). These integrated datasets allowed for observing transcriptional differences between 5xFAD + APC, 5xFAD, and WT mice and 20 unsupervised clusters were generated (Figure 1A). Additionally, a control dataset was generated to compare WT + APC versus WT (Supplementary Figure 1A). The Seurat integration workflow generates a combined dataset that achieves the maximum overlap of cells from both samples to annotate cell-types [19]. This ensured that each cluster annotated contained cells from both samples within the integrated dataset as much as possible.

The top three highly expressed conserved genes in each unsupervised cluster were used for cell-type annotation (Figure 1B, 1C). Of the 20 unsupervised clusters, only 17 expressed at least three canonical marker genes and were retained for further downstream analysis. Astrocyte cell populations were annotated by high expression of Slc1a2, Plpp3, and Clu. Endothelial cell populations were annotated with remarkable expressions of Flt1, Cldn5, and Ly6a. Microglial cell populations were chosen based on expression of markers Btg2, Kctd12, and Cx3cr1. Neuron cell populations were annotated with expression of Gria2, Nrxn1, and Syne1, though overlap with astrocyte cells were noted. Oligodendrocytes were annotated by Enpp2, Plp, and Ptgds markers and OPCs were annotated on expression of Fabp7, Apod, and Ptn. The top marker gene can be visually observed in the dimensional projection (Figure 1D).

Cell population

After cell-type annotation, the proportion of cell-types that form each sample can be observed and provide insight on critical alterations in cell population between samples (Figure 1E). Notably, the proportion of microglia in 5xFAD mice is significantly greater than their WT counterparts, making up over half the cells recovered. Decreases in astrocyte and neuron populations are also observed in 5xFAD mice, both populations reducing to ~4%. After introduction of exogenous APC, a rebound of impacted cell-types can be observed. Astrocyte and neuron cell populations return to near-WT levels as well as an incredible decrease in microglia cell population in APC-treated 5xFAD mice. The astrocyte population returns to roughly 11% of cells and neurons reached their greatest proportion after APC treatment nearly 10%. Microglial cells are cut in half, decreasing from ~53% in 5xFAD mice to 24% in APC treated mice. Notably, oligodendrocytes are also at their largest proportion following APC treatment at ~35%.

Globally, some 59,000 cells are contained in the main integrated dataset, ~15,000 from WT, ~22,000 from 5xFAD, and ~22,000 from 5xFAD with APC treatment (Figure 1C). Of these cells, 6,254 astrocyte cells were identified, making up just over 10% of the overall dataset. 5,643 endothelial cells were identified at just under 10% globally. Most of the dataset is made up of microglial cells with 21,618 cells identified and making up ~37% of total cells. The least number of cells identified were neurons with 4,331 cells that form ~7% overall. Oligodendrocytes are the second most abundant cell-type with 16,480 cells identified making up ~28% overall. Finally, OPCs make up 8% of the total cells in the dataset with 4,715 cells identified (Figure 1E). In the control dataset, APC-treated WT mice do not exhibit drastic changes in cell populations (Supplementary Figure 1B). Astrocyte, endothelial cell, microglia, neuron, oligodendrocyte, and OPC populations remain largely unaffected by exogenous APC nullifying major adverse effects in the baseline. We observed 19 DEGs (Log₂FC > ± 0.3, adj. p-value < 0.05) globally between WT with APC treatment and WT (Supplementary Figure 1C).

Differential expression

Once stable cell-types were established, a wide variety of differentially expressed genes (DEGs) in all cell-types were found. The top five DEGs per cell-type were identified and cross-referenced with bulk-RNA seq (Figure 2A). We obtained the bulk RNA sequencing data from GSE140286. The log₂FC for each gene were calculated by comparison of 6-month 5xFAD to 6-month WT using Limma Package. We mapped genes in GSE140286 with the top DEGs list from scRNA sequencing. Comparing the DEGs between bulk and single-cell sequencing reveals that genes we observed to be upregulated in single cell-types in 5xFAD mice have been found to be generally upregulated in bulk sequencing (Figure 2A). We did not observe the DEGs that were downregulated in 5xFAD mice in our dataset to broadly match in the bulk sequencing data.

Consequently, we focused on genes that are greater expressed in 5xFAD compared to WT and later strongly downregulated by APC treatment (Supplementary Figure 2A, 2B). Generally, we highlighted exceptionally expressed genes meeting a threshold of Log₂FC values ~1 resulting to 77 total DEGs across all cell-types. Many DEGs are globally upregulated in all cell-types in 5xFAD compared to WT. We identified 46 global DEGs that are all remarkably expressed in all 5xFAD cell-types and downregulated with APC treatment (Figure 2B). Many of these globally expressed genes are associated with the transgenic nature of 5xFAD and the onset of AD, their roles in each cell-type are of particular interest. The implications of these global DEGs may be more impactful in a specific cell-type versus another such as Ctsb in microglia, further investigation is needed to verify the physiological impact of specific genes in AD treatment.

Within astrocytes, Gfap particularly is abundantly found in the astrocyte cell population and profusely upregulated in 5xFAD (Figure 2C). Compared to WT, AD astrocytes have a 2.68 log₂FC increase in Gfap expression; after APC treatment, Gfap became significantly downregulated with a 2.52 log₂FC decrease in expression (Figure 2D). Many global DEGs including B2m, C1qa, Ctsb, Trem2, and Tyrobp are all impressively down regulated in 5xFAD + APC (Figure 2B). We find the DEGs found in endothelial cells appear to be non-exclusive. Plxdc2 are broadly expressed in other cell-types; however, are only strongly differentially expressed in endothelial cells (Figure 2C). Plxdc2 receives a 1.11 log₂FC expression increase followed by a 0.91 log₂FC decrease in expression with APC (Figure 2D). Like found in astrocytes, Apoe, C1qa, Ctsb, Trem2, Tyrobp, etc. were found to be especially upregulated in endothelial cells of 5xFAD mice (Figure 2B). These genes are not necessarily specific to endothelial cells; however, the effects of these DEGs on endothelial cells may provide insight on neurodegeneration. These DEGs get significantly downregulated after APC treatment in endothelial cells. The importance of any single gene in endothelial cells specifically is not well understood.

Microglial cells are widely associated with the progression and pathology of AD. Many DEGs that are globally upregulated are the most implicated in microglia regarding AD. These include genes: Apoe, Cst7, Cd63, Cd68, cathepsins (Ctsa, Ctsb, Ctsd, ect.), Irf8, Trem2, and Tyrobp. Though most of the DEGs are modulated by APC treatment, interestingly, C1qa, C1qb, C1qc, Hexa/b, Laptm5 and Ly6e do not seem to alter in expression level extensively in microglia (Figure 2B). The microglial-exclusive DEGs most significantly modulated by APC were Aplp2, Axl, Baiap2l2, Cacna1a, Ccl3, Ccl4, Ccl6, Lpl, Cd74, Cox6a2, Csf1, Fabp5, Gnas, Hif1a, Ifitm3, Igf1, Itgax, Lilrb4a, Myo1e, Nceh1, Pld3, and Spp1 (Figure 2C). The downregulation of the cathepsin family of genes by APC treatment is distinctly remarkable (Figure 2C).

In neurons, like endothelial cells, the DEGs found are expressed in various cell-types though firmly differentially expressed in neurons (Figure 2C). Robust DEGs in 5xFAD neurons with a greater than ~1.4 log₂FC increased expression include Rpl35, Rpl36, and Rpl6 (Figure 2D). It was found that APC treatment downregulated all the mentioned genes in neurons with greater than a 1 log₂FC decrease in expression (Figure 2D). Oligodendrocyte DEGs in 5xFAD also widely consist of global upregulated genes (Figure 2B). The vast majority of global DEGs in oligodendrocytes were also found to be downregulated by APC treatment (Figure 2B). This may implicate the roles of these DEGs in the large oligodendrocyte proliferation found in APC treated 5xFAD mice (Figure 1D). Global DEGs were also found in OPCs including B2m, C1qa, Cst7, cathepsins, Syngr1, Trem2, and Tyrobp (Figure 2B). We found Aif1 and Itgb5 to be remarkably differentially expressed in OPCs of 5xFAD mice (Figure 2C) Measuring a ~1 log₂FC increased expression in 5xFAD mice compared to WT and a corresponding ~1 log₂FC decreased expression with APC treatment (Figure 2D). Accordingly, wide downregulation of global DEGs were also observed in OPCs of APC treated mice akin to oligodendrocytes (Figure 2B).

Microglial alterations

Sub-setting the microglial population allowed for further transcriptome analysis of the ~22,000 cells identified in the original dataset (Figure 3A). We discovered ~900 significant DEGs (adj. p-value < 0.05) between the samples (Supplementary Figure 2C). As with our comprehensive study, we targeted DEGs that were strongly impacted by APC-treatment. We identified 36 microglial DEGs with a Log₂FC > 1 and an adjusted p-value of 0 (Figure 3B, 3C). A large portion of these DEGs were previously identified in the comprehensive differential expression testing (Figure 2B, 2C) including Aplp2, Ccl4, Fabp5.

We measure that the top genes (Log₂FC >2) that were expressed in 5xFAD microglia were Apoe, Cst7, Lyz2, and Lpl (Figure 3C–3F). These were all strongly downregulated with APC treated microglia (Figure 3C) as observed in multiple cell-types during the comprehensive results (Figure 2B). Many of the broadly regulated DEGs found previously (Figure 2B) remain largely implicated in the microglial population. Major lysosomal factors Ctsb, Ctsd, and Ctsz were found to be significantly differentially expressed in microglia (adj. p-value = 0) of 5xFAD and APC treated mice. Our subset further uniquely found additional DEGs that were not found during the differential expression testing of all cell-types. Ank, Crlf2, Rpl35, and Rpl38 were differentially expressed in the microglial cell population but were not found during previous testing comprehensively.

GO term enrichment

The data generated by differential expression testing allowed for GO term enrichment. This enrichment data provides insight towards the roles of each cell-type in AD progression as well as during APC treatment. Universally, every enriched GO term for 5xFAD cells as a reciprocal in APC treatment (Figure 4A–4C). Our enrichment data shows 5xFAD oligodendrocyte and OPCs downregulate axonogenesis followed by counter upregulation in APC treated mice. Endothelial cells and microglia upregulated axonogenesis in addition to the oligodendrocyte and OPCs in APC treatment (Figure 4A). Central nervous system and nervous system development is downregulated in astrocyte, neuron, and oligodendrocytes in 5xFAD mice (Figure 4A). Corresponding upregulated nervous system development in APC treated mice can be found in every cell-type except for microglia. Downregulation of neuron generation is found in astrocyte and OPCs in 5xFAD mice and is regulated in endothelial cells, neurons, and OPCs of APC treated mice (Figure 4A). 5xFAD astrocyte and endothelial cells were found to downregulate blood-brain barrier (BBB) maintenance and compelling upregulation of BBB maintenance was found in the microglia, neuron, and oligodendrocytes of APC treated mice (Figure 4A). AD mice exhibit considerable neuron death (Figure 4A). Neuron development is downregulated in OPCs, and neuron death is upregulated in oligodendrocytes in 5xFAD mice. Endothelial cells, microglia, neurons, oligodendrocytes, and OPCs all are implicated in the rebound of neurons of APC treated mice. Neuron development is upregulated, and neuron death is downregulated with APC treatment (Figure 4A).

We found that astrocytes, endothelial cells, and neurons play roles in inflammatory response in 5xFAD mice. Furthermore, many interleukin processes are upregulated across all cell-types in 5xFAD mice (Figure 4B). These inflammatory processes are well implicated with AD progression and receive broad downregulation in all cell-types of APC treated mice. Interleukin-1 beta is upregulated in astrocyte, endothelial cells, and neurons of AD mice. Astrocytes, endothelial cells, oligodendrocytes, and OPCs of APC treated mice contribute to the downregulation of interleukin-1 beta production (Figure 4B). The production of interleukin-2 and interleukin-10 in 5xFAD mice is implicated with endothelial cells, neurons, and oligodendrocytes. The downregulation of interleukin-2 and interleukin-10 following APC treatment are handled by the same production cell-types in 5xFAD (Figure 4B). Every cell-type was found to participate in the production of interleukin-6 in 5xFAD mice and following APC treatment, every cell-type contributed towards the downregulation of interleukin-6 (Figure 4B).

Microglial cell activation in 5xFAD mice were found to be upregulated in the microglia and OPCs. It appears that every cell-type except oligodendrocytes contributes to downregulating the microglial activation found in APC treated 5xFAD mice (Figure 4C). Myeloid cell activation in 5xFAD mice were upregulated in neurons and oligodendrocytes and downregulated in astrocytes, endothelial cells, and oligodendrocytes in APC treatment (Figure 4C). Amyloid fibril formation and amyloid-beta clearance is upregulated in 5xFAD mice in microglia, neurons, and OPCs and microglia, oligodendrocytes, and OPCs respectively (Figure 4C). Focal adhesion in 5xFAD mice is downregulated in microglia; all cell-types except for astrocytes upregulate focal adhesion assembly in APC treatment (Figure 4C).

Administration of APC attenuates the Aβ burden and improves cognitive function in AD mice

To determine whether APC treatment can reverse the amyloid plaques aggravation in AD, the 5xFAD and C57BL/6 WT mice (10 weeks old) were treated with recombinant APC daily (100 μg/kg/d i.p.) or vehicle (saline) for 5 months. We examined whether administration of recombinant APC influences the development of Aβ pathology in AD. The results demonstrated that APC treatment effectively inhibited the Aβ burden in both hippocampus and cortex (Figure 5A). Compared with vehicle, APC treatment reduced load by 53% in hippocampus and 50% in cortex, it suggests that APC efficiently delay the development of Aβ pathology and amyloid angiopathy. Considering the Aβ aggregation devastating cognitive function, we further analyzed the effects of APC treatment on the hippocampal-dependent spatial memory function with radial arm water maze. The results showed that APC treatment normalized the performance of 5xFAD mice on spatial learning and memory ability in the radial-arm water maze test (Figure 5B). Consistent with Aβ pathology, APC treatment benefits the cognitive function in AD model mice. These data indicated that APC exerts potent neuroprotective activity against AD pathology.

Selection of individual and sample preparation

Both male and female C57BL/6 wild type mice and 5xFAD C57BL/6 mice were supplied from Jackson Laboratory (Bar Harbor, ME, USA). Three biological replicates (whole brain tissue) were used per sample. All animal protocols used for this study were approved by the Institutional Animal Care and Use committee of the University of South Florida and comply to the NIH guide for the care and use of laboratory animals.

At 2 months of age, recombinant WT murine APC (100 μg/kg) was injected daily via i.p. for 5 months, collecting brain tissue in 5xFAD and WT mice. Four brain tissue samples were subject to dissociation to be sent for sequencing: WT, WT with APC, 5xFAD, 5xFAD with APC. Mice’s brains were excised and rinsed with PBS, subsequently minced with Miltenyi Biotec Adult Brian Dissociation kit for mouse and rat following the manual protocol. First, mice brain samples were minced into 2 mm by 2 mm pieces with Enzyme Mix 1, after adding Enzyme Mix 2, samples were incubated using program 37C_ABDK_01 with gentleMACS Octo Dissociator. After the completion of the program, FBS was added to stop the reaction, and the suspension was filtered by a 70 μm MACS SmartStrainer on top of a Falcon 40 μm cell strainer. Cells were pelleted at 600 × g for 5 mins at 4ºC for debris removal. Second, cold PBS and Debris Removal Solution were used following the manual protocol to remove the debris. Third, Red blood cell lysis was carried out with PEB buffer and Red Blood Cell Lysis Solution to remove erythrocytes. At last, after samples were washed with cold PBS with 0.05% BSA buffer three times and filtered through Flowmi Top Strainer, single-cell samples were counted with an automated CellCounter.

Analysis of scRNA seq data

Cells were processed using the Chromium Next Gem Single Cell 3′ Reagent V3.1 kit from 10x Genomics. A total of 10,000 cells per sample were loaded into a Chip G for GEM generation. Reverse transcription, barcoding, complementary DNA amplification and purification for library preparation were performed according to protocol. Sequencing was performed on a NovaSeq 6000 platform (Illumina) targeting 100,000 reads per cell and received as Sanger/Illumina 1.9 encoded FASTQ files. Quality scores (>28) of each FASTQ file were confirmed by FASTQC 0.11.15 (Babraham Bioinformatics, Cambridge). Processing of raw FASTQ files for bioinformatic analysis was completed by the 10x Genomics Cellranger [67] platform using the 10x Genomics mm10 reference genome using default parameters for the count function for each sample. Raw and processed data can be accessed under the GEO accession GSE227157.

The Seurat 4.0 [68] package in Rstudio 2022.12.0 + 353 IDE (R Core Team, 2022 and Posit Team, 2022) was applied to normalize and dimensionally reduce data. Raw samples contained ~21,000 features. Samples were log-normalized and top variable features were determined by the “vst” selection method for the top 2000 variable features with 30 dimensions from the canonical correlation analysis (CCA). The Seurat integration protocol was selected for this study to compare datasets while retaining biological variability and minimizing technical differences. The various tools provided by Seurat are recommended to correct for batch effects in scRNA-seq data [69]. The integration protocol was used for “anchoring” cells between datasets (batches) utilizing mutual nearest neighbors (MMN) approach [19]. These “anchors” group cells based on similar expression profiles and therefore likely to be of the same cell-type to be identified after dimensional reduction and clustering. Integration code is available in Supplementary File 1 (Integration.txt).

Cell clustering

The integrated dataset contains ~24,000 features and the top 2000 variable features which were scaled using the default linear model and a PCA reduction. Based on the amount of standard deviation (Supplementary Figure 3A), reduction to 10 principal components were used for the remaining downstream analysis to capture most variation within the dataset. UMAP was run on the top 10 principal components (PCs) (Supplementary Figure 3B). Clustering resolution (0.3) was determined by visually examining the incoming proportion of cells as resolution increased (Supplementary Figure 3B) using Clustree 0.5.0 [70].

Cell type annotation

Utilizing the (FindConservedMarkers) function in Seurat, we identified conserved marker genes in each unsupervised cluster to use for cell-type annotation. This function outputs a data frame containing all genes expressed in each sample and their corresponding Log₂FC and adjusted p-values. Adjusted p-values are generated in Seurat by statistically comparing expression with other genes in the assay and performing Bonferroni correction [71]. These conserved genes were found to be expressed in each sample and the top three genes were used for annotation (Figure 1C).

Cross-analysis with the Human Protein Atlas [72] was used to verify marker genes for annotation. Clusters that did not highly express (Log₂FC >3) canonical marker genes for Astrocytes, Endothelial cells, Microglia, Neurons, Oligodendrocytes, or Oligodendrocyte Precursor Cells (OPC) were not included for further downstream analysis to retain accuracy of annotations. Cell-type annotation code available in Supplementary File 2 (Cell-Annotation.txt).

Differentially expressed genes analysis

The Seurat vignette for differential expression testing was followed for differentially expressed gene (DEG) identification. In our testing parameters, two samples within the dataset were isolated compared with each cell-type at a time. The RNA data assay was used in contrast to the integrated assay as the RNA assay retains the original expression profiles of cells prior to “correction” that occurs during the integration protocol [19]. The original sample identities of the cells were set as the identities in the (FindMarkers) function with cell-type set as the subset identity. This function outputs a data frame containing all differentially expressed genes between the two identities set and their corresponding Log₂FC and adjusted p-values. Adjusted p-values are generated in Seurat by statistically comparing expression with other genes in the assay and performing Bonferroni correction [71]. All genes possessing an adjusted p-value of greater than or equal to 0.05 were omitted from further downstream analysis. Lists of significant DEGs (Log2FC > ± 1, adj. p-value < 0.05) per cell-type available in Excel format in Supplementary Files 3 and 4 (DEG-Dataset 1_5xFAD vs. WT.xlsx, DEG-Dataset 2_5xFAD APC vs. 5xFAD Vehicle.xlsx).

Bulk RNA-sequencing data was retrieved from the GEO accession GSE140286 [73]. Utilizing data from 5xFAD and WT mice at 6 months of age that aligns closely to the 7-month age of our mice. Differential expression testing between 6-month-old 5xFAD and WT mice bulk RNA-seq data was completed using Limma 3.16 [74] in RStudio 2022.12.0 + 353 IDE (R Core Team, 2022 and Posit Team, 2022). Enrichment data was processed by EnrichR. A data frame containing all associated GO biological processes and corresponding adjusted p-values is produced from the data obtained from the differential expression testing performed previously. EnrichR produces Benjamini-Hochberg adjusted p-values [75, 76]. All GO terms possessing an adjusted p-value of greater than or equal to 0.05 were excluded and terms representing developmental alterations and immune responses in the brain were selected. Differential expression testing and enrichment code available in Supplementary File 5 (DE_Enrichment.txt).

Amyloid plaques deposit staining

Heparin IV for anticoagulation was given by intraperitoneal injection with 1,000 units/kg 10 min before the experiment [77]. 5xFAD and C57BL/6 WT mice (10 weeks old) were treated with recombinant APC daily (100 μg/kg/d i.p.) or vehicle (saline) for 5 months. The experimental mice underwent anesthesia with 2–3% isoflurane and 100% O₂. The mice were transcardially perfused with ice-cold PBS. Brains were rapidly removed and fixed with 4% paraformaldehyde overnight at 4°C. Subsequently, the hemisphere was subjected to dehydration with 10, 20, and 30% sucrose and then embedded in a cutting temperature compound (Tissue-Tek). The fixed brains were sectioned at 25 μm thickness setting on a cryostat and postfixes. After washing with PBS, the sections were blocked with 5% normal donkey serum (Vector Lab)/0.1% Triton-X/PBS for 1 h and incubated with mouse anti-human 6E10 amyloid plaque antibody (Biolegend, CA, USA) diluted in blocking solution overnight at 4°C. After three PBS washes, sections were incubated with secondary antibodies in diluted blocking solution for 1 h at room temperature. Each stained section was incubated in 500 μM of Tiofavin S (TS, Sigma-Aldrich, MO, USA) dissolved in 50% ethanol for 7 min for TS double-staining. Finally, sections were washed with PBS three times and mounted onto slides with mounting medium, and observed on an SP8 confocal microscope (Leica). ImageJ was used to quantify the amyloid plaques load in the hippocampus and cortex.

Radial arm water maze behavior test

The experimental group mice were put into the behavior room in darkness 1 h prior to the start of the experiments. For 15 consecutive days, the experiment was started around the same time, and each mouse was in the same order. This experiment was carried out in the dark. WT or 5xFAD mice were put into the water maze one at a time at the starting arm. Each day, the platform was placed at the end of the goal arm. For each trial, mice were placed at the center of the starting arm, facing forward. As it started swimming, the experimenter timed for 60 seconds and stopped timing once it reached and climbed up the platform. Errors were counted as: (1) each time the mouse entered the wrong arm, the experimenter gently grabbed its tail and pulled it back to its starting position; (2) if the mouse stayed in the center and did not enter any arm for 15 seconds; (3) if the mouse entered the goal arm without climbing up the platform for 15 seconds; (4) if the mouse did not reach the platform within 60 seconds, while only entering 1 or 2 arms continuously. After the mouse reached the platform, it was allowed to stay on the platform for 30 seconds to gain familiarity with the surrounding. For each trial, the time the mouse took to reach the platform and the number of errors it made were recorded. After the 4th trial, the mouse waited for 30 min to run the 5th trial. After each trial, the water was stirred to avoid the remaining scent that would affect the next mouse. Male and female mice were tested separately and were given 1 h in between to let the scent dissipate.