Elena Galea
email: elena.galea@uab.es

Laura D Weinstock

Raquel Larramona-Arcas

Alyssa F Pybus

Lydia Giménez-Llort

Carole Escartin

Levi B Wood
email: levi.wood@me.gatech.edu

George W Woodruff

Multi-transcriptomic analysis points to early organelle dysfunction in human astrocytes in Alzheimer's disease

Keywords:

The phenotypic transformation of astrocytes in Alzheimer's disease (AD) is still not well understood. Recent analyses based on single-nucleus RNA sequencing of postmortem Alzheimer's disease (AD) samples are limited by the low number of sequenced astrocytes, small cohort sizes, and low number of differentially expressed genes detected. To optimize the detection of astrocytic genes, we employed a novel strategy consisting of the localization of pre-determined astrocyte and neuronal gene clusters in publicly available whole-brain transcriptomes. Specifically, we used cortical transcriptomes from 766 individuals, including cognitively normal subjects (Controls), and people diagnosed with mild cognitive impairment (MCI) or dementia due to AD. Samples came from three independent cohorts organized by the Mount Sinai Hospital, the Mayo Clinic, and the Religious Order Study/Memory and Aging Project (ROSMAP). Astrocyte-and neuron-specific gene clusters were generated from human brain cell-type specific RNAseq data using hierarchical clustering and cell-type enrichment scoring. Genes from each cluster were manually annotated according to cell-type specific functional Categories. Gene Set Variation Analysis (GSVA) and Principal Component Analysis (PCA) were used to establish changes in these functional categories among clinical cohorts. We highlight three novel findings of the study. First, individuals with the same clinical diagnosis were molecularly heterogeneous. Particularly in the Mayo Clinic and ROSMAP cohorts, over 50% of Controls presented down-regulation of genes encoding synaptic proteins typical of AD, whereas 30% of patients diagnosed with dementia due to AD presented Control-like transcriptomic profiles. Second, down-regulation of neuronal genes related to synaptic proteins coincided, in astrocytes, with upregulation of genes related to perisynaptic astrocytic processes (PAP) and down-regulation of genes encoding endolysosomal and mitochondrial proteins. Third, down-regulation of astrocytic mitochondrial genes inversely correlated with the disease stages defined by Braak and CERAD scoring. Finally, we interpreted these changes as maladaptive or adaptive from the point of view of astrocyte biology in a model of the phenotypical transformation of astrocytes in AD. The main prediction is that early malfunction of the astrocytic endolysosomal system, associated with progressive mitochondrial dysfunction, contribute to Alzheimer's disease. If this prediction is correct, therapies preventing organelle dysfunction in astrocytes may be beneficial in preclinical and clinical AD.

Introduction

'Reactive' GFAP-overexpressing astrocytes are found in the vicinity of amyloid-β plaques in postmortem brains from patients with Alzheimer's disease (AD), in both the autosomal-dominant Alzheimer's disease (ADAD) and the sporadic late-onset (AD) variants, as well as in mouse models of ADAD. The robust morphological transformation of astrocytes points to major phenotypical and, hence, functional alterations in astrocytes in AD (reviewed in [START_REF] Perez-Nievas | Deciphering the astrocyte reaction in Alzheimer's disease[END_REF]).

However, despite non-negligible research in the last decade, knowledge about phenotypical alterations in human astrocytes from AD patients is still limited due, in part, to scarce and conflicting human astrocyte-specific omics data. Prior studies include transcriptomics of laser-microdissected astrocytes [START_REF] Sekar | Alzheimer's disease is associated with altered expression of genes involved in immune response and mitochondrial processes in astrocytes[END_REF][START_REF] Simpson | Microarray analysis of the astrocyte transcriptome in the aging brain: relationship to Alzheimer's pathology and APOE genotype[END_REF], co-expression-based gene clustering of whole-brain AD transcriptomes [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], and single-cell RNA sequencing (scRNAseq) analyses [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF][START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF]. Recently, the study of phenotypic changes of astrocytes in AD has been approached with network analysis of postmortem immunohistochemical data obtained from perusal of the literature [START_REF] Viejo | Systematic review of human post-mortem immunohistochemical studies and bioinformatics analyses unveil the complexity of astrocyte reaction in Alzheimer's disease[END_REF].

Laser-microdissected astrocytes either presented an insufficient number of differentially-expressed genes (DEG) for pathway analysis [START_REF] Sekar | Alzheimer's disease is associated with altered expression of genes involved in immune response and mitochondrial processes in astrocytes[END_REF], or normalization anomalies, as suggested by the finding that 98% of the DEGs were down-regulated [START_REF] Simpson | Microarray analysis of the astrocyte transcriptome in the aging brain: relationship to Alzheimer's pathology and APOE genotype[END_REF]. The problem may lie in the low RNA yields of laser microdissection, exacerbated by the poor RNA quality after long post-mortem intervals. As an alternative to laser microdissection, clustering statistics based on gene co-expression identified function-specific gene modules using 1647 whole-brain AD transcriptomes [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF]. Cellular localization of modules was established using cell-specific markers, and causality between nodes was inferred with Bayesian statistics. Using this approach, alterations of glutamate and amino acid metabolism in an astrocytic module of 260 genes were discovered to be well-ranked for causal relevance in AD, but no further changes were unraveled.

Finally, although scRNAseq improves the cellular resolution of transcriptomics to the point of unraveling populations within a given cell-type, two challenges arise in studying astrocytes from postmortem brain samples with this technique. One is the low numbers of astrocytes per patient being sequenced, which renders astrocytic datasets underpowered (discussed in [START_REF] Liddelow | Microglia and astrocytes in disease: dynamic duo or Partners in Crime?[END_REF]). The second challenge is to analyze a sufficiently large number of subjects to ensure that patient heterogeneity, a critical factor for drug development in AD [START_REF] Devi | Heterogeneity of Alzheimer's disease: consequence for drug trials?[END_REF], manifests itself. For example, [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] and [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF] examined only 48 and 12 subjects, respectively, including patients and Controls. Step 1-3: Genes from human-brain cortex were classified as cell-specific or enriched using RNAseq data from cells isolated from aged healthy human brain cortices (Step 1), the univariate cell-type enrichment τ scoring (Step 2), and hierarchical clustering (Step 3).

Step 4: The astrocytic and neuronal gene clusters were manually categorized into functions and subfunctions.

Step 5: Control AD and MCI whole-brain transcriptomes from the MtSINAI, MAYO and ROSMAP databases were hierarchically sorted according to the cell-type specific clusters.

Step 6-7: Alterations of astrocytic and neuronal functional categories in AD and MCI groups vs Controls were statistically established using gene set variation analysis (GSVA, Step 6), and principal component analysis (PCA,Step 7).

Step 8: To compare patients with mice, the published transcriptome of APP/PS1 astrocytes was re-analyzed using our functional categorization and GSVA. The figure in Step 1 was created using cell art adapted from Servier Medical Art (https://smart.servier. com/).

The main goal of this study was to extract consensus data of astrocytic changes from bulk transcriptomic data from three large independent clinical AD databases, encompassing 766 subjects. 'Consensus data' was defined as data detected in at least two out of the three databases. Identification of astrocyte-specific genes was optimized by taking a reverse approach as compared to previous studies: instead of determining a posteriori which self-organized gene modules might be astrocytic, we localized pre-determined astrocytic gene clusters in wholebrain AD transcriptomes with two-dimensional hierarchical clustering. The workflow is depicted in Fig. 1. First, we classified brain genes as cellspecific or non-specific using RNAseq data from astrocytes, neurons, microglia, endothelial cells and oligodendrocytes isolated from aged healthy human brains [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF], and a combination of two algorithms: hierarchical clustering, and a recently proposed univariate cell-type enrichment τ scoring [START_REF] Kryuchkova-Mostacci | A benchmark of gene expression tissue-specificity metrics[END_REF]. Second, astrocytic and, for comparison, neuronal genes, were manually annotated according to functions. Third, neuronal and astrocytic gene clusters were localized in three AD whole-brain transcriptomes generated by the Mount Sinai Hospital (MtSINAI) [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF] the Mayo Clinic (MAYO) [START_REF] Allen | Human whole genome genotype and transcriptome data for Alzheimer's and other neurodegenerative diseases[END_REF], and the Religious Order Study and the Rush Memory and Aging project (ROS-MAP), which also contains subjects diagnosed with MCI [START_REF] Bennett | Religious orders study and rush memory and aging project[END_REF]. These databases present the advantages of including large cohorts (78-219 subjects), which greatly facilitates the use of multivariate systems approaches, and of having been generated by two different approaches, microarrays and RNAseq, thus diminishing techniqueassociated bias. Finally, alteration of astrocytic and neuronal functional categories in AD and MCI groups vs Controls was statistically established using gene set variation analysis (GSVA), and principal component analysis (PCA).

We present the most comprehensive transcriptomic analysis of human AD astrocytes to date. As advocated in a recent consensus article about reactive astrocytes [START_REF] Escartin | Reactive astrocyte nomenclature, definitions, and future directions[END_REF], the impact of pathway alterations suggested by omics was interpreted in the context of astrocyte biology, instead of resorting to simplistic categorizations of astrocytes as 'neuroprotective' or 'neurotoxic' from the point of view of neurons. In this vein, our analysis points to dysfunction of the endolysosomal system/mitochondrion axis as a key driver of the phenotypical transformation of reactive astrocytes in AD.

Material & methods

Transcriptome datasets

Cell-type specific RNAseq datasets generated from healthy humanbrain cortical tissue [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF] were downloaded from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) under the accession number GSE73721. Data were prealigned to gene symbols by the authors using human genome version 19. Data consisted of transcriptomic data from brain cells isolated from temporal cortical lobes of juvenile (8-18 years old) and adult (21-63 years old) non-demented individuals. Transcriptomes from 12 astrocytic, one neuronal, five oligodendrocytic, three microglial, and two endothelial-cell samples from individual patients were included (Supplementary file 1, table 2). Principal component analysis (PCA, computed in R, R Foundation for Statistical Computing, Vienna, Austria) revealed that the transcriptomic data of astrocytes isolated from young (<47 years old) and old (>47 years old) were separated along principal component 1 (PC1, Supplementary file 1, table 3). Since the goal of the study was to identify cell-specific signatures in aged subjects, data from young individuals were discarded, retaining a total of six astrocyte transcriptomes from 47 to 63-year-old subjects, together with all 11 transcriptomes from the remaining cell types. A threshold was set for gene expression such that a given transcript should present a count >1 in reads per kilo base per million mapped reads (RPKM) in at least one sample. Genes with no counts >1 were discarded from the analysis. 11,077 genes were above this threshold (Supplementary file 1, table 4).

The MtSINAI database containing microarray data of samples from the prefrontal cortex of 101 Controls and 129 AD patients obtained from the Harvard Brain Tissue Resource Center was downloaded from GEO under accession number GSE44770 [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF]. The Rosetta gene identifiers were converted to the corresponding gene symbols using the Rosetta/Merck Human 44 k 1.1 microarray platform table under accession number GPL4372. When multiple Rosetta identifiers mapped to the same gene, the identifier with the greatest variance across samples was selected. The MAYO database, which contains data from the temporal cortex of 78 Control samples and 82 ADsamples, was downloaded from Synapse.org under Synapse ID syn3163039 [START_REF] Allen | Human whole genome genotype and transcriptome data for Alzheimer's and other neurodegenerative diseases[END_REF]. The data were downloaded as MayoR-NAseq_RNAseq_TCXCounts_normalized.tsv, pre-aligned to Ensemble Gene Identifier. Ensemble Gene IDs were then converted to HGCN gene symbols using biomaRt v. 2.32.1 in R. When multiple Ensemble Gene IDs mapped to the same gene symbol, the one with greatest variance was selected. The ROSMAP database [START_REF] Bennett | Religious orders study and rush memory and aging project[END_REF] contains data from the dorsolateral prefrontal cortex obtained from autopsied nondemented individuals or patients diagnosed with AD. The data were downloaded from Synapse.org under synapse ID syn8456629 as ROS-MAP_DLPFC_netResidualExpression.tsv. The dataset contained 200 Control cases, 157 MCI cases, and 219 ADcases. In the MCI group, there was no stratification into MCI and MCI-AD. The dataset as posted by the authors had been adjusted via voom (mean-variance modelling at the observational level) normalization to remove bias associated with batch number, RNA integrity, and other data acquisition and processing covariates. Missing transcript values in the dataset were imputed using a nearest neighbor averaging method using the impute package in R.

Univariate cell-type enrichment score

To determine the cell-specificity of human genes listed in the braincell transcriptomics database [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF], we used the τ method, reported to be particularly robust compared with other scoring approaches to identify tissue specificity of a gene among different magnitudes of expression and sizes of datasets [START_REF] Kryuchkova-Mostacci | A benchmark of gene expression tissue-specificity metrics[END_REF]. The τ enrichment score was computed as:

τ = ∑ n i=1 (1 -xi ) n -1 ; xi = x i max 1≤i≤n (x i )
where x i is the mean expression of the gene within the i th cell type and n is the number of cell types. We then defined the cell type enriched in the gene as that with the greatest average expression among all cell types. Top τ scores per gene are in Supplementary file 1, table 5. A score of τ ≥ 0.8 indicates that the gene is specific for a given cell type, and a score of 0.8 > τ ≥ 0.6 that the gene is enriched in that cell type. For simplicity, henceforth we refer to genes with τ ≥ 0.6 as 'enriched', although note that the 'enriched' pool also includes specific genes. Finally, τ < 0.6 indicates that the gene is non-specific, although it is assigned to the cell type with the greatest average expression.

Hierarchical clustering

The cell-type specific data from [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF] (Supplementary file 1, table 4) were z-scored for each gene across samples. Gene and sample hierarchical clustering were conducted using a Euclidian distance measure, and the average agglomeration method in R. Hierarchical clustering offers a clear visualization of data, facilitating a transparent analysis. A recent ranking of gene-module detection methods shows that agglomerative hierarchical clustering and WGCNA perform similarly [START_REF] Saelens | A comprehensive evaluation of module detection methods for gene expression data[END_REF]. It is worth nevertheless emphasizing that the primary statistical tool used in the present study was not hierarchical clustering, but gene set variation analysis (GSVA, see below). In hierarchical clustering, gene clusters associated with neurons and astrocytes were determined by cutting the dendrogram at a height of 3.76 to achieve that at least 70% of the genes in the neuron and astrocyte clusters have the highest τ for that cell type (results from hierarchical clustering are in Supplementary file 1, tables 6-8). Corresponding cell-specific gene clusters in whole tissue data sets (MtSINAI, MAYO, and ROSMAP) were identified by sorting the genes of wholetissue data according to cell-type specific gene clusters. Whole-tissue subject samples (columns) were then clustered within each Control, MCI, or AD cohort (Euclidian distance and Ward.D2 agglomeration) to identify molecularly defined sub-groups of samples.

Data normalization within each cell type

A recognized caveat of bulk transcriptomic analyses is that DEG may reflect changes in cell-type composition rather than real changes in cellspecific transcriptomes. This is relevant to AD because there is neuronal death [START_REF] Gomez-Isla | Neuronal loss correlates with but exceeds neurofibrillary tangles in Alzheimer's disease[END_REF], and microglial proliferation [START_REF] Marlatt | Proliferation in the Alzheimer hippocampus is due to microglia, not astroglia, and occurs at sites of amyloid deposition[END_REF]. Since neurons are the most predominant cell type in the brain (40% in the human neocortex [START_REF] Pelvig | Neocortical glial cell numbers in Alzheimer's disease. A stereological study[END_REF]), a decrease of the neuronal pool, with no change in the astrocytic pool [START_REF] Marlatt | Proliferation in the Alzheimer hippocampus is due to microglia, not astroglia, and occurs at sites of amyloid deposition[END_REF], arguably results in the artificial up-regulation of astrocyte genes. Several approaches have been employed to correct for cellular composition, such as the use of the so-called 'high fidelity' genes primarily expressed by one cell type [START_REF] Kelley | Variation among intact tissue samples reveals the core transcriptional features of human CNS cell classes[END_REF]. Here, we reasoned that reliance on a subset of highly cell-specific genes could lead to artifacts should their expression change in the course of the disease, because there is no assurance that they are stably expressed as house-keeping genes. Likewise, recent de-convolution approaches such as CIBER-SORTx and MuSiC proposed to identify cell type composition and contribution to bulk RNAseq or microarray datasets [START_REF] Steen | Profiling cell type abundance and expression in bulk tissues with CIBERSORTx[END_REF][START_REF] Wang | Bulk tissue cell type deconvolution with multi-subject single-cell expression reference[END_REF] did not work well with our data, for we did not find strong agreement among methods with regards to neuron and astrocyte densities in Control brains, and changes thereof in AD, when applied to our datasets. Therefore, we elected to remove gross variance due to changes in cell density by variance stabilizing normalization using limma in R, which assumes a so-called affine transformation of gene expression within each sample, consisting of an offset and scaling factor selected to minimize inter-sample variance. What variance stabilizing normalization does is to bring genes from different samples onto the same scale in a manner that is blind to the sample group/identity. Removal of gross variance is a well-established approach to normalize transcriptomic data, based on the premise that only a negligible proportion of genes change their expression following perturbations [START_REF] Lin | Model-based variance-stabilizing transformation for Illumina microarray data[END_REF]. We reasoned that a bulk shift in the distribution of gene expressions with respect to zero in a given cell population reflected changes in the proportion of that cell type among brain cell populations. Thus, correcting gene expression according to such shift would unmask the real changes in gene expression. Gross variance was separately corrected for astrocytic and neuronal gene clusters because their respective raw data presented an opposite shift (positive and negative) in the distribution of differential gene expressions in Controls vs AD samples (Results).

Cell function annotation and gene set variation analysis

Neuronal and astrocytic gene clusters generated by combining τ (Kryuchkova-Mostacci and Robinson-Rechavi, 2017) and hierarchical clustering of transcriptomes from cells from healthy brains [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF], were organized into functional categories and subcategories (Supplementary file 2). Categories were manually curated by one expert and cross-validated by another expert using information from Gene-Cards, perusal of Medline and PubMed, and open-source platforms such as GO, KEGG, and Reactome. Manual curation was indispensable for two reasons. First, over 60% of the genes in the astrocytic cluster were not annotated in open-source platforms. This means a substantial loss of usable gene data if we only relied on genes annotated in open-source platforms. Second, annotation in such platforms represents 'canonical' pathways with little consideration for the fact that different CNS cell types are molecularly and functionally distinct [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF]. This limitation is critical for astrocytes. For example, unlike neuronal compartments such as 'synapses', 'dendrites' and 'spines', a highly specialized astrocytic compartment termed 'perisynaptic astrocyte processes' (PAP), concentrating RNA transcripts related to glutamate and GABA metabolism, energy metabolism, as well as the ribosomal machinery to perform local RNA translation [START_REF] Sakers | Astrocytes locally translate transcripts in their peripheral processes[END_REF], is not a category in current KEGG and GO databases.

The genes were organized into general categories and subcategories using a mixed set of criteria including general pathways (e.g., carbohydrate metabolism or amino acid metabolism) and subcellular compartments with a clear functional specialization (e.g., mitochondria, peroxisome, or lysosome). Such flexibility allowed us to optimize gene inclusion, while generating gene sets with enough genes for statistical purposes. Although most genes were assigned to only one category, around 100 genes were assigned to two categories (in blue, Supplementary file 2, tab 'astrocyte cluster 4'). Most of those genes (96) were PAP genes, because we reasoned that their presence in PAP does not exclude relevance in other subcellular compartments. For example, ALDOC and ACO2 are, plausibly, ubiquitous enzymes relevant for general glycolysis and mitochondrial Tricarboxylic Acid (TCA) cycle throughout the cell. Likewise, genes related to 'TNFalpha signaling' are in 'stress responses'/'cytokines' as well as in 'gliotransmission', for TNFalpha has been described to modulate glutamate exocytosis [START_REF] Santello | TNFα controls glutamatergic gliotransmission in the hippocampal dentate gyrus[END_REF]. Overlap between groups was less than 5% except for 'PAP and gliotransmission': 17.8% genes in 'gliotransmission' were in 'PAP', and 10.3% of PAP genes were in 'gliotransmission'. For neurons, each gene was assigned to only one category.

Using our novel annotation database, the GSVA package in R was used to identify the enrichment of each gene set across all samples. GSVA is a generalized gene set enrichment method that detects variations of pathway activity over a sample population in an unsupervised manner [START_REF] Hanzelmann | GSVA: gene set variation analysis for microarray and RNA-seq data[END_REF]. Statistical differences in enrichment scores for each gene set between subject groups were computed by comparing the true differences in means against a null distribution obtained by permuting the gene labels and re-computing the GSVA 1000 times. False discovery rate adjusted p-values (qFDR) were computed using the method of Benjamini & Hochberg [START_REF] Hanzelmann | GSVA: gene set variation analysis for microarray and RNA-seq data[END_REF]. Gene sets with FDR adjusted p value qFDR <0.05 were considered significant.

Relationships between clinical stages and GSVA scores

The ROSMAP dataset includes Braak and CERAD (Consortium to Establish a Registry for Alzheimer's disease) pathological scores, APOE genotype, and sex for each sample. To determine the relationship between these phenotypes and the GSVA enrichment scores for each sample, we used the stats R package to fit linear models between enrichment scores and either Braak or CERAD for each gene set. The MAYO and MtSINAI databases do not include supplementary information per subject other than the clinical classification.

Figure preparation

The graphical abstract was made using BioRender.

Results

Hierarchical clustering identifies groups of cell-type specific genes

To identify human astrocyte-specific genes, we calculated the metric τ (Kryuchkova-Mostacci and Robinson-Rechavi, 2017) for every gene in the transcriptome data from isolated brain cells generated by [START_REF] Zhang | Purification and characterization of progenitor and mature human astrocytes reveals transcriptional and functional differences with mouse[END_REF] (Supplementary file 1, table 5 'τ scores'). Cell type-specific (τ ≥ 0.8) and cell type-enriched (τ ≥ 0.6) genes represented 30% and 60% of the 11,077 genes, respectively.

In parallel, because cell-specific functions are arguably performed by co-expressed genes [START_REF] Weirauch | Gene Coexpression Networks for the Analysis of DNA Microarray Data[END_REF], we reasoned that cell-type specific genes would be identified by clustering analysis. Thus, we applied hierarchical clustering (Methods) to the cell-type specific datasets, generating 194 clusters that contained 1-1651 genes (Fig. 2A; Supplementary file 1, table 6 'Hierarchical clustering' and table 7 'Cluster description'). We found that genes from the same cell type segregated together, defining clusters of highly co-expressed genes, plausibly representing cell-specific functions. The dendrogram height of hierarchical clustering was set to identify gene clusters that were highly enriched for individual brain cell types. The criteria for 'enrichment' were that the cluster contained more than 100 genes, of which over 70% corresponded to a unique cell type with τ ≥ 0.6 (Supplementary file 1, table 8 'cluster selection'). Five clusters fulfilled these criteria, corresponding to astrocytes (cluster # 4), endothelial cells (cluster # 5), neurons (cluster # 13), microglia (cluster # 17) and oligodendrocytes (cluster # 164). There were other clusters containing over 77% of total of neuronal (e.g., clusters # 10, 116) or microglial (e.g., clusters # 30, 15) genes (Supplementary file 1, table 8 'cluster selection'), but they were excluded from the analysis because less than 54% of their genes were enriched (τ ≥ 0.6) for that cell type, suggesting lower cell type specificity.

We continued with astrocyte cluster # 4 (1651 genes) and neuronal cluster # 13 (1258 genes) for two reasons: (i) a cluster size of over 1000 genes is sufficient for functional categorization and gene set analysis, and (ii) since molecular changes in neurons are better characterized in AD than changes in astrocytes, neurons served as a positive control to assure correct data mining.

In the neuronal cluster # 13, 99.9% of the genes were neuronal according to τ score, of which 54.1% of genes were neuron-specific (τ ≥ 0.8), meaning that they encode proteins that perform functions highly specific to neurons. Indeed, the top 10 ranked genes based on τ (τ > 0.98) were related to neurotransmission (SYNPR, GABRA1, CNR1, SYT1, GABRG2 and GAD1), regulation of membrane potential (KCN2), microtubules (INA), and cell adhesion (RELN). The remaining 45.9% (τ < 0.8) were neuron-enriched or bulk genes plausibly supporting neuronspecific functions, since they were co-expressed with neuron specific genes (Supplementary file 1, table 8 'cluster selection'; Supplementary file 2, table 2 'Neurons cluster 13').

In astrocyte cluster # 4, 80% of the genes were astrocytic according to τ (Supplementary file 1, tables 7 and 8 'Cluster description' and 'Cluster selection'; Supplementary file 2, tab 'Astrocytes cluster 4').

Further, 30% had τ ≥ 0.8, including genes that encode typical astrocyte proteins such as GFAP (τ = 0.95), glutamine synthase (GLUL, τ = 0.95), glutamate decarboxylase (GLUD2, τ = 0.91), glutamate transporters (EAAT1/SLC1A3, τ = 0.83), lactate dehydrogenase (LDHB, τ = 0.91), and enzymes involved in mitochondrial fatty-acid oxidation (ACADVL, τ = 0.92). Enriched genes with τ between 0.6 and 0.8, or bulk genes below 0.6, arguably necessary for astrocyte functions, are the aquaporins AQP1 (τ = 0.739) and AQP4 (τ = 0.744), and the chaperone clusterin CLU (τ = 0.373), which was identified by GWAS as a risk factor in AD [START_REF] Han | Analyzing 74,248 samples confirms the association between CLU rs11136000 polymorphism and Alzheimer's disease in Caucasian but not Chinese population[END_REF]. These findings indicated that the group of genes gathered in cluster #4 is relevant to astrocyte biology.

Identification of cell-type specific gene clusters in AD databases

We next searched for the astrocytic and neuronal gene clusters in three AD whole-tissue transcriptomic databases (heatmaps of z-scores in Fig. 2B). Unexpectedly, hierarchical clustering of subjects (columns) revealed heterogeneity in groups with the same clinical diagnosis. Subjects were divided into at least two broad molecular groups within each of the Control, MCI, and AD cohorts across the three databases, as shown by the division of the dendrograms into two large groups in each cohort. We refer to type 2 subjects as those showing blue z-scores within the neuron gene cluster, suggesting massive down-regulation of neuronal genes, as compared to type 1 subjects, in which genes with red z-scores predominated. [START_REF] Neff | Molecular subtyping of Alzheimer's disease using RNA sequencing data reveals novel mechanisms and targets[END_REF].

The most profound down-regulation of genes within the neuronal gene cluster was observed in AD2. Because down-regulation of neuronal genes due to loss of synapses and neuronal demise is a hallmark of AD [START_REF] Gomez-Isla | Neuronal loss correlates with but exceeds neurofibrillary tangles in Alzheimer's disease[END_REF], we reasoned that AD2 was advanced AD, and Control1 a bona fide control composed of subjects that were neither demented, nor preclinical or pre-symptomatic at the time of death. Because Control1 and AD2 appeared to be the extreme groups across the continuum of AD phenotypes, they were selected for further analysis.

Since our goal was to identify changes in astrocytic functions, it was critical to eliminate gross variability among samples within the astrocytic and neuronal gene clusters due to global changes in cell composition. To do so, we applied variance stabilizing normalization (Methods) to all samples. Prior to normalization, the distributions of gene-expression differences between AD2 and Control1 shifted to the left in the neuronal cluster, suggesting depletion of neurons in AD2, while the distribution of astrocytic genes shifted to the right, suggesting a relatively higher content of astrocytes in AD2. The normalization shifted the peaks of the distributions of differences in gene expression between AD2 and Control1 to zero, thus validating this step of data processing (Fig. 3A). After normalization, the global down-regulation of all neuronal genes in cluster # 13 in type 2 groups was attenuated, but intra-cohort heterogeneity persisted such that type 2 clusters still presented a prominent down-regulation of neuronal genes as compared to type 1 clusters (Fig. 3B and Supplementary file 3). In addition, new gene clusters previously masked were unraveled by hierarchical clustering. Thus, there were three neuronal (Fig. 3B, N-a, b, c), and four astrocytic gene subgroups (Fig. 3B,b,c,d) with opposite expressions patterns in type 1 vs type 2 subjects. For example, N-b genes were globally downregulated in AD2 and Control2, as compared to Control1 and AD1. Importantly, N-b contained genes related to synaptic function, including glutamatergic and GABAergic neurotransmission (e.g., GABRB3; GABRB2; GAD1; PDYN; SYN2; GABRG2; SYN1 NAPB; GRIA2; SLC17A7; GRIK2), supporting that this cohort is a bona fide AD. The normalized expression of cluster # 4 genes in the three AD databases is in Supplementary file 4. It is worth stressing that although changes in functions follow the same trend (up or down) in Control1/MC1/AD1 and Con-trol2/MCI2/AD2 subjects, there are quantitative differences within type1 and type 2 groups; however, for simplicity, we limited herein the statistical comparison of all the sub-cohorts vs Control1.

Because the ROSMAP database included clinical covariate data, such as CERAD and Braak scores, APOE genotype and sex, we examined whether type 1 or type 2 sub-cohorts were associated with specific covariates to gain insight into why they were molecularly heterogeneous. We found that Braak and CERAD scores were significantly higher in AD groups than in Controls (Fig. 3C), confirming the clinical diagnosis. However, type 1 and type 2 subjects with the same clinical diagnoses did not significantly differ in their pathological stages according to CERAD and Braak scoring; nor were type 2 groups more enriched in APOE4 subjects or females (Fig. 3D).

GSVA reveals altered functions in AD and MCI vs Control1

We next identified changes in astrocytic and neuronal functions using the normalized expression data. Because most of the astrocytic genes were not annotated in existing pathways in open-source databases, we manually curated the astrocytic cluster into 17 functional categories (gene sets) and 135 subcategories, and the neuronal cluster in 15 functional categories and 79 subcategories (Methods, Supplementary file 2).

GSVA was then used to compare the gene sets corresponding to different functional categories among cohorts, yielding an enrichment score for each subject and gene set. We started by comparing the extreme cases, AD2 and Control1, using two criteria. First, probability values were computed using a permutation analysis across genes (cut-off for statistical significance was set at FDR-adjusted q < 0.05, Methods). Second, we considered that a function was altered in AD2 as compared to Control1 if such alteration was detected in at least two out of the three databases (consensus criterion).

For neurons, seven functions were altered in AD2 vs Control1 (Fig. 4A, statistics in Supplementary file 5, extended graphs in Supplementary file 6): 'neurotransmission' and 'gene expression' were altered in the three databases, and 'synaptic plasticity', 'regulation of membrane potential', 'neural development', 'mitochondria' and 'stress response' in MAYO and ROSMAP, but not MtSINAI. 'Synaptic plasticity', 'neurotransmission', 'neural development' and 'membrane potential' were down-regulated, and 'gene expression', 'mitochondria' and 'stress response' up-regulated. In MAYO and ROSMAP, the altered functions in Control2 and MCI2 followed trends like in AD2, and AD1 and MCI1 like in Control1 (Fig. 4B). In MtSINAI, we found no statistically significant difference between Control1 and Control2, and AD groups were similar, although changes in AD2 were more dramatic, as compared to Control1 (Fig. 4B).

For astrocytes, five functions were altered in AD2 vs Control 1 (Fig. 5A, statistics in Supplementary File 5): 'PAP', 'plasticity', 'stress response' were up-regulated in three/three databases, 'mitochondria' was down-regulated in three/three databases, while 'endolysosomal system' was down-regulated in two/three databases. Because 'PAP' are specialized astrocyte-neuron contacts and, hence, this category includes many genes related to gliotransmission, and since 'gliotransmission' was highly significantly changed in ROSMAP (qFDR<0.001), and trending in MAYO (qFDR = 0.07), henceforth 'PAP' and 'gliotransmission' were considered together. As with neurons, Control2 and MCI2 showed the same direction of change as AD2 compared to Control1 (Fig. 5B). Likewise, the enrichment of these gene sets in AD1 resembled that of Control1, particularly in MAYO and ROSMAP, while AD1 appeared to be a milder version of AD2 in MtSINAI. A. Distributions of gene expression differences between AD2 and Control1 cohorts before (yellow) and after (blue) variance stabilizing normalization, yielding zerocentered distributions (Methods). B. Hierarchical clustering of Control (blue), MCI (beige), and AD (red) subjects from Fig. 2B using normalized data in rows. For both astrocytic and neuronal genes, subjects self-organized into at least two distinct subgroups (defined as 1 and 2) within each group with the same clinical diagnosis, while genes self-organized into up to four groups (Na-c and Aa-d). In MAYO and ROSMAP, clusters showed opposite directions of change in type 1 and 2 subjects. In MtSINAI, Control2 was like Control1, and AD1 like AD2, although with less pronounced changes. C. Violin plots of pathological stages according to CERAD and Braak scoring shows no significant difference in pathological stages between sub-groups within each cohort in ROSMAP (means ± SEM, one-way ANOVA with Tukey post-hoc). D. Distribution of APOE genotypes and sex in the six groups. While 3/3 was the predominant genotype across all groups, and 3/4 was more abundant in AD, there were no differences in frequency of APOE4 alleles between type 1 and 2 subgroups of subjects with the same clinical diagnosis. Likewise, females were more abundant than males in all groups, but not in type 2 subjects as compared to type 1 (p values determined by pairwise Fisher's Exact test with Bonferroni adjustment).

Co-variance of gene sets with one another and with pathological covariates

Correlations between altered functions and pathological stages. Since CERAD and Braak scores were available for the ROSMAP data, we examined relationships between general functions and pathological stages by regressing enrichment scores for each gene set and sample against the corresponding CERAD and Braak scores, regardless of clinical diagnoses (Methods). Data from all the individuals in the three groups (Control, MCI and AD) were used in the analysis. For neurons, we found no statistically significant correlations between Braak and CERAD stages for any of the seven top dysregulated categories described in (Fig. 4A and Supplementary file 7). The lack of correlation can be interpreted to indicate similar gene dysregulation in all stages, suggesting that neuronal dysfunction is an early event in AD pathology, such that it appeared in subjects with no (i.e., Control2) or incipient (i.e., MCI2) cognitive deficits, and low CERAD and Braak scores.

For astrocytes, of the five functions that were significantly altered in AD2 vs Control1 (Fig. 5A), the following ones significantly correlated with Break/CERAD according to p-values: 'mitochondria' (0.0000203/ 0.000742), 'stress response' (0.000742/0.00196) and 'plasticity' (0.0768/0.00967) (Supplementary file 8). In contrast, 'PAP', 'gliotransmission', and 'endolysosomal system' did not correlate with pathological stages, suggesting that, as we reasoned with neurons, they represent early manifestations of alterations in astrocytic functions. For example, there was a progressive decline in the expression of nuclearencoded mitochondrial genes from CERAD 1 to 3 and from Braak 1 to 6, while the expression of genes involved in 'endolysosomal system' was altered at early CERAD and Braak stages, before the detection of cognitive deficits in subjects.

Correlations among altered functions. We examined correlations among all neuronal and astrocytic gene sets by regressing enrichment scores for each gene set against each other gene set. Hierarchical clustering of correlation coefficients and adjusted p-values are in Supplementary file 9 for ROSMAP, MAYO, and MtSINAI. In neurons, the strongest direct correlations were among highly related functions such as 'membrane potential', 'neurotransmission' and 'synaptic plasticity'. Likewise, in astrocytes, the strongest correlations were between 'PAP' and 'gliotransmission'. The strongest anti-correlations point to intertwined changes between: (i) morphological remodeling of PAP and loss of synaptic functions, and (ii) organelle dysfunction and 'plasticity' (which is enriched in survival pathways, see below).

Alteration of functional subcategories in astrocytes

We next aimed to gain further insight into which factors drive the functional alterations in AD astrocytes by searching for functional subcategories with statistically significant differences in enrichment scores between AD2 and Control1 (qFDR<0.05) in at least two out of three of the databases (statistics in Supplementary file 10).

The category 'PAP' includes different categories related to morphology, metabolism, and gene expression. The consensus subcategory significantly upregulated was 'PAP morphology' (qFDR <0.001, AD2 vs Control1, in MAYO and ROSMAP). Because PAP are specialized in the functional refinement of adjacent synapses at a micro scale [START_REF] Sakers | Astrocytes locally translate transcripts in their peripheral processes[END_REF], increased expression of PAP-related genes in AD astrocytes may be a compensatory reaction to locally optimize compromised synaptic functions.

Subcategories in 'plasticity' include annotations associated with morphological remodeling (e.g., 'adherens junction', 'ECM', 'cytoskeleton', 'ciliogenesis', 'TGFbeta signaling', 'metalloproteases'), and pathways involved in brain development whose role in adult astrocytes is not well understood, although they are characteristically up-regulated in reactive astrocytes, perhaps as a mechanism of astro-or neuroprotection (e.g., 'axon outgrowth', 'synaptogenesis', 'WNT signaling', 'hippo signaling', 'smoothened signaling', 'hedgehog signaling', 'homeobox signaling', 'notch signaling', 'FGF signaling', 'pluripotency', and 'glial cell fate commitment'). Significantly dysregulated pathways in AD2 vs Control1 were 'ECM', upregulated in the three databases, and 'hippo signaling' and 'ciliogenesis' (i.e., the process of generation of a microtubule-based and centriole-derived cilium), respectively upregulated and down-regulated in two databases-note that the fact that the general trend of a given general category is towards up-regulation does not preclude that some subcategories are down-regulated.

'Stress responses' includes 'AMPK signaling', 'antioxidant response', 'cell death', 'chaperone', 'DNA repair', 'hypoxia', 'complement system' and 'cytokines'. No subcategory was dysregulated according to qFDR in at least two databases. However, it is worth noting genes related to TNFalpha in 'cytokines' were significantly dysregulated as the subcategory 'TNFalpha' of 'gliotransmission' in MAYO and ROSMAP, and trending in MtSINAI (qFDR = 0.093).

In 'endolysosomal system', encompassing genes related to 'endoplasmic reticulum', 'Golgi', 'lysosome' and 'vesicles', no single subcategories were significantly dysregulated, while in 'mitochondria', 'ETC' and 'TCA cycle' were significantly downregulated in AD2 vs Control1 in three and two databases, respectively. For 'ETC', qFDR was <0.001 in MAYO, < 0.001 in ROSMAP, and < 0.012 in MtSINAI. For 'TCA cycle', qFDR was <0.001 in MAYO, <0.39 in ROSMAP, and < 0.01 in MtSINAI. In total, the dysregulated subcategories in astrocytes point to alterations in mitochondrial respiration concomitant with morphological remodeling.

PCA confirms main functional changes in AD astrocytes

GSVA was dependent on gene sets organized into pre-determined functions. For this reason, we also used PCA as an alternative, unsupervised approach to gain independent insight into the predominant molecular changes of AD astrocytes. First, we examined whether expression of genes contained in the astrocytic cluster # 4 segregated cohorts in ROSMAP. PCA segregated the Control1, MC1 and AD1 groups to the left, and the Control2, MCI2 and AD2 groups to the right, according to the principal component 1 (PC1, horizontal axis, Fig. 6A). Violin plots of PC1 scores per group confirmed this tendency.

Second, we examined the functions performed by the proteins encoded by the top 50 dysregulated genes in PC1, representing the maximally co-varying astrocytic genes (Supplementary file 11). Mirroring our results with GSVA, genes related to 'plasticity' and 'PAP/ gliotransmission' were over-represented in the up-regulated portion, whereas genes related to 'endolysosomal system', 'gene expression', and 'signaling' were predominant in the down-regulated portion (Fig. 6B). Genes related to 'mitochondria' and 'stress responses' appeared in both sections, indicating mixed direction of dysregulation in these functions. Unlike GSVA, PCA also detected changes in 'signaling', as suggested by the presence of genes related to this function in the up and down segments of the top PC1 genes.

Third, we asked if the changes detected in AD astrocytes in ROSMAP with PCA were reproducible. Thus, we determined the capacity of the PC1 genes to segregate cohorts in MAYO and MtSINAI. As in ROSMAP, Control1 subjects were mostly negative, and AD2 cases were mostly positive in MAYO and MtSINAI, thus validating the discriminating capacity of the maximally co-varying genes in PC1 (Fig. 6C). AD1 segregated with Control1 in MAYO, but not in MtSINAI, where AD1 aligned with AD2, confirming less intra-cohort heterogeneity in MtSINAI than in the other two databases, as shown with GSVA. Altogether, these results point to multi-factorial changes in AD astrocytes, encompassing changes in astrocyte-neuron interactions and organelle dysfunction (model in Fig. 7A and description in Discussion).

Comparison with other human studies

Comparisons among whole-tissue transcriptomic analysis [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], scRNAseq studies [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF] [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] and the present study are in (Supplementary file 11). The greatest concordance was found between our study and [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF]. Thus, 60% of the genes of the astrocytic cluster in the Zhang et al. study were in cluster # 4. Importantly, 94.9% of the common genes were First prediction: Amyloid-β drives the phenotypical change of cortical astrocytes in AD by direct and indirect mechanisms. Direct: Amyloid-β causes endolysosomal and mitochondrial dysfunction in astrocytes, resulting in an adaptive astroprotective shift of gliotransmission to a GABA-predominant mode as a strategy to preserve mitochondrial respiration. Indirect: Synaptic dysfunction caused by amyloid-β in neurons augments the production of the PAP machinery and synaptogenic factors in astrocytes-a phenomenon of enhanced astrocytic plasticity to protect neurons. Second prediction: The phenotypical transformation of astrocytes in AD happens in early and late stages. Control2 subjects, cognitively normal, would be at stage a, and AD2 patients, diagnosed with dementia, at stage b. In stage a, the astroprotective change in glutamate/GABA fluxes may preserve neural-circuit homeostasis, because the GABA produced by MAO-B may locally compensate for the deficits in inhibitory tone that disrupts circuitry synchronization [START_REF] Lee | Slow wave sleep is a promising intervention target for Alzheimer's disease[END_REF]. By contrast, deficits in the endolysosomal system would result in partial phagocytosis of amyloid-β [START_REF] Sollvander | Accumulation of amyloid-beta by astrocytes result in enlarged endosomes and microvesicle-induced apoptosis of neurons[END_REF] and dystrophic neurites [START_REF] Gomez-Arboledas | Phagocytic clearance of presynaptic dystrophies by reactive astrocytes in Alzheimer's disease[END_REF], thus exacerbating amyloid-β accumulation, and hindering neuronal repair due to lack of debris elimination. In stage b, the mechanisms triggered to protect mitochondria may fail, and cortical astrocytes become deeply malfunctional, as suggested by the reduction in TCA cycle and ETC transcripts, despite the increase in protective stress responses and pro-survival pathways aimed at preventing astrocyte demise. For simplicity, only predominantly altered functions are depicted. ECM, extracellular matrix; ETC: electron transport chain, GLU: glutamate; oAβ: oligomeric amyloid-β; PAP: perisynaptic astrocyte processes; TCA: tricarboxylic acid. Figure adapted from Servier Medical Art (https://smart.servier.com/). B. GSVA of the transcriptome of astrocytes isolated from APP/PS1 mouse cortices using our manually curated functional categorization shows down-regulation of astrocyte-enriched functions. C. Left, GSVA of the transcriptome of APP/PS1 astrocytes using gene sets consisting of the top/bottom 50 dysregulated human astrocyte genes in PC1 from the ROSMAP database, and the top/bottom 50 DEGs from APP/PS1 astrocytes [START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF]. Human astrocyte-enriched genes dysregulated in AD (up and down) are down-regulated in APP/PS1 astrocytes, confirming B. Right, GSVA of PC1 AD and mouse APP/PS1 DEG gene sets in ROSMAP, MAYO and MtSINAI whole-brain transcriptomes shows that both the up and down-regulated mouse gene pools tend to be up-regulated in AD databases. Discussion of the result of this comparison can be found in the main text. concentrated in the UP portion of our PC1 list, and were particularly enriched in the 50 top UP PC1, which presented a 50% overlap with the Zhang et al. list. Common genes were related to 'gliotransmission' (e.g., glutamate transporter SLC1A3), 'plasticity' (e.g., hippo signaling TP53BP2), 'stress response '(e.g., antioxidant PON2), and 'mitochondria' (e.g., amino acid metabolism HADHB). By contrast, there was scarce representation of astrocytic genes from [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF] in the downregulated portion of our cluster # 4, where most of the genes related to mitochondrial metabolism and the endolysosomal system lie. Since the astrocytic cluster was generated by co-expression within AD transcriptomes in [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], we reason that down-regulated astrocytic genes may not co-express any longer with up-regulated ones, thus explaining why down-regulation of organelle-related genes was not detected in [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF]. In summary, our study confirms the relevance of glutamate and amino acid metabolism in AD pathogenesis reported in large-scale network analysis [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], and links it to mitochondrial dysfunction, which, we stress, presented a marked negative correlation with pathological stages according to regressions with Braak/CERAD.

As to scRNASeq studies, only 1.1% of Cluster # 4, and 0.3% of the top 50 PC1 genes, overlap with the 69 astrocytic DEG listed in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] (Supplementary file 11). Although the low number of DEG hinders statistically sound pathway analysis, it is worth noting that the top GO Biological Process of astrocyte DEGs in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] is ECM (GO:0030198), which coincides with our detection of consensus upregulation of the ECM functional category. Indeed, common genes between top 50 PC1 and [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] include PAP genes PREX2 and PLXNB1.

The overlap between the astrocyte cluster listed in [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF] and cluster # 4 is 12.7%, and 1.3% considering only the top 50 PC1 genes, although the glutamate transporter SLC1A3 is common to both lists. Considering functions, among the most dysregulated GO Biological Process in [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF], including 'central nervous system myelination' (GO:0022010), 'axonogenesis' (GO:0007409) and 'chemical synaptic transmission' (GO:0007268), only the latter two functions are equivalent to our detected up-regulation of developmental pathways and gliotransmission-related genes.

It is worth noting that the astrocytic DEG listed in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] and [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF] only share 32 genes, representing 47% of the DEG in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] and 4.5% of those in [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF] (Supplementary file 11). Moreover, only 2.3% and 14.6% of the astrocytic cluster from [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF] overlap with [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] and [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF].

The only partial equivalence among studies may be due to several factors, including: (i) differences in methodologies, (ii) classifications of disease (e.g., clinical in [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF] and our study, but associated to amyloid-β pathology regardless of clinical diagnosis in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF], (iii) brain regions sampled (prefrontal or temporal cortex in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF] [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF] and in our case, but entorhinal cortex in [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF]), and (iv) the total number of subjects analyzed (1647 in [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], 48 in [START_REF] Mathys | Single-cell transcriptomic analysis of Alzheimer's disease[END_REF], 12 in [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF], and 766 herein.)

Finally, the main consistent finding between our study and the network analysis of proteins detected immunohistochemically in AD brains [START_REF] Viejo | Systematic review of human post-mortem immunohistochemical studies and bioinformatics analyses unveil the complexity of astrocyte reaction in Alzheimer's disease[END_REF] is the detection of up-regulation of ECM and oxidative stress (included in 'stress responses' in our case'). As with the rest of studies, the main difference is that the study by [START_REF] Viejo | Systematic review of human post-mortem immunohistochemical studies and bioinformatics analyses unveil the complexity of astrocyte reaction in Alzheimer's disease[END_REF] does not include down-regulation of proteins related to endolysosomal systems and mitochondrial functions such as TCA and ETC associated with AD astrocytes, because their analysis is centered on upregulated proteins. Nor is subject heterogeneity reported in any of the studies discussed above.

Comparison with cortical APP/PS1 astrocytes

We concluded our analysis by asking whether changes identified in cortical AD astrocytes are recapitulated in cortical astrocytes isolated from 15-18-month-old APP/PS1 mice [START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF]. Comparison with recent snRNAseq mouse data [START_REF] Habib | Disease-associated astrocytes in Alzheimer's disease and aging[END_REF] was not appropriate because this latter study was carried out in the hippocampus.

We performed the following head-to-head comparisons between human and mouse data. First, we examined the changes of the functional categories identified in the human astrocytic cluster # 4 in APP/PS1 vs WT mice (n = 4 per group). GSVA revealed down-regulation of all functional categories (Fig. 7B). That is, APP/PS1 astrocytes mimic the down-regulation in 'endolysosomal system' and 'mitochondria' detected in human AD astrocytes, but not the up-regulation of 'PAP', 'stress responses' and 'plasticity'.

Second, we examined the expression of the top ranked down-and upregulated PC1 genes of human AD astrocytes in the mouse transcriptomes (gene sets 1 and 3, Fig. 7C-left panel), and the expression of the 50 most downregulated and upregulated DEG in APP/PS1 astrocytes [START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF] in the clinical databases (gene sets 2 and 4, Fig. 7Cright panel). Fig. 7C shows the hierarchical clustering of GSVAgenerated enrichment scores. To help interpret these comparisons, it is worth stressing that, in addition to a difference in species, human gene sets 1 and 3 consisted of astrocyte-enriched genes, while the mouse data consisted of the entire mouse astrocyte transcriptome, encompassing enriched and ubiquitous genes. Thus, gene set 2 is enriched in typical homeostatic astrocytic genes related to neuronal support (e.g., FZD10, HES5, SOX9), and gene set 4 in inflammation-related genes (e.g., TREM2, IL1β, CCL3, CCL4, CCL5, C14). These genes are predominantly microglial, according to τ scores, and except for C1S, did not coincide with functionally related human astrocyte-enriched genes in cluster # 4 subcategories such as 'cytokines' and 'complement' (Supplementary file 12).

Both human clusters segregated with the down-regulated mouse gene set 2 (Fig. 7C left panel), confirming the global down-regulation of astrocyte-enriched genes in APP/PS1 astrocytes, as opposed to the upregulation of generic inflammation-related genes. As also concluded in [START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF], this divergence points to a loss of homeostatic astrocytic functions and the adoption of a microglia-like phenotype, suggesting that cortical APP/PS1 astrocytes may undergo a phenotypical involution with time.

As expected, AD astrocyte gene sets 1 and 3 were respectively downregulated and upregulated in AD2 cohorts (Fig. 7C right panel), while the mouse gene set 2, containing neuron support genes, downregulated in APP/PS1 mice tended to be upregulated in AD2 cohorts. This might mean that the loss of homeostatic genes may be more advanced in APP/PS1 mice at 15-18 months of age than in AD2 patients. On the other hand, gene set 4 containing inflammation-related genes was largely upregulated in AD2, suggesting a predominant expression of typical microglial genes in AD2, plausibly upregulated in both microglia and astrocytes. Collectively, the comparison of human and mouse data suggests that there are both shared and distinctive changes in astrocytes in human AD and 15-18-month APP/PS1 mice. This factor should be taken into consideration when using animal models of AD.

Discussion

Because of the challenges associated with the isolation of astrocytes or their nuclei, and to profit from the wealth of large publicly available human-brain transcriptomes, the goal of this study was to gain insight into what happens to astrocytes in AD using systems-biology tools. The approach consisted of the re-compartmentalization of whole-brain transcriptomes from Control, MCI, and AD subjects into contributions from individual cell types by using pre-clustered cell-specific gene sets. The study has four main findings.

First, unbiased hierarchical clustering of patient cohorts reveals major intra-cohort heterogeneity, i.e., relationships are not robust between gene profiles, clinical diagnosis, and pathological stage according to CERAD and Braak scoring, as recently reported in [START_REF] Neff | Molecular subtyping of Alzheimer's disease using RNA sequencing data reveals novel mechanisms and targets[END_REF].

In other words, some non-demented subjects present AD-like astrocytes and neurons, while some demented subjects present healthy-like ones, particularly in MAYO and ROSMAP.

Second, GSVA and PCA show that AD astrocytes-so defined because they appear in patients with the typical down-regulation of neuronal synaptic functions documented in AD-undergo a profound phenotypical transformation encompassing up-regulation of genes related to 'PAP', 'gliotransmission', 'plasticity' and 'stress responses', and downregulation of genes encoding organelle proteins (e.g., nuclear-encoded mitochondrial genes related to TCA and ETC). The detection of downregulated astrocytic genes that may have gone unnoticed in large-scale clustering [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], for their expression change in opposite direction to most genes in the original cluster, is an asset of the study.

Third, according to regression analyses between functions and CERAD or Braak scores, some astrocytic functions may be altered early in the AD continuum (e.g., 'PAP', 'gliotransmission' and 'endolysosomal system'), while other functions may be increasingly dysregulated in parallel with the deposition of neuritic plaques and neurofibrillary tangles (e.g., mitochondrial functions, plasticity-related pathways, and stress responses).

Fourth, multivariate statistics and hierarchical clustering show that cortical astrocytes from human AD and from a 15-18-month-old mouse model of ADAD [START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF] are distinct transcriptional entities.

Limitations

A limitation of the study is the potential differences in gene expression due to different cortical regions being used in the three databases. However, the fact that consensus differences were found in astrocytic functional categories associated with AD supports the possibility that reactive astrocytes share common changes across brain areas, and that our method is robust enough to detect them. Likewise, our analysis with a single astrocytic gene cluster does not reveal astrocyte populations, as recently described in a mouse model of AD using snRNAseq [START_REF] Habib | Disease-associated astrocytes in Alzheimer's disease and aging[END_REF], although it is tempting to speculate that the four astrocytic gene sub-clusters identified in the clinical cohorts (Aa-d, Fig. 3) might represent distinct transcriptomic states. Arguably, future studies and approaches such as spatial transcriptomics will reveal local and regional differences of human astrocytes in AD.

Another limitation of our study is that pathway dysregulation during AD may result in the remodeling of cell-type specific gene clusters, such that gene clusters identified in healthy brain transcriptomes may no longer exist in AD. However, quantification of network reorganization in AD subjects using a metric called modular differential connectivity showed equal or enhanced connectivity in 95.5% of the modules [START_REF] Zhang | Integrated systems approach identifies genetic nodes and networks in lateonset Alzheimer's disease[END_REF], suggesting that clusters identified by co-expression analyses in healthy cells are, for the most part, preserved in AD. Nevertheless, we recognize that our analysis may overlook changes in transcripts ubiquitous in all cell types (e.g., glycolysis-related transcripts), as well as the up-regulation in astrocytes of transcripts typical of non-astrocytic cells, namely microglia, including the inflammatory response common to both glial cells described in 15-18-month-old APP/ PS1 transgenic mice [START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF]. Hence, our results may be just the tip of the iceberg in terms of defining the totality of changes in human AD astrocytes.

Because transcriptomic data are cross-sectional and descriptive, they do not demonstrate cause-effect relationships between pathway dysregulation and pathology, or unequivocally prove that changes are adaptive (meant to compensate for dysregulated functions and maintain or restore homeostasis), maladaptive (irreversibly contributing to dyshomeostasis and neurodegeneration), or epiphenomena with scarce bearing on disease onset and progression. Finally, genetic data does not directly inform about changes in proteins and metabolites. Hence, predictions, based upon these data need validation in appropriate models.

Key implications and predictions

Use of molecular data for subject stratification in clinical trials. A disconnect between presence of AD pathology and cognition has been long recognized (Neuropathology Group. Medical Research Council Cognitive, F., and Aging, S, 2001). Neurocognitive measurements, atrophy patterns assessed with neuroimaging, neurofibrillary tangles assessed postmortem, CSF levels of amyloid-β1-42 and Tau, as well as Tau neuroimaging, have revealed subtypes of AD [START_REF] Devi | Heterogeneity of Alzheimer's disease: consequence for drug trials?[END_REF][START_REF] Dujardin | Tau molecular diversity contributes to clinical heterogeneity in Alzheimer's disease[END_REF]. The clinical and biological heterogeneity in AD, which is increasingly being recognized as a major obstacle to establishing drug efficacy in clinical trials, starts to be addressed with molecular tools. An example is the recent detection of three molecular classes in an extended MtSINAI cohort (MMBB-AD), using RNAseq transcriptomics of four brain regions [START_REF] Neff | Molecular subtyping of Alzheimer's disease using RNA sequencing data reveals novel mechanisms and targets[END_REF]. They found the largest global change between AD patients and non-demented controls in the parahippocampal gyrus (PHG), with AD patients presenting three molecular subtypes termed A, B and C. Although our analysis was performed in cortical samples rather than the PHG, our AD2 group might correspond to their 'C typical' AD class, and our AD1 to their 'A atypical' AD class. The reason is that both AD2 and C present down-regulation of synaptic networks as compared to AD1 and A. No further comparison between our study and [START_REF] Neff | Molecular subtyping of Alzheimer's disease using RNA sequencing data reveals novel mechanisms and targets[END_REF]) can be performed with respect to astrocytes because, among the only three astrocytic genes reported to be upregulated in the 'C typical' AD class', only GNA12 appears in our astrocytic cluster # 4-although it is indeed upregulated in ROSMAP (Supplementary data 13).

What are Control2 and AD1? The terms 'resistance' and 'resilience' have been coined to describe two distinct clinical scenarios in AD consisting of avoiding pathology ('resistance') vs coping with pathology ('resilience') (reviewed in (Arenaza-Urquijo and Vemuri, 2018)). Thus, resistant subjects would be cognitively normal without significant amyloid-β or Tau pathology, despite their being at high-risk considering factors such as age or APOE genotype. By contrast, resilient subjects would remain cognitively normal despite significant AD pathology. We speculate that our AD1 subjects were perhaps 'resilient' to AD pathology, which they had, according to high CERAD and Braak scores, but that they developed dementia due to comorbidities (e.g., Lewy bodies, vascular disease, TDP43 deposits). However, the possibility exists that the tissues used for transcriptomic analysis and neuropathological assessment of neuritic plaques and neurofibrillary tangles were not the same. This means that the transcriptomic analysis in demented AD1 subjects with a healthy-like molecular profile might have been performed in cortical tissue spared from ongoing pathology, such that AD1 subjects were not truly resilient to AD pathology. By contrast, regardless of where omics and neuropathological assessments were performed, Control2 subjects with AD-like molecular profiles in ROSMAP and MAYO were cognitively normal. Thus, Control2 subjects might be highrisk, pre-symptomatic, resistant subjects in whom the appearance of neuritic plaques, Tau tangles, and cognitive impairment was delayed despite the decreased expression of synaptic genes, plausibly caused by the incipient accumulation of soluble oligomeric amyloid-β [START_REF] Marsh | Synaptic dysfunction in Alzheimer's disease: the effects of amyloid beta on synaptic vesicle dynamics as a novel target for therapeutic intervention[END_REF]. Whatever the case, our study supports the use of omicsbased molecular phenotyping and clustering statistics as unbiased tools to stratify patients in clinical trials in the spirit of personalized medicine.

Early dysregulation of the endolysosomal system in astrocytes. The category 'endolysosomal system' was found down-regulated in AD2 vs Control1 by GSVA. Further, genes related to endoplasmic reticulum (e. g., FAF2, PIGX), endocytosis (e.g., ATXN2, MEGF8, RAB11FIP2), and autophagy (ATG9A) ranked among the top 50 PC1 downregulated in ROSMAP (Fig. 6B), and were significantly down-regulated in at least two of the three databases according to the Wilcoxon test (Supplementary file 4). These findings support cumulative evidence of inefficient phagocytosis and degradation of amyloid-β protofibrils [START_REF] Funato | Astrocytes containing amyloid beta-protein (Abeta)-positive granules are associated with Abeta40-positive diffuse plaques in the aged human brain[END_REF][START_REF] Sanchez-Mico | Amyloid-beta impairs the phagocytosis of dystrophic synapses by astrocytes in Alzheimer's disease[END_REF][START_REF] Sollvander | Accumulation of amyloid-beta by astrocytes result in enlarged endosomes and microvesicle-induced apoptosis of neurons[END_REF] and dystrophic neurites [START_REF] Gomez-Arboledas | Phagocytic clearance of presynaptic dystrophies by reactive astrocytes in Alzheimer's disease[END_REF] in astrocytes in mouse and in vitro models of AD. The decreased production of exosomes from astrocytes exposed to amyloid-β [START_REF] Abdullah | Amyloid-beta reduces exosome release from astrocytes by enhancing JNK phosphorylation[END_REF]) also points to malfunction of the astrocytic endolysosomal system in AD models. Thus, our study and others do not support the popular notion that astrocytes may efficiently clear amyloid-β plaques, as initially suggested by evidence in non-primary cultured astrocytes [START_REF] Wyss-Coray | Adult mouse astrocytes degrade amyloid-β in vitro and in situ[END_REF], which, arguably, are endowed with greater motility and morphological plasticity than primary or in situ astrocytes. Moreover, targeted enhancement of lysosomal biogenesis by overexpression of the transcription factor EB (TFEB) in astrocytes enhances Tau clearance, and reduces pathological hallmarks in the hippocampus of PS19 tauopathy mice [START_REF] Martini-Stoica | TFEB enhances astroglial uptake of extracellular tau species and reduces tau spreading[END_REF], suggesting that the improvement of lysosomal function in astrocytes may reverse Tau pathology in AD. Along these lines, TFEB has been reported to regulate AD risk genes in an astrocyte subpopulation detected in AD brains by snRNAseq [START_REF] Grubman | A single-cell atlas of entorhinal cortex from individuals with Alzheimer's disease reveals cell-type-specific gene expression regulation[END_REF].

Early and concomitant alteration of neuron-astrocyte contacts. This prediction is based on two findings: the inverse correlation of changes in synaptic functions and PAP, and the lack of correlation between these functional categories with pathological stages, suggesting that they change early in AD. Accordingly, both PAP and synaptic genes are dysregulated in non-demented subjects with a type 2 molecular profile (Control2), which, as noted earlier, may represent preclinical AD. Specifically, GSVA showed up-regulation of 'PAP morphology', whose leading gene EZR encodes ezrin, a protein involved in PAP motility, and in the anchoring of the astrocyte membrane to neurons or the extracellular matrix [START_REF] Derouiche | Perspectives for Ezrin and radixin in astrocytes: kinases, functions and pathology[END_REF]. Also, genes related to integrin and cell motility (e.g., FERMT2, PREX2) were part of the top genes detected by PCA as defining AD astrocytes (Fig. 6B). We posit that up-regulation of PAP genes is an adaptive change to improve the coverage of neurons to counteract the down-regulation of synaptic elements directly caused by oligomeric amyloid-β [START_REF] Marsh | Synaptic dysfunction in Alzheimer's disease: the effects of amyloid beta on synaptic vesicle dynamics as a novel target for therapeutic intervention[END_REF]. This scenario is further supported by the up-regulation of genes encoding for synaptogenic factors such as thrombospondin and glypicans (THSD1, GPC4) in the three databases (Supplementary file 4).

Mitochondrion dysfunction in AD astrocytes may impact gliotransmission. It is generally assumed that the well-established mitochondrion dysfunction [START_REF] Swerdlow | Mitochondria and mitochondrial cascades in Alzheimer's disease[END_REF] and reduced glucose metabolism in AD (often referred to as 'hypometabolism' [START_REF] Hoffman | FDG PET imaging in patients with pathologically verified dementia[END_REF]) are neuronal phenomena. However, our analysis reveals down-regulation of nuclear-encoded mitochondrial genes encoding TCA cycle and ETC components in astrocytes, suggesting deficiencies in energy-generating mitochondrial pathways. For example, one of the consensus downregulated genes, ATPAF1, is a component of the ATP synthase (complex V), and it is known that impairment of the ETC beyond complexes I-III increases reactive-oxygen species generation while decreasing ATP production [START_REF] Shi | Oxidative stress and transcriptional regulation in Alzheimer disease[END_REF]. In vitro and ex vivo studies link mitochondrial dysfunction in astrocytes to amyloid-β. Thus, in co-cultures of rat hippocampal neurons and astrocytes, amyloid-β causes loss of mitochondrial potential and transient mitochondrial depolarization, concomitant to activation of NADPH oxidase and decrease glutathione production, in astrocytes, but not neurons [START_REF] Abramov | β-Amyloid peptides induce mitochondrial dysfunction and oxidative stress in astrocytes and death of neurons through activation of NADPH oxidase[END_REF]. The latter study also shows that loss of membrane potential and depolarization of astrocytic mitochondria is prevented by antioxidants, and by the administration of glutamate as a substrate for the mitochondrial complex I. This suggests that amyloid-β-elicited mitochondrial dysfunction in astrocytes is caused by oxidative stress and deficits in substrate supply, such that increasing the supply of fuels amenable for mitochondrial oxidation including aromatic amino acids such as glutamate may be astroprotective. Based on this evidence, we interpret the consensus up-regulation of the mitochondrial glutamate transporter SLC25A18, and the plasma membrane glutamate transporters SLC1A3, SLC7A11 and SLC6A11 (Supplementary file 4) to indicate that AD astrocytes may take up glutamate to counteract the damage caused by amyloid-β to mitochondrial respiration. Along these lines, HADH8, which encodes an enzyme involved in mitochondrial β-oxidation, was found upregulated in AD astrocytes by PCA.

Likewise, the consensus up-regulation in AD astrocytes of the monoamine oxidase MAO-B, a dopamine-degrading enzyme located in the outer mitochondrial membrane, may be interpreted as an attempt to reverse respiratory deficits. This is supported by the recent discovery that MAO-B increases the polarization of the inner mitochondrial membrane and ATP production by shuttling electrons through the intermembrane space [START_REF] Graves | Dopamine metabolism by a monoamine oxidase mitochondrial shuttle activates the electron transport chain[END_REF]. MAO-B up-regulation in AD is relevant for two reasons. First, the MAO-B ligand (11)C-deuterium-Ldeprenyl may serve as a PET-based biomarker of early astrocyte dysfunction in AD, as documented in ADAD [START_REF] Rodriguez-Vieitez | Diverging longitudinal changes in astrocytosis and amyloid PET in autosomal dominant Alzheimer's disease[END_REF]. Second, MAO-B may affect the excitatory/inhibitory balance of neural circuits by producing GABA from putrescine, as shown in APP/ PS1 mouse hippocampus [START_REF] Jo | GABA from reactive astrocytes impairs memory in mouse models of Alzheimer's disease[END_REF]. In support of this scenario, the genes encoding the GABA transporter SLC6A11/GAT3, which may extrude GABA coupled to glutamate uptake, and SLC7A2, which imports the putrescine precursor L-arginine, are consensus dysregulated genes (Supplementary file 4).

Taken together, the data support that impaired mitochondrial respiration in astrocytes switches gliotransmission to a GABA-dominant mode. What may the clinical impact of this switch be? A detrimental effect of astrocytic GABA has been proposed based on the beneficial effects of the MAO-B inhibitor selegiline on electrophysiological readouts in APP/PS1 mouse hippocampi [START_REF] Jo | GABA from reactive astrocytes impairs memory in mouse models of Alzheimer's disease[END_REF]; however, there is no benefit of selegiline in patients [START_REF] Birks | Selegiline for Alzheimer's disease[END_REF], thus raising doubts as to whether APP/PS1 mice appropriately model the impact of astrocytic GABA in AD. Alternatively, astrocytic GABA might counteract the loss of GABAergic tone that has been causally related with neuronal hyperactivity, desynchronization of neural circuits, and amyloid-β production [START_REF] Lee | Slow wave sleep is a promising intervention target for Alzheimer's disease[END_REF]. Since our regression analyses reveal that mitochondrial dysregulation in astrocytes progressively increases in advanced pathological stages, and since the up-regulation of MAO-B detected by PET in prodromal ADAD is transient [START_REF] Rodriguez-Vieitez | Diverging longitudinal changes in astrocytosis and amyloid PET in autosomal dominant Alzheimer's disease[END_REF], we reason that the two strategies used by the astrocytic mitochondria to preserve membrane potential and ETC (i.e., enhanced supply of substrates such as aromatic amino acids like glutamate, and of electrons via MAO-B) eventually fail. As a consequence, glutathione production might be decreased and NADPH oxidase activity increased [START_REF] Shi | Oxidative stress and transcriptional regulation in Alzheimer disease[END_REF], perhaps exacerbating oxidative stress and overall astrocyte malfunction in advanced AD.

Pathological-stage dependent up-regulation of 'plasticity' and 'stress response' pathways might explain why, unlike neurons, astrocytes survive in AD. The progressive up-regulation of the subcategory 'ECM', which includes integrins and proteoglycans, is not surprising considering the indisputable adoption of a 'reactive' morphology by astrocytes in AD, consisting of engrossment of primary and secondary processes due to over-expression of the intermediate filament protein GFAP, and increased production of ECM proteins [START_REF] Escartin | Questions and (some) answers on reactive astrocytes[END_REF]. More intriguing is the consensus up-regulation of 'hippo signaling', an evolutionarily conserved network with a central role in the regulation of cell proliferation, cell fate, organ growth and regeneration [START_REF] Misra | The hippo signaling network and its biological functions[END_REF], and no reported role in healthy adult astrocytes, despite the fact that some of the genes in this pathway are highly specific to human adult astrocytes according to τ (e.g., WWC1 τ = 0.91; WWOX τ = 0.86, and YAP1 τ = 0.88). YAP1 was, moreover, detected by PCA as a relevant gene in the astrocytic phenotype in AD. The up-regulation of 'hippo signaling' in AD astrocytes, as well as the global up-regulation of genes involved in brain development may be a manifestation of the well documented phenomenon of re-activation of developmental pathways in reactive astrocytes that some authors interpret as a failed attempt at reprogramming into neural stem cells [START_REF] Torper | Brain repair from intrinsic cell sources: turning reactive glia into neurons[END_REF]), but we interpret to be a protective response to facilitate the survival of defective astrocytes. In addition, the recent observation that YAP1 regulates scarborder formation in spinal cord injury in mice [START_REF] Xie | Astrocytic YAP promotes the formation of glia scars and neural regeneration after spinal cord injury[END_REF], suggests that hippo signaling is coordinated with the increased production of ECM in AD astrocytes. Finally, the consensus dysregulation of the subcategory 'TNF-alpha' (Supplementary file 10), and the consensus upregulation of members of 'TNFalpha' TNFRSF11B, TNFSF13 and TRA-F3IP2, as well as of members of the IL17 family (Il17D, IL17RB, and IL17RD), and of complement factors (CD59, C1S) (Supplementary file 4), might be interpreted as detrimental 'neuroinflammation', although evidence also supports modulatory or protective actions of these factors [START_REF] Masgrau | Should we stop saying 'Glia' and 'Neuroinflammation'?[END_REF]. First, TNFα may facilitate the release of glutamate by improving vesicle docking at the astrocyte plasma membrane prior to exocytosis [START_REF] Santello | TNFα controls glutamatergic gliotransmission in the hippocampal dentate gyrus[END_REF], such that TNFα signaling may help adjust gliotransmission in AD astrocytes. Second, the finding that Nrf2/IL17D axis acts as an antioxidant protective pathway in stress responses induced by tumorigenic stimuli and viral infections [START_REF] Saddawi-Konefka | Nrf2 induces IL-17D to mediate tumor and virus surveillance[END_REF], supports an astroprotective role of members of the IL17 family in AD. Also, members of the complement system have been shown to facilitate amyloid-β phagocytosis [START_REF] Iram | Astrocytes from old Alzheimer's disease mice are impaired in Abeta uptake and in neuroprotection[END_REF].

Model of the phenotypical transformation of cortical astrocytes in AD

The key predictions from our analysis are summarized in Fig. 7A. We posit that amyloid-β drives the phenotypical change of astrocytes in AD in several stages by causing dysfunction of the endolysosomal/mitochondrial axis, which, in turn, prompts a change in the balance of excitatory/inhibitory neurotransmission. Impairment of mitophagy due to endolysosomal malfunction may further exacerbate mitochondrial dysregulation. Control2 subjects (cognitively normal/AD-like molecular profile) would be at stage a, and AD2 patients (demented, AD molecular profile) at stage b. In both a and b, the phenotypical change is complex and mixed with regards to possible clinical impact, for changes that may disrupt the homeostasis of astrocytes themselves and/or astrocyteneuron communications coexist with changes that may protect astrocytes and/or their interactions with neurons. The model thus represents a departure from simplistic neuroprotective/neurotoxic classifications of reactive astrocytes from the point of view of neurons. As recently discussed [START_REF] Escartin | Reactive astrocyte nomenclature, definitions, and future directions[END_REF], an 'astrocyte-centric' stance focused on the clarification of the maladaptive or adaptive nature of pathway changes in the complex astrocyte biology should be adopted. What is then the net impact of the mixed functional changes of reactive astrocytes in AD? In short, we posit that, despite organelle dysfunction, astrocytes manage in early disease stages to partially preserve their functions, including the modulation of neural circuits, through adaptive changes, while they become globally malfunctional in advanced stages. Thus, the early changes in astrocyte-neuron interplay may delay the onset of clinical symptoms, while detrimental phenomena would prevail in stage b. If these predictions are correct, therapies aimed at protecting and restoring the functions of the endolysosomal system and mitochondria to halt the transformation of astrocytes to stage b might help arrest the progression of AD. Even if GABA were detrimental, as reported in [START_REF] Jo | GABA from reactive astrocytes impairs memory in mouse models of Alzheimer's disease[END_REF], preservation of mitochondrial respiration is therapeutically indicated because increased GABA production, whatever its effect, would be the result of mitochondrial impairment. In summary, the present study points to prevention of organelle dysfunction in astrocytes as a therapeutic strategy in AD. These predictions need to be tested in transgenic mice or human cellular models that appropriately recapitulate the complex transformation of astrocytes in the AD continuum.

Conclusions

Our systems-biology based identification of astrocytic clusters in three independent datasets was aimed to extract consensus data of the main pathway changes in reactive astrocytes in AD. Arguably, the large number of subjects analyzed allowed us to detect molecular heterogeneity of subjects with the same clinical diagnosis by hierarchical clustering. We highlight that manually curated functional annotations were implemented to circumvent deficits in astrocyte-specific annotations in current platforms. The findings, supported by two independent statistical analysis, have led to a model of stage-dependent astrocyte dysfunction in AD caused by amyloid-β-elicited damage of the endolysosomal and mitochondrial systems. Defective mitophagy due to lysosomal malfunction may further exacerbate the mitochondrial dysfunction caused by Aβ. Future studies may clarify whether the highly significant negative correlation between down-regulation of nuclearencoded mitochondrial genes and Braak and CERAD scores means that mitochondrion stress in astrocytes contribute to the progression of AD.

CRediT authorship contribution statement

Declaration of Competing Interest

None.

Fig. 1 .

 1 Fig. 1. Workflow for the identification of altered astrocytic functions in AD and MCI (See Methods for details).Step 1-3: Genes from human-brain cortex were classified as cell-specific or enriched using RNAseq data from cells isolated from aged healthy human brain cortices (Step 1),

Fig. 2 .

 2 Fig. 2. Identification of cell-specific gene clusters in whole-brain transcriptomes by two-step hierarchical clustering. A. Hierarchical clustering was performed on transcriptomes of brain cell populations isolated from individual healthy subjects.Each column is a sample from a single subject and each row is a z-scored gene. The squares annotate the number and location of clusters with the highest z scores, that coincide with the highest τ scores for a given cell type (cluster information in Supplementary file 1, tables 6-8). The turquoise and purple squares are the main astrocytic and neuronal clusters (cluster # 4 and cluster # 13), respectively. Microglia, oligodendrocyte, and endothelial-cell clusters are 17, 164 and 5, respectively (shown in yellow). Clusters 30, 15, 10 and 116 are highly enriched in microglial or neuronal genes, but fall below the criteria established for cell-specificity (Methods). B. Whole-brain transcriptomes from individual subjects of Control (blue), MCI (beige) and AD (red) groups from MAYO, ROSMAP, MtSINAI were hierarchically clustered. The rows from the cell-type clustering in A. were kept fixed to identify cellspecific clusters. The figure was generated with raw data before correction for cell composition. The clustering segregated all groups with the same clinical diagnosis into at least 2 subgroups (defined as 1 and 2) showing distinct changes in gene-expression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .

 3 Fig. 3. Neuronal and astrocytic gene sub-clusters in the three clinical databases.A. Distributions of gene expression differences between AD2 and Control1 cohorts before (yellow) and after (blue) variance stabilizing normalization, yielding zerocentered distributions (Methods). B. Hierarchical clustering of Control (blue), MCI (beige), and AD (red) subjects from Fig.2Busing normalized data in rows. For both astrocytic and neuronal genes, subjects self-organized into at least two distinct subgroups (defined as 1 and 2) within each group with the same clinical diagnosis, while genes self-organized into up to four groups (Na-c and Aa-d). In MAYO and ROSMAP, clusters showed opposite directions of change in type 1 and 2 subjects. In MtSINAI, Control2 was like Control1, and AD1 like AD2, although with less pronounced changes. C. Violin plots of pathological stages according to CERAD and Braak scoring shows no significant difference in pathological stages between sub-groups within each cohort in ROSMAP (means ± SEM, one-way ANOVA with Tukey post-hoc). D. Distribution of APOE genotypes and sex in the six groups. While 3/3 was the predominant genotype across all groups, and 3/4 was more abundant in AD, there were no differences in frequency of APOE4 alleles between type 1 and 2 subgroups of subjects with the same clinical diagnosis. Likewise, females were more abundant than males in all groups, but not in type 2 subjects as compared to type 1 (p values determined by pairwise Fisher's Exact test with Bonferroni adjustment).

Fig. 4 .

 4 Fig. 4. GSVA unravels neuron-enriched functions altered in AD. Differences in the general functional categories among cohorts in MAYO, ROSMAP and MtSINAI databases were examined by GSVA. A. Comparison between AD2 and Control1. Functions are ranked according to the difference in enrichment scores (ES). Significantly changed functions (qFDR <0.05) are labeled in red. B. Enrichment scores of functional categories significantly changed in at least 2 out of 3 databases. (#) qFDR<0.05 (&) qFDR <0.01, and (*) qFDR <0.001, vs Control1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .

 5 Fig. 5. GSVA unravels astrocyte-enriched functions altered in AD. Differences in the functional categories among cohorts in MAYO, ROSMAP and MtSINAI databases were examined by GSVA. A. Comparison between AD2 and Control1. Functions are ranked according to the difference in enrichment scores (ES). Significantly changed functions (qFDR <0.05) are labeled in red. B. Enrichment scores of functional categories significantly changed in at least 2 out of 3 databases. (#) qFDR <0.05, (&) qFDR <0.01, and (*) qFDR <0.001, vs Control1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .

 6 Fig. 6. PCA identifies top altered functions in AD astrocytes. A. PCA analysis of astrocyte-enriched genes in ROSMAP. Dot plots and violin plots along PC1 (mean ± SEM) show differential expression of astrocyte genes in type 2 vs type 1 subjects. B. Top varying genes in PC1 (up-regulated and down-regulated), colour-coded for function according to the circos plot, inform about functions altered in AD astrocytes. The vertical axis represents the relative contribution (importance/weight) of each gene for the scoring C. of individual samples in panel A along PC1. C. Scoring of MtSINAI and MAYO databases using genes overlapping with ROSMAP PCA in A. discriminates Control and AD cohorts along PC1, as shown by the dot and violin plots (means ± SEM).

Fig. 7 .

 7 Fig. 7. Model of cortical human AD astrocytes and comparison with mice. A. The model revolves around two predictions.First prediction: Amyloid-β drives the phenotypical change of cortical astrocytes in AD by direct and indirect mechanisms. Direct: Amyloid-β causes endolysosomal and mitochondrial dysfunction in astrocytes, resulting in an adaptive astroprotective shift of gliotransmission to a GABA-predominant mode as a strategy to preserve mitochondrial respiration. Indirect: Synaptic dysfunction caused by amyloid-β in neurons augments the production of the PAP machinery and synaptogenic factors in astrocytes-a phenomenon of enhanced astrocytic plasticity to protect neurons. Second prediction: The phenotypical transformation of astrocytes in AD happens in early and late stages. Control2 subjects, cognitively normal, would be at stage a, and AD2 patients, diagnosed with dementia, at stage b. In stage a, the astroprotective change in glutamate/GABA fluxes may preserve neural-circuit homeostasis, because the GABA produced by MAO-B may locally compensate for the deficits in inhibitory tone that disrupts circuitry synchronization[START_REF] Lee | Slow wave sleep is a promising intervention target for Alzheimer's disease[END_REF]. By contrast, deficits in the endolysosomal system would result in partial phagocytosis of amyloid-β[START_REF] Sollvander | Accumulation of amyloid-beta by astrocytes result in enlarged endosomes and microvesicle-induced apoptosis of neurons[END_REF] and dystrophic neurites[START_REF] Gomez-Arboledas | Phagocytic clearance of presynaptic dystrophies by reactive astrocytes in Alzheimer's disease[END_REF], thus exacerbating amyloid-β accumulation, and hindering neuronal repair due to lack of debris elimination. In stage b, the mechanisms triggered to protect mitochondria may fail, and cortical astrocytes become deeply malfunctional, as suggested by the reduction in TCA cycle and ETC transcripts, despite the increase in protective stress responses and pro-survival pathways aimed at preventing astrocyte demise. For simplicity, only predominantly altered functions are depicted. ECM, extracellular matrix; ETC: electron transport chain, GLU: glutamate; oAβ: oligomeric amyloid-β; PAP: perisynaptic astrocyte processes; TCA: tricarboxylic acid. Figure adapted from Servier Medical Art (https://smart.servier.com/). B. GSVA of the transcriptome of astrocytes isolated from APP/PS1 mouse cortices using our manually curated functional categorization shows down-regulation of astrocyte-enriched functions. C. Left, GSVA of the transcriptome of APP/PS1 astrocytes using gene sets consisting of the top/bottom 50 dysregulated human astrocyte genes in PC1 from the ROSMAP database, and the top/bottom 50 DEGs from APP/PS1 astrocytes[START_REF] Orre | Isolation of glia from Alzheimer's mice reveals inflammation and dysfunction[END_REF]. Human astrocyte-enriched genes dysregulated in AD (up and down) are down-regulated in APP/PS1 astrocytes, confirming B. Right, GSVA of PC1 AD and mouse APP/PS1 DEG gene sets in ROSMAP, MAYO and MtSINAI whole-brain transcriptomes shows that both the up and down-regulated mouse gene pools tend to be up-regulated in AD databases. Discussion of the result of this comparison can be found in the main text.

Elena

  Galea: Conceptualization, Methodology, Data curation, Funding acquisition, Supervision, Writingoriginal draft, Writingreview & editing. Laura D. Weinstock: Methodology, Software, Data curation, Formal analysis, Writingreview & editing. Raquel Larramona-Arcas: Data curation, Writingreview & editing. Alyssa F. Pybus: Methodology, Software, Data curation, Formal analysis, Writing review & editing. Lydia Giménez-Llort: Writingreview & editing. Carole Escartin: Writingoriginal draft, Writingreview & editing. Levi B. Wood: Conceptualization, Methodology, Software, Data curation, Formal analysis, Funding acquisition, Supervision, Writingreview & editing.

  Type 1 sub-cohorts included Control1, MCI1 and AD1, while type 2 included Control2, MCI2 and AD2. Numbers of Con-trol1/Control2 patients were 33/45 in MAYO, 51/149 in ROSMAP, and 81/20 in MtSINAI. Numbers of AD1/AD2 patients were 28/54 in MAYO, 70/149 in ROSMAP, and 66/63 in MtSINAI. MCI1/MC2 numbers were 99/58 in ROSMAP. Heterogeneity in MtSINAI data has been recently reported and replicated in ROSMAP

Acknowledgements

The authors thank Dr. Alberto Lleó of the Unitat de Memòria at the Hospital Sant Pau, Barcelona, Spain, for critical reading of the manuscript. This work was supported by Fundació Marato-TV3 (grant number 201414-30), and Grup de Recerca en demèncias: Sant Pau SGR2017547, Generalitat de Catalunya (E.G.), by funds from The Rotary Coins for Alzheimer's Research Trust Fund (L.B.W.) and startup funds from the Woodruff School of Mechanical Engineering at Georgia Tech (L.B.W.). L. D.W. and A.F.P were supported in part by the National Institutes of Health Cell and Tissue Engineering Biotechnology Training Grant (T32-GM008433). CNRS, CEA and France Alzheimer (C.E.). The results published here are in whole or in part based on data obtained from the AD Knowledge Portal (https://adknowledgeportal.synapse.org). The authors thank the researchers in charge of the portal for their kind and patient assistance. Study data were provided by the following sources: The Mayo Clinic Alzheimers Disease Genetic Studies, led by Dr. Nilufer Ertekin-Taner and Dr. Steven G. Younkin, Mayo Clinic, Jacksonville, FL using samples from the Mayo Clinic Study of Aging, the Mayo Clinic Alzheimers Disease Research Center, and the Mayo Clinic Brain Bank. Data collection was supported through funding by NIA grants P50 AG016574, R01 AG032990, U01 AG046139, R01 AG018023, U01 AG006576, U01 AG006786, R01 AG025711, R01 AG017216, R01 AG003949, NINDS grant R01 NS080820, CurePSP Foundation, and support from Mayo Foundation. Study data includes samples collected through the Sun Health Research Institute Brain and Body Donation Program of Sun City, Arizona. The Brain and Body Donation Program is supported by the National Institute of Neurological Disorders and Stroke (U24 NS072026 National Brain and Tissue Resource for Parkinson's Disease and Related Disorders), the National Institute on Aging (P30 AG19610 Arizona Alzheimer's Disease Core Center), the Arizona Department of Health Services (contract 211002, Arizona Alzheimer's Research Center), the Arizona Biomedical Research Commission (contracts 4001, 0011, 05-901 and 1001 to the Arizona Parkinson's Disease Consortium) and the Michael J. Fox Foundation for Parkinson's

Research. Study data were also provided by the Rush Alzheimer's Disease Center, Rush University Medical Center, Chicago. Data collection was supported through funding by NIA grants P30AG10161 (ROS), R01AG15819 (ROSMAP; genomics and RNAseq), R01AG17917 (MAP),

R01AG30146, R01AG36042 (5hC methylation, ATACseq), RC2AG036547 (H3K9Ac), R01AG36836 (RNAseq), R01AG48015 (monocyte RNAseq) RF1AG57473 (single nucleus RNAseq), U01AG32984 (genomic and whole exome sequencing), U01AG46152 (ROSMAP AMP-AD, targeted proteomics), U01AG46161(TMT proteomics), U01AG61356 (whole genome sequencing, targeted proteomics, ROSMAP AMP-AD), the Illinois Department of Public Health (ROSMAP), and the Translational Genomics Research Institute (genomic). ROSMAP data can be requested at www.radc.rush.edu. All data generated during this study are included in this published article and its supplementary information files. Codes are available upon reasonable request.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi. org/10.1016/j.nbd.2022.105655.